criticalKernelsThinning3D.cpp

 
#include <iostream>
#include <chrono>
#include <unordered_map>
#ifdef WITH_QGLVIEWER
#include <DGtal/io/viewers/Viewer3D.h>
#endif
#include <DGtal/base/Common.h>
#include <DGtal/helpers/StdDefs.h>
#include <DGtal/io/Color.h>
#include <DGtal/io/readers/GenericReader.h>
#include <DGtal/io/writers/GenericWriter.h>
#include "DGtal/images/imagesSetsUtils/SetFromImage.h"
#include "DGtal/images/SimpleThresholdForegroundPredicate.h"
#include "DGtal/images/ImageSelector.h"
#include "DGtal/images/imagesSetsUtils/ImageFromSet.h"
 
#include <DGtal/topology/SurfelAdjacency.h>
#include <DGtal/topology/CubicalComplex.h>
#include <DGtal/topology/CubicalComplexFunctions.h>
#include <DGtal/io/boards/Board3D.h>
 
#include <DGtal/topology/VoxelComplex.h>
#include <DGtal/topology/VoxelComplexFunctions.h>
#include "DGtal/topology/NeighborhoodConfigurations.h"
#include "DGtal/topology/tables/NeighborhoodTables.h"
 
// Distance Map
#include "DGtal/geometry/volumes/distance/ExactPredicateLpSeparableMetric.h"
#include "DGtal/geometry/volumes/distance/VoronoiMap.h"
#include "DGtal/geometry/volumes/distance/DistanceTransformation.h"
 
#include "CLI11.hpp"
 
 
using namespace DGtal;
using namespace std;
using namespace DGtal::Z3i;
 
int main(int argc, char* const argv[]){
  
  // parse command line using CLI ----------------------------------------------
  CLI::App app;
  std::string inputFileName;
  std::string sk_string {"1isthmus"};
  string select_string {"dmax"};
  string foreground {"black"};
  string outputFilenameImg;
  string outputFilenameSDP;
  string outputFilenameOBJ;
  string outputFilenameInputOBJ;
  
  int thresholdMin {0};
  int thresholdMax {255};
  int persistence {0};
  bool profile {false};
  bool verbose  {false};
  bool visualize {false};
  bool useInputImgToExp {false};
  
  app.description("Compute the thinning of a volume using the CriticalKernels framework\nBasic usage: criticalKernelsThinning3D --input <volFileName> --skel <ulti,end, 1isthmus, isthmus> --select [ -f <white,black> -m <minlevel> -M <maxlevel> -v ] [--persistence <value> ] --persistence <value> ] \n options for --skel {ulti end 1isthmus isthmus} \n options for --select = {dmax random first} \n Example: \n criticalKernelsThinning3D --input ${DGtal}/examples/samples/Al.100.vol --select dmax --skel 1isthmus --persistence 1 --visualize --verbose --exportImage ./Al100_dmax_1isthmus_p1.vol \n");
  app.add_option("-i,--input,1", inputFileName, "Input vol file." )
  ->required()
  ->check(CLI::ExistingFile);
 
  app.add_option("--skel,-s", sk_string,"Type of skeletonization. Options: 1isthmus, isthmus, end, ulti.", true )
   -> check(CLI::IsMember({"ulti", "end","isthmus", "1isthmus"}));
  app.add_option("--select,-c", select_string, "Select the ordering for skeletonization. Options: dmax, random, first", true)
   -> check(CLI::IsMember({"random", "dmax", "first"}));
  app.add_option("--foreground,-f",foreground, "Foreground color in binary image", true )
   -> check(CLI::IsMember({"white", "black"}));
  app.add_option("--thresholdMin,-m", thresholdMin, "Threshold min (excluded) to define binary shape", true );
  app.add_option("--thresholdMax,-M", thresholdMax, "Threshold max (included) to define binary shape", true );
  app.add_option("--persistence,-p",persistence,"Persistence value, implies use of persistence algorithm if p>=1", true )
  ->check(CLI::PositiveNumber);
  app.add_flag("--profile", profile, "Profile algorithm");
  app.add_flag("--verbose,-v",verbose, "Verbose output");
  app.add_option("--exportImage,-o",outputFilenameImg, "Export the resulting set of points to a image compatible with GenericWriter.");
  app.add_option("--exportSDP,-e",outputFilenameSDP, "Export the resulting set of points in a simple (sequence of discrete point (sdp))." );
  app.add_option("--exportOBJ,-O",outputFilenameOBJ, "Export the resulting set of points in an OBJ file." );
  app.add_option("--exportInputOBJ,-I",outputFilenameInputOBJ, "Export the input set of points in an OBJ file." );
#ifdef WITH_QGLVIEWER
  app.add_flag("--visualize,-t", visualize, "Visualize result in viewer");
#endif
  app.add_flag("--useInputImgToExp,-k", useInputImgToExp, "Use input image type to export result (allowing to keep same domain (and same image spacing when using ITK)).");
    
  app.get_formatter()->column_width(40);
  CLI11_PARSE(app, argc, argv);
  // END parse command line using CLI ----------------------------------------------
 
  
  
  if(verbose){
    std::cout << "Skel: " << sk_string << std::endl;
    std::cout << "Select: " << select_string << std::endl;
    std::cout << "Persistence: " << persistence << std::endl;
    std::cout << "Input: " << inputFileName << std::endl;
  }
  trace.beginBlock("Reading input");
  using Domain = Z3i::Domain ;
 
#ifdef WITH_ITK
  using Image = ImageSelector < Z3i::Domain, unsigned char, ITKIMAGEDATA_CONTAINER_I>::Type ;
#else
  using Image = ImageSelector < Z3i::Domain, unsigned char>::Type ;
#endif
 
  
  Image image = GenericReader<Image>::import(inputFileName);
  trace.endBlock();
 
  DigitalSet image_set (image.domain());
  SetFromImage<Z3i::DigitalSet>::append<Image>(
      image_set, image,
      thresholdMin, thresholdMax);
 
 
  // Create a VoxelComplex from the set
 
  using DigitalTopology = DT26_6;
  using DigitalSet =
    DGtal::DigitalSetByAssociativeContainer<Domain ,
      std::unordered_set< typename Domain::Point> >;
  using Complex = DGtal::VoxelComplex<KSpace>;
 
  auto & sk = sk_string;
  KSpace ks;
  KSpace::Point d1( KSpace::Point::diagonal( 1 ) );
  ks.init(image.domain().lowerBound() - d1 ,
      image.domain().upperBound() + d1 , true);
 
  trace.beginBlock("construct with table");
  Complex vc(ks);
  vc.construct(image_set, functions::loadTable(simplicity::tableSimple26_6 ));
  trace.endBlock();
  trace.beginBlock("load isthmus table");
  boost::dynamic_bitset<> isthmus_table;
  if (sk == "isthmus")
    isthmus_table = *functions::loadTable(isthmusicity::tableIsthmus);
  else if (sk == "1isthmus")
    isthmus_table = *functions::loadTable(isthmusicity::tableOneIsthmus);
  auto pointMap = *functions::mapZeroPointNeighborhoodToConfigurationMask<Point>();
 
  trace.endBlock();
  using namespace DGtal::functions ;
  // SKEL FUNCTION:
  std::function< bool(const Complex&, const Cell&) > Skel ;
  if (sk == "ulti") Skel = skelUltimate<Complex>;
  else if (sk == "end") Skel = skelEnd<Complex>;
  else if (sk == "isthmus" || sk == "1isthmus")
      Skel = [&isthmus_table, &pointMap](const Complex & fc,
               const Complex::Cell & c){
        return skelWithTable(isthmus_table, pointMap, fc, c);
      };
  else throw std::runtime_error("Invalid skel string");
  auto start = std::chrono::system_clock::now();
 
  // SELECT FUNCTION
  /*
   * Calculate distance map even if not requested:
   */
  trace.beginBlock("Create Distance Map");
  using L3Metric = ExactPredicateLpSeparableMetric<Z3i::Space, 3>;
  using DT       = DistanceTransformation<Z3i::Space, DigitalSet, L3Metric>;
  L3Metric l3;
  DT dt(image.domain(), image_set, l3);
  trace.endBlock();
 
  std::function< std::pair<typename Complex::Cell, typename Complex::Data>(const Complex::Clique&) > Select ;
  auto & sel = select_string;
  if (sel == "random") Select = selectRandom<Complex>;
  else if (sel == "first") Select = selectFirst<Complex>;
  else if (sel == "dmax"){
    Select =
      [&dt](const Complex::Clique & clique){
        return selectMaxValue<DT, Complex>(dt,clique);
      };
  } else throw std::runtime_error("Invalid skel string");
 
  trace.beginBlock("Thinning");
  Complex vc_new(ks);
  if (persistence == 0)
    vc_new = asymetricThinningScheme< Complex >(
        vc, Select,  Skel, verbose);
  else
    vc_new = persistenceAsymetricThinningScheme< Complex >(
        vc, Select,  Skel, persistence, verbose);
  trace.endBlock();
 
  auto end = std::chrono::system_clock::now();
  auto elapsed = std::chrono::duration_cast<std::chrono::seconds> (end - start) ;
  if (profile) std::cout <<"Time elapsed: " << elapsed.count() << std::endl;
 
 
  DigitalSet thin_set(image.domain());
  vc_new.dumpVoxels(thin_set);
  const auto & all_set = image_set;
 
  if (outputFilenameSDP != "")
  {
    std::ofstream out;
    out.open(outputFilenameSDP.c_str());
    for (auto &p : thin_set)
    {
      out << p[0] << " " << p[1] << " " << p[2] << std::endl;
    }
  }
 
 
  //-------------- export OBJ -------------------------------------------
  if ( outputFilenameInputOBJ != "" )
    {
      Board3D< Space,KSpace > board( ks );
      // Display lines that are not in the mesh.
      SCell surfel;
      board << SetMode3D( surfel.className(), "Basic");
      board.setLineColor( Color::Black );
      board.setFillColor( Color( 200, 200, 255, 255 ) );
      for ( auto p : all_set )
        {
          SCell voxel  = ks.sSpel( p );
          for ( Dimension k = 0; k < 3; k++ ) {
            Point q = p;
            q[ k ] += 1;
            if ( all_set.find( q ) == all_set.end() )
              board << ks.sIncident( voxel, k, true );
            q[ k ] -= 2;
            if ( all_set.find( q ) == all_set.end() )
              board << ks.sIncident( voxel, k, false );
          }
        }
      board.saveOBJ( outputFilenameInputOBJ );
    }
 
  //-------------- export OBJ -------------------------------------------
  if ( outputFilenameOBJ != "" )
    {
      Board3D< Space,KSpace > board( ks );
      // Display lines that are not in the mesh.
      SCell surfel;
      board << SetMode3D( surfel.className(), "Basic");
      board.setLineColor( Color::Black );
      board.setFillColor( Color::Red );
      for ( auto p : thin_set )
        {
          SCell voxel  = ks.sSpel( p );
          for ( Dimension k = 0; k < 3; k++ ) {
            Point q = p;
            q[ k ] += 1;
            if ( thin_set.find( q ) == thin_set.end() )
              board << ks.sIncident( voxel, k, true );
            q[ k ] -= 2;
            if ( thin_set.find( q ) == thin_set.end() )
              board << ks.sIncident( voxel, k, false );
          }
        }
      board.saveOBJ( outputFilenameOBJ );
    }
 
if (outputFilenameImg != "")
  {
    if(verbose)
      std::cout << "outputFilename" << outputFilenameImg << std::endl;
 
    unsigned int foreground_value = 255;
    if (useInputImgToExp){
      for(auto p: image.domain()){image.setValue(p, 0);}
      ImageFromSet<Image>::append(image, thin_set, foreground_value);
      image >> outputFilenameImg;
    }else{
      auto thin_image = ImageFromSet<Image>::create(thin_set, foreground_value, false, useInputImgToExp);
      thin_image >> outputFilenameImg;
    }
  }
#ifdef WITH_QGLVIEWER
if(visualize)
  {
    int argc(1);
    char** argv(nullptr);
    QApplication app(argc, argv);
    Viewer3D<> viewer( ks );
    viewer.setWindowTitle("criticalKernelsThinning3D");
    viewer.show();
 
    viewer.setFillColor(Color(255, 255, 255, 255));
    viewer << thin_set;
 
    // All kspace voxels
    viewer.setFillColor(Color(40, 200, 55, 10));
    viewer << all_set;
 
    viewer << Viewer3D<>::updateDisplay;
 
    app.exec();
  }
#endif
  
}