3dCompSurfelData.cpp

 
 
#include <iostream>
#include <sstream>
#include <stdio.h>
 
#include "DGtal/base/Common.h"
#include "DGtal/helpers/StdDefs.h"
#include "DGtal/io/viewers/PolyscopeViewer.h"
#include "DGtal/io/readers/PointListReader.h"
#include "DGtal/io/readers/MeshReader.h"
 
#include "DGtal/topology/helpers/Surfaces.h"
#include "DGtal/topology/SurfelAdjacency.h"
#include "DGtal/topology/CanonicCellEmbedder.h"
#include "DGtal/math/Statistic.h"
 
#include "DGtal/io/Color.h"
#include "DGtal/io/colormaps/GradientColorMap.h"
#include "DGtal/io/readers/GenericReader.h"
 
#include "CLI11.hpp"
 
 
#include <limits>
 
using namespace std;
using namespace DGtal;
using namespace Z3i;
 
 
 
template < typename Point>
void
getBoundingUpperAndLowerPoint(const std::vector<Point> &vectorPt, Point &ptLower, Point &ptUpper){
  for(unsigned int i =1; i<vectorPt.size(); i++){
    if(vectorPt.at(i)[0] < ptLower[0]){
      ptLower[0] = vectorPt.at(i)[0];
    }
    if(vectorPt.at(i)[1] < ptLower[1]){
      ptLower[1] = vectorPt.at(i)[1];
    }
   if(vectorPt.at(i)[2] < ptLower[2]){
      ptLower[2] =vectorPt.at(i)[2];
    }
   if(vectorPt.at(i)[0] > ptUpper[0]){
      ptUpper[0] = vectorPt.at(i)[0];
   }
   if(vectorPt.at(i)[1] > ptUpper[1]){
     ptUpper[1] = vectorPt.at(i)[1];
    }
   if(vectorPt.at(i)[2] > ptUpper[2]){
      ptUpper[2] =vectorPt.at(i)[2];
    }
  }
}
 
 
int main( int argc, char** argv )
{
  typedef PointVector<4, double> Point4D;
  typedef PointVector<1, int> Point1D;
  
  // parse command line using CLI ----------------------------------------------
  CLI::App app;
  std::string inputFileName;
  std::string referenceFileName;
  std::string fileMeasureOutput;
  std::string snapShotName;
  bool useLabels {false};
  bool drawSurfelAssociations {false};
  unsigned int labelIndex {4};
  std::vector<unsigned int> vectSDPindex;
  
  app.description("Computes generic scalar surfel data comparisons (squared error) (given from an input data file and from a reference one).\n Each surfels are associated to the nearest one of the reference surfels (computed by a 'brut force' search). This association can also be limited to surfel of same label (if available in the data and by using the --compAccordingLabels option). The comparison and surfel association can be displayed and result statistics are saved on output file (--fileMeasureOutput). You can also remove the interactive 3d view by just doing a snapshot and exit with option --doSnapShotAndExit.\n \t  Example of use: \n \t \t 3dCompSurfelData -i surfelCurvatureInput.dat -r surfelCurvatureRef.dat  --fixMaxColorValue 0.8 -d visuSEcurvature.png -o  statMeasures.dat \n \n  \t \t => From the two compared files you should obtain the result of the comparison (statMeasures.dat) with the associated visualisation (visuSEcurvature.png).\n");
  
  app.add_option("-i,--input,1", inputFileName, "input file: sdpa (sequence of discrete points with attribute)" )
  ->required()
  ->check(CLI::ExistingFile);
  
  app.add_option("-r,--reference,2", referenceFileName, "input file: sdpa (sequence of discrete points with attribute)" )
  ->required()
  ->check(CLI::ExistingFile);
  app.add_flag("-l,--compAccordingLabels", useLabels, "apply the comparisos only on points with same labels (by default fifth colomn)" );
  app.add_flag("-a,--drawSurfelAssociations", drawSurfelAssociations, "Draw the surfel association.");
  app.add_option("-o,--fileMeasureOutput",fileMeasureOutput,  "specify the output file to store (append) the error stats else the result is given to std output. ");
  
  auto labOptio = app.add_option("--labelIndex", labelIndex,"set the index of the label (by default set to 4)  " );
  app.add_option("--SDPindex",vectSDPindex, "specify the sdp index (by default 0,1,2,3).")
   ->expected(4);
  
  app.get_formatter()->column_width(40);
  CLI11_PARSE(app, argc, argv);
  // END parse command line using CLI ----------------------------------------------
 
 
 
 
  Z3i::KSpace K;
 
  std::vector<Point4D> surfelAndScalarInput;
  std::vector<Point4D> surfelAndScalarReference;
  std::vector<Point1D> vectLabelsInput;
  std::vector<Point1D> vectLabelsReference;
 
  if(useLabels){
    if(labOptio->count() > 0){
      std::vector<unsigned int > vectIndex;
      vectIndex.push_back(labelIndex);
      vectLabelsInput = PointListReader<Point1D>::getPointsFromFile(inputFileName, vectIndex);
      vectLabelsReference = PointListReader<Point1D>::getPointsFromFile(referenceFileName, vectIndex);
    }else{
      vectLabelsInput = PointListReader<Point1D>::getPointsFromFile(inputFileName);
      vectLabelsReference = PointListReader<Point1D>::getPointsFromFile(referenceFileName);
    }
  }
 
 
  if(vectSDPindex.size() == 4) {
    surfelAndScalarInput = PointListReader<Point4D>::getPointsFromFile(inputFileName, vectSDPindex);
    surfelAndScalarReference = PointListReader<Point4D>::getPointsFromFile(referenceFileName, vectSDPindex);
  }else{
    surfelAndScalarInput = PointListReader<Point4D>::getPointsFromFile(inputFileName);
    surfelAndScalarReference = PointListReader<Point4D>::getPointsFromFile(referenceFileName);
  }
 
 
  Point4D ptLower = surfelAndScalarReference.at(0);
  Point4D ptUpper = surfelAndScalarReference.at(0);
  getBoundingUpperAndLowerPoint(surfelAndScalarReference,  ptLower, ptUpper);
  getBoundingUpperAndLowerPoint(surfelAndScalarInput,  ptLower, ptUpper);
 
  K.init(Z3i::Point(2*(ptLower[0]-1)+1, 2*(ptLower[1]-1)+1, 2*(ptLower[2]-1)+1),
         Z3i::Point(2*(ptUpper[0]+1)+1, 2*(ptUpper[1]+1)+1, 2*(ptUpper[2]+1)+1), true);
 
 
  std::vector<Cell> vectSurfelsReference;
  std::vector<Cell> vectSurfelsInput;
 
  // Construction of the set of surfels
  for(unsigned int i =0; i<surfelAndScalarReference.size(); i++){
    Point4D pt4d = surfelAndScalarReference.at(i);
    Cell c = K.uCell(Z3i::Point(pt4d[0], pt4d[1], pt4d[2]));
    vectSurfelsReference.push_back(c);
  }
  // Construction of the set of surfels
  for(unsigned int i =0; i<surfelAndScalarInput.size(); i++){
    Point4D pt4d = surfelAndScalarInput.at(i);
    Cell c = K.uCell(Z3i::Point(pt4d[0], pt4d[1], pt4d[2]));
    vectSurfelsInput.push_back(c);
  }
 
 
  // Brut force association of each surfel of the input to the reference with the minimal distance.
  // For each surfel of the input we associate an index to the nearest surfel of the reference.
  // if the option --compAccordingLabels is selected then we search only surfel of same label (more precise comparisons in some cases)
 
 
  CanonicCellEmbedder<KSpace> embeder(K);
  std::vector<unsigned int> vectIndexMinToReference;
 
  trace.info() << "Associating each input surfel to reference:" << std::endl;
  trace.progressBar(0, surfelAndScalarInput.size());
 
  for(unsigned int i=0;i <vectSurfelsInput.size(); i++){
    trace.progressBar(i, vectSurfelsInput.size());
    Z3i::RealPoint ptCenterInput = embeder(vectSurfelsInput.at(i));
    unsigned int indexDistanceMin = 0;
    double distanceMin = std::numeric_limits<int>::max() ;
    for(unsigned int j=0; j <vectSurfelsReference.size(); j++){
      if(useLabels && vectLabelsReference.at(j) != vectLabelsInput.at(i)){
        continue;
      }
      Z3i::RealPoint ptCenterRef = embeder(vectSurfelsReference.at(j));
      double distance = (ptCenterRef - ptCenterInput).norm();
      if(distance < distanceMin){
        distanceMin = distance;
        indexDistanceMin = j;
      }
    }
    vectIndexMinToReference.push_back(indexDistanceMin);
 
  }
  trace.progressBar(surfelAndScalarInput.size(), surfelAndScalarInput.size());
  trace.info() << std::endl;
 
 
 
  //-------------------------
  // Displaying input with error and computing statistics
  stringstream s;
  s << "3dCompSurfelData - DGtalTools: ";
  string name = inputFileName.substr(inputFileName.find_last_of("/")+1,inputFileName.size()) ;
  s << " " <<  name << " (W key to display settings)";
  polyscope::options::programName = s.str();
  polyscope::view::setNavigateStyle(polyscope::NavigateStyle::Free);
 
 
  typedef PolyscopeViewer<> Viewer;
  Viewer viewer(K);
  viewer.allowReuseList = true;
 
  Statistic<double> statErrors(true);
  std::ofstream outputStatStream;
  
  if(fileMeasureOutput != ""){
    outputStatStream.open(fileMeasureOutput.c_str(), ios::app );
  }
 
  double maxSqError=0;
  for(unsigned int i=0;i <surfelAndScalarInput.size(); i++){
    double scalarInput = surfelAndScalarInput.at(i)[3];
    double scalarRef = surfelAndScalarReference.at(vectIndexMinToReference.at(i))[3];
    double sqError = (scalarRef-scalarInput)*(scalarRef-scalarInput);
    statErrors.addValue(sqError);
    if(sqError> maxSqError){
      maxSqError =sqError;
    }
  }
 
  trace.info() << "Maximal error:" << maxSqError << std::endl;
  for(unsigned int i=0; i <surfelAndScalarInput.size(); i++){
    double scalarInput = surfelAndScalarInput.at(i)[3];
    double scalarRef = surfelAndScalarReference.at(vectIndexMinToReference.at(i))[3];
    double sqError = (scalarRef-scalarInput)*(scalarRef-scalarInput);
 
    viewer << WithQuantity(vectSurfelsInput.at(i), "Squared Error", sqError);
  }
 
  if(drawSurfelAssociations){
    // Second loop to draw associations
    // Allow for easier groupping with viewer
    for(unsigned int i=0; i <surfelAndScalarInput.size(); i++){
      viewer.drawLine(embeder(vectSurfelsInput.at(i)),embeder(vectSurfelsReference.at(vectIndexMinToReference.at(i))));
    }
  }
 
  statErrors.terminate();
  if(fileMeasureOutput != ""){
    outputStatStream << "input= " <<  inputFileName << " reference=" << referenceFileName << " " ;
    outputStatStream << statErrors << std::endl;
  }else{
    trace.info()  << statErrors;
  }
 
  viewer.show();
  return 0;
}