GoogleGit

blob: a8bdd7ceceda353e504c693a48d6b0a2238cd056 [file] [log] [blame]
  1. // Ceres Solver - A fast non-linear least squares minimizer
  2. // Copyright 2015 Google Inc. All rights reserved.
  3. // http://ceres-solver.org/
  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. #include <algorithm>
  41. #include <cmath>
  42. #include <iostream>
  43. #include <vector>
  44. #include <sstream>
  45. #include <string>
  46. #include "ceres/ceres.h"
  47. #include "gflags/gflags.h"
  48. #include "glog/logging.h"
  49. #include "fields_of_experts.h"
  50. #include "pgm_image.h"
  51. DEFINE_string(input, "", "File to which the output image should be written");
  52. DEFINE_string(foe_file, "", "FoE file to use");
  53. DEFINE_string(output, "", "File to which the output image should be written");
  54. DEFINE_double(sigma, 20.0, "Standard deviation of noise");
  55. DEFINE_bool(verbose, false, "Prints information about the solver progress.");
  56. DEFINE_bool(line_search, false, "Use a line search instead of trust region "
  57. "algorithm.");
  58. namespace ceres {
  59. namespace examples {
  60. // This cost function is used to build the data term.
  61. //
  62. // f_i(x) = a * (x_i - b)^2
  63. //
  64. class QuadraticCostFunction : public ceres::SizedCostFunction<1, 1> {
  65. public:
  66. QuadraticCostFunction(double a, double b)
  67. : sqrta_(std::sqrt(a)), b_(b) {}
  68. virtual bool Evaluate(double const* const* parameters,
  69. double* residuals,
  70. double** jacobians) const {
  71. const double x = parameters[0][0];
  72. residuals[0] = sqrta_ * (x - b_);
  73. if (jacobians != NULL && jacobians[0] != NULL) {
  74. jacobians[0][0] = sqrta_;
  75. }
  76. return true;
  77. }
  78. private:
  79. double sqrta_, b_;
  80. };
  81. // Creates a Fields of Experts MAP inference problem.
  82. void CreateProblem(const FieldsOfExperts& foe,
  83. const PGMImage<double>& image,
  84. Problem* problem,
  85. PGMImage<double>* solution) {
  86. // Create the data term
  87. CHECK_GT(FLAGS_sigma, 0.0);
  88. const double coefficient = 1 / (2.0 * FLAGS_sigma * FLAGS_sigma);
  89. for (unsigned index = 0; index < image.NumPixels(); ++index) {
  90. ceres::CostFunction* cost_function =
  91. new QuadraticCostFunction(coefficient,
  92. image.PixelFromLinearIndex(index));
  93. problem->AddResidualBlock(cost_function,
  94. NULL,
  95. solution->MutablePixelFromLinearIndex(index));
  96. }
  97. // Create Ceres cost and loss functions for regularization. One is needed for
  98. // each filter.
  99. std::vector<ceres::LossFunction*> loss_function(foe.NumFilters());
  100. std::vector<ceres::CostFunction*> cost_function(foe.NumFilters());
  101. for (int alpha_index = 0; alpha_index < foe.NumFilters(); ++alpha_index) {
  102. loss_function[alpha_index] = foe.NewLossFunction(alpha_index);
  103. cost_function[alpha_index] = foe.NewCostFunction(alpha_index);
  104. }
  105. // Add FoE regularization for each patch in the image.
  106. for (int x = 0; x < image.width() - (foe.Size() - 1); ++x) {
  107. for (int y = 0; y < image.height() - (foe.Size() - 1); ++y) {
  108. // Build a vector with the pixel indices of this patch.
  109. std::vector<double*> pixels;
  110. const std::vector<int>& x_delta_indices = foe.GetXDeltaIndices();
  111. const std::vector<int>& y_delta_indices = foe.GetYDeltaIndices();
  112. for (int i = 0; i < foe.NumVariables(); ++i) {
  113. double* pixel = solution->MutablePixel(x + x_delta_indices[i],
  114. y + y_delta_indices[i]);
  115. pixels.push_back(pixel);
  116. }
  117. // For this patch with coordinates (x, y), we will add foe.NumFilters()
  118. // terms to the objective function.
  119. for (int alpha_index = 0; alpha_index < foe.NumFilters(); ++alpha_index) {
  120. problem->AddResidualBlock(cost_function[alpha_index],
  121. loss_function[alpha_index],
  122. pixels);
  123. }
  124. }
  125. }
  126. }
  127. // Solves the FoE problem using Ceres and post-processes it to make sure the
  128. // solution stays within [0, 255].
  129. void SolveProblem(Problem* problem, PGMImage<double>* solution) {
  130. // These parameters may be experimented with. For example, ceres::DOGLEG tends
  131. // to be faster for 2x2 filters, but gives solutions with slightly higher
  132. // objective function value.
  133. ceres::Solver::Options options;
  134. options.max_num_iterations = 100;
  135. if (FLAGS_verbose) {
  136. options.minimizer_progress_to_stdout = true;
  137. }
  138. if (FLAGS_line_search) {
  139. options.minimizer_type = ceres::LINE_SEARCH;
  140. }
  141. options.linear_solver_type = ceres::SPARSE_NORMAL_CHOLESKY;
  142. options.function_tolerance = 1e-3; // Enough for denoising.
  143. ceres::Solver::Summary summary;
  144. ceres::Solve(options, problem, &summary);
  145. if (FLAGS_verbose) {
  146. std::cout << summary.FullReport() << "\n";
  147. }
  148. // Make the solution stay in [0, 255].
  149. for (int x = 0; x < solution->width(); ++x) {
  150. for (int y = 0; y < solution->height(); ++y) {
  151. *solution->MutablePixel(x, y) =
  152. std::min(255.0, std::max(0.0, solution->Pixel(x, y)));
  153. }
  154. }
  155. }
  156. } // namespace examples
  157. } // namespace ceres
  158. int main(int argc, char** argv) {
  159. using namespace ceres::examples;
  160. std::string
  161. usage("This program denoises an image using Ceres. Sample usage:\n");
  162. usage += argv[0];
  163. usage += " --input=<noisy image PGM file> --foe_file=<FoE file name>";
  164. CERES_GFLAGS_NAMESPACE::SetUsageMessage(usage);
  165. CERES_GFLAGS_NAMESPACE::ParseCommandLineFlags(&argc, &argv, true);
  166. google::InitGoogleLogging(argv[0]);
  167. if (FLAGS_input.empty()) {
  168. std::cerr << "Please provide an image file name.\n";
  169. return 1;
  170. }
  171. if (FLAGS_foe_file.empty()) {
  172. std::cerr << "Please provide a Fields of Experts file name.\n";
  173. return 1;
  174. }
  175. // Load the Fields of Experts filters from file.
  176. FieldsOfExperts foe;
  177. if (!foe.LoadFromFile(FLAGS_foe_file)) {
  178. std::cerr << "Loading \"" << FLAGS_foe_file << "\" failed.\n";
  179. return 2;
  180. }
  181. // Read the images
  182. PGMImage<double> image(FLAGS_input);
  183. if (image.width() == 0) {
  184. std::cerr << "Reading \"" << FLAGS_input << "\" failed.\n";
  185. return 3;
  186. }
  187. PGMImage<double> solution(image.width(), image.height());
  188. solution.Set(0.0);
  189. ceres::Problem problem;
  190. CreateProblem(foe, image, &problem, &solution);
  191. SolveProblem(&problem, &solution);
  192. if (!FLAGS_output.empty()) {
  193. CHECK(solution.WriteToFile(FLAGS_output))
  194. << "Writing \"" << FLAGS_output << "\" failed.";
  195. }
  196. return 0;
  197. }