eulerCharacteristic.cpp

#include <iostream>
#include <DGtal/base/Common.h>
#include <DGtal/helpers/StdDefs.h>
#include <DGtal/io/readers/VolReader.h>

#include "CLI11.hpp"

#include <DGtal/images/ImageContainerBySTLVector.h>
#include <DGtal/images/IntervalForegroundPredicate.h>

using namespace std;
using namespace DGtal;
using namespace Z3i;

int main(int argc, char**argv)
{

  // parse command line CLI ----------------------------------------------
  CLI::App app;
  std::string filename;
  int thresholdMin {0};
  int thresholdMax {255};

  app.description("Computes the Euleur Characteristic of  a vol to a 8-bit raw file.\n Typical use example:\n \t eulerCharacteristic <volFileName> -m <minlevel> -M <maxlevel>\n");
  app.add_option("--input,-i,1", filename, "Input vol file." )->required()->check(CLI::ExistingFile);
  app.add_option("--thresholdMin,-m", thresholdMin, "threshold min (excluded) to define binary shape (default 0)", true);
  app.add_option("--thresholdMax,-M", thresholdMax, "threshold max (included) to define binary shape (default 255)", true);
  
  app.get_formatter()->column_width(40);
  CLI11_PARSE(app, argc, argv);
  // END parse command line using CLI ----------------------------------------------
 
  //Importing the Vol
  trace.beginBlock("Loading the vol file");
  typedef ImageContainerBySTLVector<Z3i::Domain, unsigned char>  MyImageC;
  MyImageC  imageC = VolReader< MyImageC >::importVol ( filename );
  trace.info()<<imageC<<std::endl;
  trace.endBlock();  

  //Constructing the cubical complex
  trace.beginBlock("Construting the cubical complex");
  KSpace::CellSet myCellSet;
  KSpace  ks;
  bool space_ok = ks.init( imageC.domain().lowerBound(), imageC.domain().upperBound(), true );
  if (!space_ok)
    {
      trace.error() << "Error in the Khamisky space construction."<<std::endl;
      return 2;
    }
  functors::IntervalForegroundPredicate<MyImageC> interval(imageC, thresholdMin,thresholdMax);  
  for(MyImageC::Domain::ConstIterator it =imageC.domain().begin(), itend= imageC.domain().end();
      it != itend; ++it)
    {
      if (interval( *it ))
        {
          Domain dom( 2*(*it), 2*(*it) + Point::diagonal(2));
          for(Domain::ConstIterator itdom = dom.begin(), itdomend = dom.end(); itdom != itdomend; ++itdom)
            myCellSet.insert( ks.uCell( *itdom) );
        }
    }
  trace.info() << "Got "<< myCellSet.size()<< " cells"<<std::endl;
  trace.endBlock();

  trace.beginBlock("Computing the characteristics");
  std::vector<int> cells(4,0);
  
  for(KSpace::CellSet::const_iterator it = myCellSet.begin(), itend = myCellSet.end(); it !=itend; ++it)
    cells[ ks.uDim(*it) ] ++; 
  
  trace.info() << "Got "<< cells[0]<< " pointels "<<cells[1]<<" linels  "<< cells[2]<<" surfels and "<<cells[3]<<"  bells"<<std::endl;
  trace.endBlock();

  trace.info() << "Volumetric Euler Characteristic = "<<cells[0] - cells[1] + cells[2] - cells[3]<<std::endl;

  return 0;
}