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 61 namespace ceres { 62 namespace examples { 63 64 // This cost function is used to build the data term. 65 // 66 // f_i(x) = a * (x_i - b)^2 67 // 68 class QuadraticCostFunction : public ceres::SizedCostFunction<1, 1> { 69 public: 70 QuadraticCostFunction(double a, double b) 71 : sqrta_(std::sqrt(a)), b_(b) {} 72 virtual bool Evaluate(double const* const* parameters, 73 double* residuals, 74 double** jacobians) const { 75 const double x = parameters[0][0]; 76 residuals[0] = sqrta_ * (x - b_); 77 if (jacobians != NULL && jacobians[0] != NULL) { 78 jacobians[0][0] = sqrta_; 79 } 80 return true; 81 } 82 private: 83 double sqrta_, b_; 84 }; 85 86 // Creates a Fields of Experts MAP inference problem. 87 void CreateProblem(const FieldsOfExperts& foe, 88 const PGMImage<double>& image, 89 Problem* problem, 90 PGMImage<double>* solution) { 91 // Create the data term 92 CHECK_GT(FLAGS_sigma, 0.0); 93 const double coefficient = 1 / (2.0 * FLAGS_sigma * FLAGS_sigma); 94 for (unsigned index = 0; index < image.NumPixels(); ++index) { 95 ceres::CostFunction* cost_function = 96 new QuadraticCostFunction(coefficient, 97 image.PixelFromLinearIndex(index)); 98 problem->AddResidualBlock(cost_function, 99 NULL, 100 solution->MutablePixelFromLinearIndex(index)); 101 } 102 103 // Create Ceres cost and loss functions for regularization. One is needed for 104 // each filter. 105 std::vector<ceres::LossFunction*> loss_function(foe.NumFilters()); 106 std::vector<ceres::CostFunction*> cost_function(foe.NumFilters()); 107 for (int alpha_index = 0; alpha_index < foe.NumFilters(); ++alpha_index) { 108 loss_function[alpha_index] = foe.NewLossFunction(alpha_index); 109 cost_function[alpha_index] = foe.NewCostFunction(alpha_index); 110 } 111 112 // Add FoE regularization for each patch in the image. 113 for (int x = 0; x < image.width() - (foe.Size() - 1); ++x) { 114 for (int y = 0; y < image.height() - (foe.Size() - 1); ++y) { 115 // Build a vector with the pixel indices of this patch. 116 std::vector<double*> pixels; 117 const std::vector<int>& x_delta_indices = foe.GetXDeltaIndices(); 118 const std::vector<int>& y_delta_indices = foe.GetYDeltaIndices(); 119 for (int i = 0; i < foe.NumVariables(); ++i) { 120 double* pixel = solution->MutablePixel(x + x_delta_indices[i], 121 y + y_delta_indices[i]); 122 pixels.push_back(pixel); 123 } 124 // For this patch with coordinates (x, y), we will add foe.NumFilters() 125 // terms to the objective function. 126 for (int alpha_index = 0; alpha_index < foe.NumFilters(); ++alpha_index) { 127 problem->AddResidualBlock(cost_function[alpha_index], 128 loss_function[alpha_index], 129 pixels); 130 } 131 } 132 } 133 } 134 135 // Solves the FoE problem using Ceres and post-processes it to make sure the 136 // solution stays within [0, 255]. 137 void SolveProblem(Problem* problem, PGMImage<double>* solution) { 138 // These parameters may be experimented with. For example, ceres::DOGLEG tends 139 // to be faster for 2x2 filters, but gives solutions with slightly higher 140 // objective function value. 141 ceres::Solver::Options options; 142 options.max_num_iterations = 100; 143 if (FLAGS_verbose) { 144 options.minimizer_progress_to_stdout = true; 145 } 146 options.trust_region_strategy_type = ceres::LEVENBERG_MARQUARDT; 147 options.linear_solver_type = ceres::SPARSE_NORMAL_CHOLESKY; 148 options.function_tolerance = 1e-3; // Enough for denoising. 149 150 ceres::Solver::Summary summary; 151 ceres::Solve(options, problem, &summary); 152 if (FLAGS_verbose) { 153 std::cout << summary.FullReport() << "\n"; 154 } 155 156 // Make the solution stay in [0, 255]. 157 for (int x = 0; x < solution->width(); ++x) { 158 for (int y = 0; y < solution->height(); ++y) { 159 *solution->MutablePixel(x, y) = 160 std::min(255.0, std::max(0.0, solution->Pixel(x, y))); 161 } 162 } 163 } 164 } // namespace examples 165 } // namespace ceres 166 167 int main(int argc, char** argv) { 168 using namespace ceres::examples; 169 std::string 170 usage("This program denoises an image using Ceres. Sample usage:\n"); 171 usage += argv[0]; 172 usage += " --input=<noisy image PGM file> --foe_file=<FoE file name>"; 173 google::SetUsageMessage(usage); 174 google::ParseCommandLineFlags(&argc, &argv, true); 175 google::InitGoogleLogging(argv[0]); 176 177 if (FLAGS_input.empty()) { 178 std::cerr << "Please provide an image file name.\n"; 179 return 1; 180 } 181 182 if (FLAGS_foe_file.empty()) { 183 std::cerr << "Please provide a Fields of Experts file name.\n"; 184 return 1; 185 } 186 187 // Load the Fields of Experts filters from file. 188 FieldsOfExperts foe; 189 if (!foe.LoadFromFile(FLAGS_foe_file)) { 190 std::cerr << "Loading \"" << FLAGS_foe_file << "\" failed.\n"; 191 return 2; 192 } 193 194 // Read the images 195 PGMImage<double> image(FLAGS_input); 196 if (image.width() == 0) { 197 std::cerr << "Reading \"" << FLAGS_input << "\" failed.\n"; 198 return 3; 199 } 200 PGMImage<double> solution(image.width(), image.height()); 201 solution.Set(0.0); 202 203 ceres::Problem problem; 204 CreateProblem(foe, image, &problem, &solution); 205 206 SolveProblem(&problem, &solution); 207 208 if (!FLAGS_output.empty()) { 209 CHECK(solution.WriteToFile(FLAGS_output)) 210 << "Writing \"" << FLAGS_output << "\" failed."; 211 } 212 213 return 0; 214 } 215