1// Ceres Solver - A fast non-linear least squares minimizer
2// Copyright 2012 Google Inc. All rights reserved.
3// http://code.google.com/p/ceres-solver/
4//
5// Redistribution and use in source and binary forms, with or without
6// modification, are permitted provided that the following conditions are met:
7//
8// * Redistributions of source code must retain the above copyright notice,
9//   this list of conditions and the following disclaimer.
10// * Redistributions in binary form must reproduce the above copyright notice,
11//   this list of conditions and the following disclaimer in the documentation
12//   and/or other materials provided with the distribution.
13// * Neither the name of Google Inc. nor the names of its contributors may be
14//   used to endorse or promote products derived from this software without
15//   specific prior written permission.
16//
17// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
18// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
19// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
20// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
21// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
22// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
23// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
24// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
25// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
26// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
27// POSSIBILITY OF SUCH DAMAGE.
28//
29// Author: strandmark@google.com (Petter Strandmark)
30//
31// Denoising using Fields of Experts and the Ceres minimizer.
32//
33// Note that for good denoising results the weighting between the data term
34// and the Fields of Experts term needs to be adjusted. This is discussed
35// in [1]. This program assumes Gaussian noise. The noise model can be changed
36// by substituing another function for QuadraticCostFunction.
37//
38// [1] S. Roth and M.J. Black. "Fields of Experts." International Journal of
39//     Computer Vision, 82(2):205--229, 2009.
40
41#include <algorithm>
42#include <cmath>
43#include <iostream>
44#include <vector>
45#include <sstream>
46#include <string>
47
48#include "ceres/ceres.h"
49#include "gflags/gflags.h"
50#include "glog/logging.h"
51
52#include "fields_of_experts.h"
53#include "pgm_image.h"
54
55DEFINE_string(input, "", "File to which the output image should be written");
56DEFINE_string(foe_file, "", "FoE file to use");
57DEFINE_string(output, "", "File to which the output image should be written");
58DEFINE_double(sigma, 20.0, "Standard deviation of noise");
59DEFINE_bool(verbose, false, "Prints information about the solver progress.");
60DEFINE_bool(line_search, false, "Use a line search instead of trust region "
61            "algorithm.");
62
63namespace ceres {
64namespace examples {
65
66// This cost function is used to build the data term.
67//
68//   f_i(x) = a * (x_i - b)^2
69//
70class QuadraticCostFunction : public ceres::SizedCostFunction<1, 1> {
71 public:
72  QuadraticCostFunction(double a, double b)
73    : sqrta_(std::sqrt(a)), b_(b) {}
74  virtual bool Evaluate(double const* const* parameters,
75                        double* residuals,
76                        double** jacobians) const {
77    const double x = parameters[0][0];
78    residuals[0] = sqrta_ * (x - b_);
79    if (jacobians != NULL && jacobians[0] != NULL) {
80      jacobians[0][0] = sqrta_;
81    }
82    return true;
83  }
84 private:
85  double sqrta_, b_;
86};
87
88// Creates a Fields of Experts MAP inference problem.
89void CreateProblem(const FieldsOfExperts& foe,
90                   const PGMImage<double>& image,
91                   Problem* problem,
92                   PGMImage<double>* solution) {
93  // Create the data term
94  CHECK_GT(FLAGS_sigma, 0.0);
95  const double coefficient = 1 / (2.0 * FLAGS_sigma * FLAGS_sigma);
96  for (unsigned index = 0; index < image.NumPixels(); ++index) {
97    ceres::CostFunction* cost_function =
98        new QuadraticCostFunction(coefficient,
99                                  image.PixelFromLinearIndex(index));
100    problem->AddResidualBlock(cost_function,
101                              NULL,
102                              solution->MutablePixelFromLinearIndex(index));
103  }
104
105  // Create Ceres cost and loss functions for regularization. One is needed for
106  // each filter.
107  std::vector<ceres::LossFunction*> loss_function(foe.NumFilters());
108  std::vector<ceres::CostFunction*> cost_function(foe.NumFilters());
109  for (int alpha_index = 0; alpha_index < foe.NumFilters(); ++alpha_index) {
110    loss_function[alpha_index] = foe.NewLossFunction(alpha_index);
111    cost_function[alpha_index] = foe.NewCostFunction(alpha_index);
112  }
113
114  // Add FoE regularization for each patch in the image.
115  for (int x = 0; x < image.width() - (foe.Size() - 1); ++x) {
116    for (int y = 0; y < image.height() - (foe.Size() - 1); ++y) {
117      // Build a vector with the pixel indices of this patch.
118      std::vector<double*> pixels;
119      const std::vector<int>& x_delta_indices = foe.GetXDeltaIndices();
120      const std::vector<int>& y_delta_indices = foe.GetYDeltaIndices();
121      for (int i = 0; i < foe.NumVariables(); ++i) {
122        double* pixel = solution->MutablePixel(x + x_delta_indices[i],
123                                               y + y_delta_indices[i]);
124        pixels.push_back(pixel);
125      }
126      // For this patch with coordinates (x, y), we will add foe.NumFilters()
127      // terms to the objective function.
128      for (int alpha_index = 0; alpha_index < foe.NumFilters(); ++alpha_index) {
129        problem->AddResidualBlock(cost_function[alpha_index],
130                                  loss_function[alpha_index],
131                                  pixels);
132      }
133    }
134  }
135}
136
137// Solves the FoE problem using Ceres and post-processes it to make sure the
138// solution stays within [0, 255].
139void SolveProblem(Problem* problem, PGMImage<double>* solution) {
140  // These parameters may be experimented with. For example, ceres::DOGLEG tends
141  // to be faster for 2x2 filters, but gives solutions with slightly higher
142  // objective function value.
143  ceres::Solver::Options options;
144  options.max_num_iterations = 100;
145  if (FLAGS_verbose) {
146    options.minimizer_progress_to_stdout = true;
147  }
148
149  if (FLAGS_line_search) {
150    options.minimizer_type = ceres::LINE_SEARCH;
151  }
152
153  options.linear_solver_type = ceres::SPARSE_NORMAL_CHOLESKY;
154  options.function_tolerance = 1e-3;  // Enough for denoising.
155
156  ceres::Solver::Summary summary;
157  ceres::Solve(options, problem, &summary);
158  if (FLAGS_verbose) {
159    std::cout << summary.FullReport() << "\n";
160  }
161
162  // Make the solution stay in [0, 255].
163  for (int x = 0; x < solution->width(); ++x) {
164    for (int y = 0; y < solution->height(); ++y) {
165      *solution->MutablePixel(x, y) =
166          std::min(255.0, std::max(0.0, solution->Pixel(x, y)));
167    }
168  }
169}
170}  // namespace examples
171}  // namespace ceres
172
173int main(int argc, char** argv) {
174  using namespace ceres::examples;
175  std::string
176      usage("This program denoises an image using Ceres.  Sample usage:\n");
177  usage += argv[0];
178  usage += " --input=<noisy image PGM file> --foe_file=<FoE file name>";
179  google::SetUsageMessage(usage);
180  google::ParseCommandLineFlags(&argc, &argv, true);
181  google::InitGoogleLogging(argv[0]);
182
183  if (FLAGS_input.empty()) {
184    std::cerr << "Please provide an image file name.\n";
185    return 1;
186  }
187
188  if (FLAGS_foe_file.empty()) {
189    std::cerr << "Please provide a Fields of Experts file name.\n";
190    return 1;
191  }
192
193  // Load the Fields of Experts filters from file.
194  FieldsOfExperts foe;
195  if (!foe.LoadFromFile(FLAGS_foe_file)) {
196    std::cerr << "Loading \"" << FLAGS_foe_file << "\" failed.\n";
197    return 2;
198  }
199
200  // Read the images
201  PGMImage<double> image(FLAGS_input);
202  if (image.width() == 0) {
203    std::cerr << "Reading \"" << FLAGS_input << "\" failed.\n";
204    return 3;
205  }
206  PGMImage<double> solution(image.width(), image.height());
207  solution.Set(0.0);
208
209  ceres::Problem problem;
210  CreateProblem(foe, image, &problem, &solution);
211
212  SolveProblem(&problem, &solution);
213
214  if (!FLAGS_output.empty()) {
215    CHECK(solution.WriteToFile(FLAGS_output))
216        << "Writing \"" << FLAGS_output << "\" failed.";
217  }
218
219  return 0;
220}
221