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
diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
index 4b2938b..2307a03 100644
--- a/examples/CMakeLists.txt
+++ b/examples/CMakeLists.txt
@@ -54,6 +54,11 @@
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
diff --git a/examples/denoising.cc b/examples/denoising.cc
new file mode 100644
index 0000000..086be00
--- /dev/null
+++ b/examples/denoising.cc
@@ -0,0 +1,214 @@
+// 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;
+}
diff --git a/examples/fields_of_experts.cc b/examples/fields_of_experts.cc
new file mode 100644
index 0000000..0cee40b
--- /dev/null
+++ b/examples/fields_of_experts.cc
@@ -0,0 +1,152 @@
+// 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
diff --git a/examples/fields_of_experts.h b/examples/fields_of_experts.h
new file mode 100644
index 0000000..845a4cf
--- /dev/null
+++ b/examples/fields_of_experts.h
@@ -0,0 +1,145 @@
+// 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_
diff --git a/examples/pgm_image.h b/examples/pgm_image.h
new file mode 100644
index 0000000..15e99e4
--- /dev/null
+++ b/examples/pgm_image.h
@@ -0,0 +1,319 @@
+// 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_
