Adding a denoising example using Fields of Experts.

We have permission from Stefan Roth to use the coefficients from his
Matlab toolbox. They have been added as *.foe files.

Change-Id: Ice529e5cab0302b9f27648dd3c8e5ed7b9662aba

data/2x2.foe

+2 3
+0 1 0 1
+0 0 1 1
+0.586612685392731 1.157638405566669 0.846059486257292
+-0.0582774013402734 0.0339010363051084 -0.0501593018104054 0.0745568557931712
+0.0492112815304123 -0.0307820846538285 -0.123247230948424 0.104812330861557
+0.0562633568728865 0.0152832583489560 -0.0576215592718086 -0.0139673758425540
+

data/3x3.foe

+3 8
+0 1 2 0 1 2 0 1 2
+0 0 0 1 1 1 2 2 2
+1.201143e-01 7.520515e-02 9.078330e-02 1.280545e-01 6.276734e-02 1.201840e-01 1.092460e-01 1.217102e-01
+-1.06290 1.81622 -0.80592 0.27489 0.13790 -0.49836 0.76669 -1.95329 1.32468
+-1.39695 1.62887 -0.27482 1.46226 -2.21903 0.82016 -0.12290 0.67906 -0.57683
+0.89899 0.27696 0.04805 -0.07191 0.00777 -0.00488 -0.80267 -0.31249 -0.03957
+-0.86612 1.66525 -1.03168 1.87671 -3.28802 1.92990 -1.22515 2.01814 -1.08079
+-0.37343 -0.23977 0.56955 -0.41741 0.05711 0.42774 -0.14881 0.02528 0.10072
+0.62783 0.27746 0.44180 -1.22190 -0.27753 -1.23532 0.54919 0.15855 0.68658
+0.46483 -1.12896 0.57926 0.22646 -0.25658 0.18881 0.57866 -1.13277 0.48707
+0.86662 0.19780 -1.06759 -1.92170 -0.10882 2.02231 1.14262 -0.09817 -1.03353
+

data/5x5.foe

+5 24
+0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4
+0 0 0 0 0 1 1 1 1 1 2 2 2 2 2 3 3 3 3 3 4 4 4 4 4
+
+5.472286e-02 6.425581e-02 5.031088e-02 5.989018e-02 4.424956e-02 5.119065e-02 6.743059e-02 2.339211e-02 5.973362e-02 6.281576e-02 1.495131e-02 5.693016e-02 3.189429e-02 7.672358e-02 4.524088e-02 3.681065e-02 6.383650e-02 4.347603e-02 2.264134e-02 2.561474e-02 5.886950e-02 4.114462e-02 4.893348e-02 5.834900e-02
+
+   -0.31071 1.04120 -0.78749 -0.47895 0.56812 -0.02713 0.12199 -0.19871 0.01017 0.08976 -0.05437 0.19638 -0.11093 -0.02853 0.00427 -0.19578 0.29112 -0.05108 -0.08012 0.02944 -0.87865 1.54098 -0.60090 -0.14424 0.05388 
+   0.59362 0.07378 0.01163 0.11076 0.38838 -1.56819 -0.12448 -0.26775 -0.14495 -1.26480 1.69550 0.23647 0.50640 0.20682 1.71141 -0.93659 -0.27274 -0.40550 -0.30937 -1.27168 0.21164 0.09704 0.14615 0.14513 0.43322 
+   0.84855 -1.15931 -0.04677 1.55189 -1.20232 -0.84787 0.88787 0.05491 -0.59425 0.44273 0.10771 0.15608 -0.20982 -0.40088 0.46844 0.00920 -0.34651 0.07722 1.10449 -1.00405 -0.00980 0.27614 0.32393 -1.83075 1.34408 
+   -1.03323 1.88111 -1.43100 0.78609 -0.25928 -0.03655 0.21945 -0.13154 -0.10918 0.11203 0.09740 -0.02567 -0.23481 0.18294 -0.07056 0.20877 -0.12025 -0.27859 0.34269 -0.10385 0.70720 -1.63594 1.39715 -0.56761 0.10165 
+   0.53443 -0.08024 -0.11745 -0.03402 0.25033 -0.96043 0.15798 0.15981 0.10961 -0.24634 0.77044 -0.18866 -0.27562 0.18315 0.07681 -0.74267 0.51002 0.40930 -0.00294 -1.63639 0.40120 -0.42511 -0.15842 -0.28824 1.59764 
+   -0.88573 -0.33458 -0.01419 -0.01103 0.00077 0.47044 0.43260 0.04198 -0.00840 0.00150 0.57390 0.08128 0.04955 -0.00643 0.01470 0.50375 0.07196 0.02376 0.03349 0.03517 -0.67904 -0.28139 -0.05936 -0.03884 -0.02249 
+   -0.01311 -0.01988 0.10721 -0.08376 0.01636 -0.32442 0.98192 -1.44473 1.21759 -0.47859 0.89529 -2.36566 3.17358 -2.70415 1.10629 -0.96782 2.47203 -3.22965 2.76881 -1.14343 0.42498 -1.09401 1.42223 -1.22142 0.50653 
+   0.21862 -0.20464 0.03462 -0.03753 0.01527 0.13402 -0.07069 -0.00552 -0.04434 -0.02653 0.15630 0.02185 0.04572 -0.11760 -0.06958 0.15633 -0.00898 0.00952 -0.06227 -0.12147 0.37032 0.05372 0.10450 -0.22764 -0.32445 
+   -0.02092 -0.35567 0.54827 0.50360 -0.95729 0.42842 0.61450 -1.06896 -0.96641 1.79327 -1.49830 -0.13333 0.96113 0.66182 -1.14234 2.35571 -0.68686 -0.72675 -0.09191 0.23092 -1.29917 0.56449 0.33405 -0.11554 0.06663 
+   -0.15212 -0.31661 -0.10629 -0.68078 1.64419 0.60883 0.48936 0.35629 0.85845 -2.71321 -1.02463 -0.40494 -0.33908 -0.03011 1.32162 0.87303 0.38931 0.15034 -0.46469 -0.16997 -0.30057 -0.21769 -0.03504 0.29551 -0.03051 
+   -0.63745 1.09815 0.07905 -1.29684 0.78552 0.31671 -0.60369 -0.31422 1.37878 -0.84617 -0.08181 0.25602 0.09426 -0.69533 0.50771 -0.26215 0.08620 0.18111 0.24247 -0.35641 0.18171 0.29815 -0.90155 0.49852 -0.00995 
+   -0.24394 0.32079 -0.18885 -0.43086 0.61484 0.95412 -1.08795 0.41381 1.25966 -1.62800 -1.63362 1.77016 -0.41047 -1.74982 2.02167 1.69856 -1.82451 0.27727 1.56121 -1.66173 -0.80206 0.85485 -0.10948 -0.64836 0.67266 
+   0.60487 -0.01646 0.13186 0.13391 0.66691 -0.16619 0.00981 0.02774 -0.14688 -0.23655 -0.03876 0.01428 -0.05848 0.00855 0.01915 0.28643 0.01215 0.00615 0.02837 -0.27025 -0.67548 -0.06341 -0.09492 -0.01467 -0.16753 
+   0.77534 -1.90839 2.47496 -2.26067 1.01011 -0.36223 0.89421 -1.24538 1.23108 -0.59609 -0.44538 1.04742 -1.25543 1.08209 -0.44929 -0.10218 0.51811 -0.92194 0.90930 -0.42024 0.18932 -0.73810 1.25011 -1.24150 0.56562 
+   0.93435 -0.03480 0.04710 -0.40291 -0.92490 -0.66962 0.13901 0.00497 0.31148 0.63533 0.30957 -0.23678 -0.07453 -0.24971 0.25021 0.68393 0.03453 -0.00104 -0.20619 -0.21134 -1.23377 0.08326 0.02804 0.53578 0.24696 
+   -0.31333 -0.06452 -0.14492 0.32917 0.18505 0.37680 -0.05739 -0.20598 -0.39758 0.33795 0.15776 -0.10590 -0.08163 -0.26879 0.26603 0.13942 0.14388 0.15602 -0.85931 0.46755 -0.37693 0.14924 0.18844 1.29502 -1.31627 
+   -0.57261 -0.15101 -0.03142 0.04828 -0.09470 1.30126 0.29676 0.02642 -0.10865 0.12241 -1.04666 -0.17702 0.10930 0.48806 0.79384 0.35790 0.00519 -0.09722 -0.76968 -1.93132 -0.05982 0.03488 -0.03692 0.34823 1.14489 
+   0.02861 -0.04563 -0.39979 -0.34141 0.84996 -0.11996 0.09395 0.10138 0.81867 -0.96541 0.17159 -0.18584 0.25071 -0.22578 -0.06749 1.50317 -1.28229 -0.16426 -0.02308 0.15558 -1.66312 1.45438 0.22480 -0.25807 0.09117 
+   -0.55490 0.19579 0.22691 0.60510 -0.46819 0.06905 -0.09310 -0.00511 -0.26410 0.29786 -0.15810 0.19022 0.02614 -0.00490 -0.06838 0.11087 -0.14454 -0.04154 -0.16050 0.23955 -0.60009 0.80946 0.13003 0.34930 -0.69227 
+   0.02280 0.50837 0.09286 -0.70868 0.16414 -0.80301 0.17835 0.36944 0.75854 -0.61169 0.44842 -0.99949 0.04847 0.08457 0.35292 0.60861 -0.02401 -0.48740 0.41152 -0.41065 -0.34356 0.41056 -0.01730 -0.58536 0.53149 
+   1.54887 -0.44899 -0.27397 -0.36683 -0.36138 -1.38762 0.36466 0.10235 0.35630 0.40105 -0.17410 -0.05898 0.13052 -0.04032 0.38137 -0.82904 0.31802 -0.02420 -0.20938 0.57263 0.95569 -0.26814 0.08144 0.25318 -1.02387 
+   0.26344 -0.49099 0.61113 -1.33142 0.89974 0.05990 -0.07519 0.14265 -0.33453 0.30248 0.08324 -0.12872 0.07988 -0.24786 0.14599 0.17826 -0.18072 0.08816 -0.10930 0.12304 0.45104 -0.70917 0.41756 -0.41199 0.18098 
+   -0.46895 0.79795 -0.38564 0.04986 -0.04034 0.11811 -0.23006 0.26011 -0.13551 0.03185 -0.01839 -0.09632 0.18000 -0.13208 0.03061 0.11165 -0.25590 0.42916 -0.49657 0.24341 0.51922 -1.52509 2.17579 -1.86952 0.70906 
+   -1.59203 1.37214 -0.36807 0.42638 -0.00689 2.86095 -2.02733 0.40858 -0.61961 -0.28165 -1.77476 0.55742 -0.12732 0.31007 0.76569 0.19805 0.83530 -0.18053 0.29590 -1.07493 0.27739 -0.69004 0.21195 -0.37112 0.59251 
+
+

data/README.foe

+The *.foe files contain coefficients provided by Stefan Roth, who agreed to
+release them under a BSD license. See his home page:
+
+  http://www.gris.informatik.tu-darmstadt.de/~sroth/research/foe/index.html
+
+The coefficients in the *.foe files have been obtained by extracting the
+matrices from the MATLAB files and performing matrix multiplication.
+
+The format of the files is ASCII:
+  <s = filter size> <K = number of filters>
+  <alpha_1> ... <alpha_K>
+  <f_1,1> ... <f_1,s^2>
+  ...
+  <f_K,1> ... <f_K,s^2>

examples/CMakeLists.txt

@@ -54,6 +54,11 @@ IF (${GFLAGS})
                  bundle_adjuster.cc
                  bal_problem.cc)
   TARGET_LINK_LIBRARIES(bundle_adjuster ceres)
+
+  ADD_EXECUTABLE(denoising
+                 denoising.cc
+                 fields_of_experts.cc)
+  TARGET_LINK_LIBRARIES(denoising ceres)
 ENDIF (${GFLAGS})
 
 ADD_EXECUTABLE(simple_bundle_adjuster

examples/denoising.cc

+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2012 Google Inc. All rights reserved.
+// http://code.google.com/p/ceres-solver/
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions are met:
+//
+// * Redistributions of source code must retain the above copyright notice,
+//   this list of conditions and the following disclaimer.
+// * Redistributions in binary form must reproduce the above copyright notice,
+//   this list of conditions and the following disclaimer in the documentation
+//   and/or other materials provided with the distribution.
+// * Neither the name of Google Inc. nor the names of its contributors may be
+//   used to endorse or promote products derived from this software without
+//   specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+// POSSIBILITY OF SUCH DAMAGE.
+//
+// Author: strandmark@google.com (Petter Strandmark)
+//
+// Denoising using Fields of Experts and the Ceres minimizer.
+//
+// Note that for good denoising results the weighting between the data term
+// and the Fields of Experts term needs to be adjusted. This is discussed
+// in [1]. This program assumes Gaussian noise. The noise model can be changed
+// by substituing another function for QuadraticCostFunction.
+//
+// [1] S. Roth and M.J. Black. "Fields of Experts." International Journal of
+//     Computer Vision, 82(2):205--229, 2009.
+
+#include <algorithm>
+#include <cmath>
+#include <iostream>
+#include <vector>
+#include <sstream>
+#include <string>
+
+#include "ceres/ceres.h"
+#include "gflags/gflags.h"
+#include "glog/logging.h"
+
+#include "fields_of_experts.h"
+#include "pgm_image.h"
+
+DEFINE_string(input, "", "File to which the output image should be written");
+DEFINE_string(foe_file, "", "FoE file to use");
+DEFINE_string(output, "", "File to which the output image should be written");
+DEFINE_double(sigma, 20.0, "Standard deviation of noise");
+DEFINE_bool(verbose, false, "Prints information about the solver progress.");
+
+namespace ceres {
+namespace examples {
+
+// This cost function is used to build the data term.
+//
+//   f_i(x) = a * (x_i - b)^2
+//
+class QuadraticCostFunction : public ceres::SizedCostFunction<1, 1> {
+ public:
+  QuadraticCostFunction(double a, double b)
+    : sqrta_(std::sqrt(a)), b_(b) {}
+  virtual bool Evaluate(double const* const* parameters,
+                        double* residuals,
+                        double** jacobians) const {
+    const double x = parameters[0][0];
+    residuals[0] = sqrta_ * (x - b_);
+    if (jacobians != NULL && jacobians[0] != NULL) {
+      jacobians[0][0] = sqrta_;
+    }
+    return true;
+  }
+ private:
+  double sqrta_, b_;
+};
+
+// Creates a Fields of Experts MAP inference problem.
+void CreateProblem(const FieldsOfExperts& foe,
+                   const PGMImage<double>& image,
+                   Problem* problem,
+                   PGMImage<double>* solution) {
+  // Create the data term
+  CHECK_GT(FLAGS_sigma, 0.0);
+  const double coefficient = 1 / (2.0 * FLAGS_sigma * FLAGS_sigma);
+  for (unsigned index = 0; index < image.NumPixels(); ++index) {
+    ceres::CostFunction* cost_function =
+        new QuadraticCostFunction(coefficient,
+                                  image.PixelFromLinearIndex(index));
+    problem->AddResidualBlock(cost_function,
+                              NULL,
+                              solution->MutablePixelFromLinearIndex(index));
+  }
+
+  // Create Ceres cost and loss functions for regularization. One is needed for
+  // each filter.
+  std::vector<ceres::LossFunction*> loss_function(foe.NumFilters());
+  std::vector<ceres::CostFunction*> cost_function(foe.NumFilters());
+  for (int alpha_index = 0; alpha_index < foe.NumFilters(); ++alpha_index) {
+    loss_function[alpha_index] = foe.NewLossFunction(alpha_index);
+    cost_function[alpha_index] = foe.NewCostFunction(alpha_index);
+  }
+
+  // Add FoE regularization for each patch in the image.
+  for (int x = 0; x < image.width() - (foe.Size() - 1); ++x) {
+    for (int y = 0; y < image.height() - (foe.Size() - 1); ++y) {
+      // Build a vector with the pixel indices of this patch.
+      std::vector<double*> pixels;
+      const std::vector<int>& x_delta_indices = foe.GetXDeltaIndices();
+      const std::vector<int>& y_delta_indices = foe.GetYDeltaIndices();
+      for (int i = 0; i < foe.NumVariables(); ++i) {
+        double* pixel = solution->MutablePixel(x + x_delta_indices[i],
+                                               y + y_delta_indices[i]);
+        pixels.push_back(pixel);
+      }
+      // For this patch with coordinates (x, y), we will add foe.NumFilters()
+      // terms to the objective function.
+      for (int alpha_index = 0; alpha_index < foe.NumFilters(); ++alpha_index) {
+        problem->AddResidualBlock(cost_function[alpha_index],
+                                  loss_function[alpha_index],
+                                  pixels);
+      }
+    }
+  }
+}
+
+// Solves the FoE problem using Ceres and post-processes it to make sure the
+// solution stays within [0, 255].
+void SolveProblem(Problem* problem, PGMImage<double>* solution) {
+  // These parameters may be experimented with. For example, ceres::DOGLEG tends
+  // to be faster for 2x2 filters, but gives solutions with slightly higher
+  // objective function value.
+  ceres::Solver::Options options;
+  options.max_num_iterations = 100;
+  if (FLAGS_verbose) {
+    options.minimizer_progress_to_stdout = true;
+  }
+  options.trust_region_strategy_type = ceres::LEVENBERG_MARQUARDT;
+  options.linear_solver_type = ceres::SPARSE_NORMAL_CHOLESKY;
+  options.function_tolerance = 1e-3;  // Enough for denoising.
+
+  ceres::Solver::Summary summary;
+  ceres::Solve(options, problem, &summary);
+  if (FLAGS_verbose) {
+    std::cout << summary.FullReport() << "\n";
+  }
+
+  // Make the solution stay in [0, 255].
+  for (int x = 0; x < solution->width(); ++x) {
+    for (int y = 0; y < solution->height(); ++y) {
+      *solution->MutablePixel(x, y) =
+          std::min(255.0, std::max(0.0, solution->Pixel(x, y)));
+    }
+  }
+}
+}  // namespace examples
+}  // namespace ceres
+
+int main(int argc, char** argv) {
+  using namespace ceres::examples;
+  std::string
+      usage("This program denoises an image using Ceres.  Sample usage:\n");
+  usage += argv[0];
+  usage += " --input=<noisy image PGM file> --foe_file=<FoE file name>";
+  google::SetUsageMessage(usage);
+  google::ParseCommandLineFlags(&argc, &argv, true);
+  google::InitGoogleLogging(argv[0]);
+
+  if (FLAGS_input.empty()) {
+    std::cerr << "Please provide an image file name.\n";
+    return 1;
+  }
+
+  if (FLAGS_foe_file.empty()) {
+    std::cerr << "Please provide a Fields of Experts file name.\n";
+    return 1;
+  }
+
+  // Load the Fields of Experts filters from file.
+  FieldsOfExperts foe;
+  if (!foe.LoadFromFile(FLAGS_foe_file)) {
+    std::cerr << "Loading \"" << FLAGS_foe_file << "\" failed.\n";
+    return 2;
+  }
+
+  // Read the images
+  PGMImage<double> image(FLAGS_input);
+  if (image.width() == 0) {
+    std::cerr << "Reading \"" << FLAGS_input << "\" failed.\n";
+    return 3;
+  }
+  PGMImage<double> solution(image.width(), image.height());
+  solution.Set(0.0);
+
+  ceres::Problem problem;
+  CreateProblem(foe, image, &problem, &solution);
+
+  SolveProblem(&problem, &solution);
+
+  if (!FLAGS_output.empty()) {
+    CHECK(solution.WriteToFile(FLAGS_output))
+        << "Writing \"" << FLAGS_output << "\" failed.";
+  }
+
+  return 0;
+}

examples/fields_of_experts.cc

+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2012 Google Inc. All rights reserved.
+// http://code.google.com/p/ceres-solver/
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions are met:
+//
+// * Redistributions of source code must retain the above copyright notice,
+//   this list of conditions and the following disclaimer.
+// * Redistributions in binary form must reproduce the above copyright notice,
+//   this list of conditions and the following disclaimer in the documentation
+//   and/or other materials provided with the distribution.
+// * Neither the name of Google Inc. nor the names of its contributors may be
+//   used to endorse or promote products derived from this software without
+//   specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+// POSSIBILITY OF SUCH DAMAGE.
+//
+// Author: strandmark@google.com (Petter Strandmark)
+//
+// Class for loading the data required for descibing a Fields of Experts (FoE)
+// model.
+
+#include "fields_of_experts.h"
+
+#include <fstream>
+#include <cmath>
+
+#include "pgm_image.h"
+
+namespace ceres {
+namespace examples {
+
+FieldsOfExpertsCost::FieldsOfExpertsCost(const std::vector<double>& filter)
+    : filter_(filter) {
+  set_num_residuals(1);
+  for (int i = 0; i < filter_.size(); ++i) {
+    mutable_parameter_block_sizes()->push_back(1);
+  }
+}
+
+// This is a dot product between a the scalar parameters and a vector of filter
+// coefficients.
+bool FieldsOfExpertsCost::Evaluate(double const* const* parameters,
+                                   double* residuals,
+                                   double** jacobians) const {
+  int num_variables = filter_.size();
+  residuals[0] = 0;
+  for (int i = 0; i < num_variables; ++i) {
+    residuals[0] += filter_[i] * parameters[i][0];
+  }
+
+  if (jacobians != NULL) {
+    for (int i = 0; i < num_variables; ++i) {
+      if (jacobians[i] != NULL) {
+        jacobians[i][0] = filter_[i];
+      }
+    }
+  }
+
+  return true;
+}
+
+// This loss function builds the FoE terms and is equal to
+//
+//   f(x) = alpha_i * log(1 + (1/2)s)
+//
+void FieldsOfExpertsLoss::Evaluate(double sq_norm, double rho[3]) const {
+  const double c = 0.5;
+  const double sum = 1.0 + sq_norm * c;
+  const double inv = 1.0 / sum;
+  // 'sum' and 'inv' are always positive, assuming that 's' is.
+  rho[0] = alpha_ *  log(sum);
+  rho[1] = alpha_ * c * inv;
+  rho[2] = - alpha_ * c * c * inv * inv;
+}
+
+FieldsOfExperts::FieldsOfExperts()
+    :  size_(0), num_filters_(0) {
+}
+
+bool FieldsOfExperts::LoadFromFile(const std::string& filename) {
+  std::ifstream foe_file(filename.c_str());
+  foe_file >> size_;
+  foe_file >> num_filters_;
+  if (size_ < 0 || num_filters_ < 0) {
+    return false;
+  }
+  const int num_variables = NumVariables();
+
+  x_delta_indices_.resize(num_variables);
+  for (int i = 0; i < num_variables; ++i) {
+    foe_file >> x_delta_indices_[i];
+  }
+
+  y_delta_indices_.resize(NumVariables());
+  for (int i = 0; i < num_variables; ++i) {
+    foe_file >> y_delta_indices_[i];
+  }
+
+  alpha_.resize(num_filters_);
+  for (int i = 0; i < num_filters_; ++i) {
+    foe_file >> alpha_[i];
+  }
+
+  filters_.resize(num_filters_);
+  for (int i = 0; i < num_filters_; ++i) {
+    filters_[i].resize(num_variables);
+    for (int j = 0; j < num_variables; ++j) {
+      foe_file >> filters_[i][j];
+    }
+  }
+
+  // If any read failed, return failure.
+  if (!foe_file) {
+    size_ = 0;
+    return false;
+  }
+
+  // There cannot be anything else in the file. Try reading another number and
+  // return failure if that succeeded.
+  double temp;
+  foe_file >> temp;
+  if (foe_file) {
+    size_ = 0;
+    return false;
+  }
+
+  return true;
+}
+
+ceres::CostFunction* FieldsOfExperts::NewCostFunction(int alpha_index) const {
+  return new FieldsOfExpertsCost(filters_[alpha_index]);
+}
+
+ceres::LossFunction* FieldsOfExperts::NewLossFunction(int alpha_index) const {
+  return new FieldsOfExpertsLoss(alpha_[alpha_index]);
+}
+
+
+}  // namespace examples
+}  // namespace ceres

examples/fields_of_experts.h

+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2012 Google Inc. All rights reserved.
+// http://code.google.com/p/ceres-solver/
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions are met:
+//
+// * Redistributions of source code must retain the above copyright notice,
+//   this list of conditions and the following disclaimer.
+// * Redistributions in binary form must reproduce the above copyright notice,
+//   this list of conditions and the following disclaimer in the documentation
+//   and/or other materials provided with the distribution.
+// * Neither the name of Google Inc. nor the names of its contributors may be
+//   used to endorse or promote products derived from this software without
+//   specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+// POSSIBILITY OF SUCH DAMAGE.
+//
+// Author: strandmark@google.com (Petter Strandmark)
+//
+// Class for loading the data required for descibing a Fields of Experts (FoE)
+// model. The Fields of Experts regularization consists of terms of the type
+//
+//   alpha * log(1 + (1/2)*sum(F .* X)^2),
+//
+// where F is a d-by-d image patch and alpha is a constant. This is implemented
+// by a FieldsOfExpertsSum object which represents the dot product between the
+// image patches and a FieldsOfExpertsLoss which implements the log(1 + (1/2)s)
+// part.
+//
+// [1] S. Roth and M.J. Black. "Fields of Experts." International Journal of
+//     Computer Vision, 82(2):205--229, 2009.
+
+#ifndef CERES_EXAMPLES_FIELDS_OF_EXPERTS_H_
+#define CERES_EXAMPLES_FIELDS_OF_EXPERTS_H_
+
+#include <iostream>
+#include <vector>
+
+#include "ceres/loss_function.h"
+#include "ceres/cost_function.h"
+#include "ceres/sized_cost_function.h"
+
+#include "pgm_image.h"
+
+namespace ceres {
+namespace examples {
+
+// One sum in the FoE regularizer. This is a dot product between a filter and an
+// image patch. It simply calculates the dot product between the filter
+// coefficients given in the constructor and the scalar parameters passed to it.
+class FieldsOfExpertsCost : public ceres::CostFunction {
+ public:
+  explicit FieldsOfExpertsCost(const std::vector<double>& filter);
+  // The number of scalar parameters passed to Evaluate must equal the number of
+  // filter coefficients passed to the constructor.
+  virtual bool Evaluate(double const* const* parameters,
+                        double* residuals,
+                        double** jacobians) const;
+
+ private:
+  const std::vector<double>& filter_;
+};
+
+// The loss function used to build the correct regularization. See above.
+//
+//   f(x) = alpha_i * log(1 + (1/2)s)
+//
+class FieldsOfExpertsLoss : public ceres::LossFunction {
+ public:
+  explicit FieldsOfExpertsLoss(double alpha) : alpha_(alpha) { }
+  virtual void Evaluate(double, double*) const;
+
+ private:
+  const double alpha_;
+};
+
+// This class loads a set of filters and coefficients from file. Then the users
+// obtains the correct loss and cost functions through NewCostFunction and
+// NewLossFunction.
+class FieldsOfExperts {
+ public:
+  // Creates an empty object with size() == 0.
+  FieldsOfExperts();
+  // Attempts to load filters from a file. If unsuccessful it returns false and
+  // sets size() == 0.
+  bool LoadFromFile(const std::string& filename);
+
+  // Side length of a square filter in this FoE. They are all of the same size.
+  int Size() const {
+    return size_;
+  }
+
+  // Total number of pixels the filter covers.
+  int NumVariables() const {
+    return size_ * size_;
+  }
+
+  // Number of filters used by the FoE.
+  int NumFilters() const {
+    return num_filters_;
+  }
+
+  // Creates a new cost function. The caller is responsible for deallocating the
+  // memory. alpha_index specifies which filter is used in the cost function.
+  ceres::CostFunction* NewCostFunction(int alpha_index) const;
+  // Creates a new loss function. The caller is responsible for deallocating the
+  // memory. alpha_index specifies which filter this loss function is for.
+  ceres::LossFunction* NewLossFunction(int alpha_index) const;
+
+  // Gets the delta pixel indices for all pixels in a patch.
+  const std::vector<int>& GetXDeltaIndices() const {
+    return x_delta_indices_;
+  }
+  const std::vector<int>& GetYDeltaIndices() const {
+    return y_delta_indices_;
+  }
+
+ private:
+  // The side length of a square filter.
+  int size_;
+  // The number of different filters used.
+  int num_filters_;
+  // Pixel offsets for all variables.
+  std::vector<int> x_delta_indices_, y_delta_indices_;
+  // The coefficients in front of each term.
+  std::vector<double> alpha_;
+  // The filters used for the dot product with image patches.
+  std::vector<std::vector<double> > filters_;
+};
+
+}  // namespace examples
+}  // namespace ceres
+
+#endif  // CERES_EXAMPLES_FIELDS_OF_EXPERTS_H_

examples/pgm_image.h

+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2012 Google Inc. All rights reserved.
+// http://code.google.com/p/ceres-solver/
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions are met:
+//
+// * Redistributions of source code must retain the above copyright notice,
+//   this list of conditions and the following disclaimer.
+// * Redistributions in binary form must reproduce the above copyright notice,
+//   this list of conditions and the following disclaimer in the documentation
+//   and/or other materials provided with the distribution.
+// * Neither the name of Google Inc. nor the names of its contributors may be
+//   used to endorse or promote products derived from this software without
+//   specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+// POSSIBILITY OF SUCH DAMAGE.
+//
+// Author: strandmark@google.com (Petter Strandmark)
+//
+// Simple class for accessing PGM images.
+
+#ifndef CERES_EXAMPLES_PGM_IMAGE_H_
+#define CERES_EXAMPLES_PGM_IMAGE_H_
+
+#include <algorithm>
+#include <cstring>
+#include <fstream>
+#include <iostream>
+#include <sstream>
+#include <string>
+#include <vector>
+
+#include "glog/logging.h"
+
+namespace ceres {
+namespace examples {
+
+template<typename Real>
+class PGMImage {
+ public:
+  // Create an empty image
+  PGMImage(int width, int height);
+  // Load an image from file
+  explicit PGMImage(std::string filename);
+  // Sets an image to a constant
+  void Set(double constant);
+
+  // Reading dimensions
+  int width() const;
+  int height() const;
+  int NumPixels() const;
+
+  // Get individual pixels
+  Real* MutablePixel(int x, int y);
+  Real  Pixel(int x, int y) const;
+  Real* MutablePixelFromLinearIndex(int index);
+  Real  PixelFromLinearIndex(int index) const;
+  int LinearIndex(int x, int y) const;
+
+  // Adds an image to another
+  void operator+=(const PGMImage& image);
+  // Adds a constant to an image
+  void operator+=(Real a);
+  // Multiplies an image by a constant
+  void operator*=(Real a);
+
+  // File access
+  bool WriteToFile(std::string filename) const;
+  bool ReadFromFile(std::string filename);
+
+  // Accessing the image data directly
+  bool SetData(const std::vector<Real>& new_data);
+  const std::vector<Real>& data() const;
+
+ protected:
+  int height_, width_;
+  std::vector<Real> data_;
+};
+
+// --- IMPLEMENTATION
+
+template<typename Real>
+PGMImage<Real>::PGMImage(int width, int height)
+  : height_(height), width_(width), data_(width*height, 0.0) {
+}
+
+template<typename Real>
+PGMImage<Real>::PGMImage(std::string filename) {
+  if (!ReadFromFile(filename)) {
+    height_ = 0;
+    width_  = 0;
+  }
+}
+
+template<typename Real>
+void PGMImage<Real>::Set(double constant) {
+  for (int i = 0; i < data_.size(); ++i) {
+    data_[i] = constant;
+  }
+}
+
+template<typename Real>
+int PGMImage<Real>::width() const {
+  return width_;
+}
+
+template<typename Real>
+int PGMImage<Real>::height() const {
+  return height_;
+}
+
+template<typename Real>
+int PGMImage<Real>::NumPixels() const {
+  return width_ * height_;
+}
+
+template<typename Real>
+Real* PGMImage<Real>::MutablePixel(int x, int y) {
+  return MutablePixelFromLinearIndex(LinearIndex(x, y));
+}
+
+template<typename Real>
+Real PGMImage<Real>::Pixel(int x, int y) const {
+  return PixelFromLinearIndex(LinearIndex(x, y));
+}
+
+template<typename Real>
+Real* PGMImage<Real>::MutablePixelFromLinearIndex(int index) {
+  CHECK(index >= 0);
+  CHECK(index < width_ * height_);
+  CHECK(index < data_.size());
+  return &data_[index];
+}
+
+template<typename Real>
+Real  PGMImage<Real>::PixelFromLinearIndex(int index) const {
+  CHECK(index >= 0);
+  CHECK(index < width_ * height_);
+  CHECK(index < data_.size());
+  return data_[index];
+}
+
+template<typename Real>
+int PGMImage<Real>::LinearIndex(int x, int y) const {
+  return x + width_*y;
+}
+
+// Adds an image to another
+template<typename Real>
+void PGMImage<Real>::operator+= (const PGMImage<Real>& image) {
+  CHECK(data_.size() == image.data_.size());
+  for (int i = 0; i < data_.size(); ++i) {
+    data_[i] += image.data_[i];
+  }
+}
+
+// Adds a constant to an image
+template<typename Real>
+void PGMImage<Real>::operator+= (Real a) {
+  for (int i = 0; i < data_.size(); ++i) {
+    data_[i] += a;
+  }
+}
+
+// Multiplies an image by a constant
+template<typename Real>
+void PGMImage<Real>::operator*= (Real a) {
+  for (int i = 0; i < data_.size(); ++i) {
+    data_[i] *= a;
+  }
+}
+
+template<typename Real>
+bool PGMImage<Real>::WriteToFile(std::string filename) const {
+  std::ofstream outputfile(filename.c_str());
+  outputfile << "P2" << std::endl;
+  outputfile << "# PGM format" << std::endl;
+  outputfile << " # <width> <height> <levels> " << std::endl;
+  outputfile << " # <data> ... " << std::endl;
+  outputfile << width_ << ' ' << height_ << " 255 " << std::endl;
+
+  // Write data
+  int num_pixels = width_*height_;
+  for (int i = 0; i < num_pixels; ++i) {
+    // Convert to integer by rounding when writing file
+    outputfile << static_cast<int>(data_[i] + 0.5) << ' ';
+  }
+
+  return outputfile;  // Returns true/false
+}
+
+namespace  {
+
+// Helper function to read data from a text file, ignoring "#" comments.
+template<typename T>
+bool GetIgnoreComment(std::istream* in, T& t) {
+  std::string word;
+  bool ok;
+  do {
+    ok = true;
+    (*in) >> word;
+    if (word.length() > 0 && word[0] == '#') {
+      // Comment; read the whole line
+      ok = false;
+      std::getline(*in, word);
+    }
+  } while (!ok);
+
+  // Convert the string
+  std::stringstream sin(word);
+  sin >> t;
+
+  // Check for success
+  if (!in || !sin) {
+    return false;
+  }
+  return true;
+}
+}  // namespace
+
+template<typename Real>
+bool PGMImage<Real>::ReadFromFile(std::string filename) {
+  std::ifstream inputfile(filename.c_str());
+
+  // File must start with "P2"
+  char ch1, ch2;
+  inputfile >> ch1 >> ch2;
+  if (!inputfile || ch1 != 'P' || (ch2 != '2' && ch2 != '5')) {
+    return false;
+  }
+
+  // Read the image header
+  int two_fifty_five;
+  if (!GetIgnoreComment(&inputfile, width_)  ||
+      !GetIgnoreComment(&inputfile, height_) ||
+      !GetIgnoreComment(&inputfile, two_fifty_five) ) {
+    return false;
+  }
+  // Assert that the number of grey levels is 255.
+  if (two_fifty_five != 255) {
+    return false;
+  }
+
+  // Now read the data
+  int num_pixels = width_*height_;
+  data_.resize(num_pixels);
+  if (ch2 == '2') {
+    // Ascii file
+    for (int i = 0; i < num_pixels; ++i) {
+      int pixel_data;
+      bool res = GetIgnoreComment(&inputfile, pixel_data);
+      if (!res) {
+        return false;
+      }
+      data_[i] = pixel_data;
+    }
+    // There cannot be anything else in the file (except comments). Try reading
+    // another number and return failure if that succeeded.
+    int temp;
+    bool res = GetIgnoreComment(&inputfile, temp);
+    if (res) {
+      return false;
+    }
+  } else {
+    // Read the line feed character
+    if (inputfile.get() != '\n') {
+      return false;
+    }
+    // Binary file
+    // TODO(strandmark): Will not work on Windows (linebreak conversion).
+    for (int i = 0; i < num_pixels; ++i) {
+      unsigned char pixel_data = inputfile.get();
+      if (!inputfile) {
+        return false;
+      }
+      data_[i] = pixel_data;
+    }
+    // There cannot be anything else in the file. Try reading another byte
+    // and return failure if that succeeded.
+    inputfile.get();
+    if (inputfile) {
+      return false;
+    }
+  }
+
+  return true;
+}
+
+template<typename Real>
+bool PGMImage<Real>::SetData(const std::vector<Real>& new_data) {
+  // This function cannot change the dimensions
+  if (new_data.size() != data_.size()) {
+    return false;
+  }
+  std::copy(new_data.begin(), new_data.end(), data_.begin());
+  return true;
+}
+
+template<typename Real>
+const std::vector<Real>& PGMImage<Real>::data() const {
+  return data_;
+}
+
+}  // namespace examples
+}  // namespace ceres
+
+
+#endif  // CERES_EXAMPLES_PGM_IMAGE_H_