Home | History | Annotate | Download | only in examples 
      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 (at) 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:
     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.
     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].
    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 
    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