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 
55 DEFINE_string(input, "", "File to which the output image should be written");
56 DEFINE_string(foe_file, "", "FoE file to use");
57 DEFINE_string(output, "", "File to which the output image should be written");
58 DEFINE_double(sigma, 20.0, "Standard deviation of noise");
59 DEFINE_bool(verbose, false, "Prints information about the solver progress.");
60 DEFINE_bool(line_search, false, "Use a line search instead of trust region "
61             "algorithm.");
62 
63 namespace ceres {
64 namespace examples {
65 
66 // This cost function is used to build the data term.
67 //
68 //   f_i(x) = a * (x_i - b)^2
69 //
70 class QuadraticCostFunction : public ceres::SizedCostFunction<1, 1> {
71  public:
QuadraticCostFunction(double a,double b)72   QuadraticCostFunction(double a, double b)
73     : sqrta_(std::sqrt(a)), b_(b) {}
Evaluate(double const * const * parameters,double * residuals,double ** jacobians) const74   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.
CreateProblem(const FieldsOfExperts & foe,const PGMImage<double> & image,Problem * problem,PGMImage<double> * solution)89 void 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].
SolveProblem(Problem * problem,PGMImage<double> * solution)139 void 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 
main(int argc,char ** argv)173 int 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