2dCompImage.cpp

#include "DGtal/base/Common.h"
#include "DGtal/helpers/StdDefs.h"
#include <DGtal/io/readers/GenericReader.h>
#include <DGtal/io/writers/GenericWriter.h>


#include <DGtal/images/Image.h>
#include <DGtal/images/ImageContainerBySTLVector.h>
#include <DGtal/images/imagesSetsUtils/SetFromImage.h>
#include <DGtal/io/writers/PPMWriter.h>
#include <DGtal/math/Statistic.h>
#include "DGtal/io/colormaps/GradientColorMap.h"
#include <DGtal/math/Statistic.h>

#include "CLI11.hpp"


using namespace std;
using namespace DGtal;

typedef ImageContainerBySTLVector < Z2i::Domain,  unsigned char > Image2D;
typedef ImageContainerBySTLVector < Z2i::Domain,  unsigned int > Image2DErr;
typedef GradientColorMap<unsigned int, CMAP_JET, 1 > JetMap;



Statistic< int>
getMAEstats(const Image2D & imageA, const Image2D &imageB, Image2DErr &imageMAE)
{
  Statistic< int> stat(false);
  for(auto const &point: imageA.domain())
  {
    unsigned int error = abs((imageA(point)-imageB(point)));
    stat.addValue(error);
    imageMAE.setValue(point, error);     
  }
  stat.terminate();
  return stat;
}


Statistic< int>
getMSEstats(const Image2D & imageA, const Image2D &imageB, Image2DErr &imageMSE)
{
  Statistic< int> stat(false);
  for(Image2D::Domain::ConstIterator it = imageA.domain().begin();
      it!=imageA.domain().end(); it++)
  {
    int error = (imageA(*it)-imageB(*it))*(imageA(*it)-imageB(*it));
    stat.addValue(error);
    imageMSE.setValue(*it, error);     
  }
  stat.terminate();
  return stat;
}



template<typename StatisticT>
void
displayStats(const StatisticT &aStat, const string &type)
{
  std::cout << "# Stats on " << type << ": mean max min unbiased_variance nb_samples  " << std::endl; 
  std::cout << aStat.mean() << " " << aStat.max() << " " << aStat.min() << " "
            << " " << aStat.unbiasedVariance() << " " << aStat.samples() << std::endl;  
}


int main( int argc, char** argv )
{
  
  // parse command line using CLI ----------------------------------------------
  CLI::App app;
  std::string inputFileNameA;
  std::string inputFileNameB;
  std::string basenameOutput;
  int maxValueMAE;
  int maxValueMSE;
  
  app.description("Compare images and displays differences (squared and absolute differences).\n Typical use example:\n \t 2dCompImage imageA.pgm imageB.pgm -e errorImage -S 100 \n");
  app.add_option("-a,--imageA,1", inputFileNameA, "Input filename of image A." )
  ->required()
  ->check(CLI::ExistingFile);
  app.add_option("-b,--imageB,2", inputFileNameA, "Input filename of image B." )
   ->required()
   ->check(CLI::ExistingFile);

  
  app.add_option("--imageError,-e",basenameOutput, "Output error image basename (will generate two images <basename>MSE.ppm and <basename>MAE.ppm).");
  
  auto setMaxMSEOpt = app.add_option("--setMaxColorValueMSE,-S",maxValueMSE, "Set the maximal color value for the scale display of MSE (else the scale is set the maximal MSE value).");
  auto setMaxMAEOpt = app.add_option("--setMaxColorValueMAE,-A",maxValueMAE, "Set the maximal color value for the scale display of MAE (else the scale is set from the maximal MAE value).");

  
  app.get_formatter()->column_width(40);
  CLI11_PARSE(app, argc, argv);
  // END parse command line using CLI ----------------------------------------------

  
  //  Input images ----------------------------------------------------
  Image2D imageA = GenericReader<Image2D>::import(inputFileNameA);
  Image2D imageB = GenericReader<Image2D>::import(inputFileNameB);
  
  Image2DErr imageErr (imageA.domain());  


  //  Absolute Error between imageA and imageB -------------------------
  Statistic<int> statMA = getMAEstats(imageA, imageB, imageErr); 
  int maxVal = statMA.max();
  if(setMaxMAEOpt->count() > 0)
  {
    maxVal = maxValueMAE;
  }
  JetMap jmapMA(0, maxVal);
  displayStats(statMA, "Absolute errror");
  stringstream maeName; maeName << basenameOutput;
  maeName << "MAE.ppm";
  PPMWriter<Image2DErr, JetMap>::exportPPM(maeName.str(), imageErr,  jmapMA);

  //  Squared Error between imageA and imageB -------------------------
  Statistic<int> statSE = getMSEstats(imageA, imageB, imageErr); 
  maxVal = statMA.max();
  if(setMaxMSEOpt->count()>0)
  {
    maxVal = maxValueMSE;
  }    
  JetMap jmapSE(0, maxVal);
  displayStats(statSE, "Squared error");
  stringstream mseName; mseName << basenameOutput;
  mseName << "MSE.ppm";
  PPMWriter<Image2DErr, JetMap>::exportPPM(mseName.str(), imageErr,  jmapSE);

  return 0;
}