DGtalTools  1.2.0
3dVolBoundaryViewer.cpp
1 
30 #include <iostream>
31 
32 #include "DGtal/base/Common.h"
33 #include "DGtal/base/BasicFunctors.h"
34 #include "DGtal/kernel/SpaceND.h"
35 #include "DGtal/kernel/domains/HyperRectDomain.h"
36 #include "DGtal/images/ImageSelector.h"
37 #include "DGtal/images/IntervalForegroundPredicate.h"
38 #include "DGtal/topology/KhalimskySpaceND.h"
39 #include "DGtal/topology/DigitalSurface.h"
40 #include "DGtal/topology/SetOfSurfels.h"
41 #include "DGtal/io/viewers/Viewer3D.h"
42 #include "DGtal/io/readers/PointListReader.h"
43 #include "DGtal/io/readers/GenericReader.h"
44 #ifdef WITH_ITK
45 #include "DGtal/io/readers/DicomReader.h"
46 #endif
47 #include "DGtal/io/Color.h"
48 #include "DGtal/io/colormaps/GradientColorMap.h"
49 
50 #include "CLI11.hpp"
51 
52 
53 using namespace std;
54 using namespace DGtal;
55 
56 
102 int main( int argc, char** argv )
103 {
104  typedef SpaceND<3,int> Space;
105  typedef KhalimskySpaceND<3,int> KSpace;
106  typedef HyperRectDomain<Space> Domain;
107  typedef ImageSelector<Domain, unsigned char>::Type Image;
108  typedef DigitalSetSelector< Domain, BIG_DS+HIGH_BEL_DS >::Type DigitalSet;
109  typedef SurfelAdjacency<KSpace::dimension> MySurfelAdjacency;
110 
111  // parse command line using CLI ----------------------------------------------
112  CLI::App app;
113  app.description("Display the boundary of a volume file by using QGLviewer. The mode specifies if you wish to see surface elements (BDRY), the inner voxels (INNER) or the outer voxels (OUTER) that touch the boundary. \n \t Example: 3dVolBoundaryViewer $DGtal/examples/samples/lobster.vol -m 60");
114  std::string inputFileName;
115  DGtal::int64_t rescaleInputMin {0};
116  DGtal::int64_t rescaleInputMax {255};
117  int dicomMin {-1000};
118  int dicomMax {3000};
119  int thresholdMin {0};
120  int thresholdMax {255};
121  std::string mode {"INNER"};
122 
123  app.add_option("-i,--input,1", inputFileName, "vol file (.vol, .longvol .p3d, .pgm3d and if WITH_ITK is selected: dicom, dcm, mha, mhd). For longvol, dicom, dcm, mha or mhd formats, the input values are linearly scaled between 0 and 255." )
124  ->required()
125  ->check(CLI::ExistingFile);
126 
127  app.add_option("--thresholdMin,-m", thresholdMin, "threshold min (excluded) to define binary shape.", true);
128  app.add_option("--thresholdMax,-M", thresholdMax, "threshold max (included) to define binary shape.", true);
129 #ifdef WITH_ITK
130  app.add_option("--dicomMin",dicomMin,"set minimum density threshold on Hounsfield scale", true );
131  app.add_option("--dicomMax",dicomMin,"set maximum density threshold on Hounsfield scale", true );
132 #endif
133  app.add_option("--mode", mode,"set mode for display: INNER: inner voxels, OUTER: outer voxels, BDRY: surfels", true )
134  -> check(CLI::IsMember({"INNER", "OUTER", "BDRY"}));
135 
136 
137  app.get_formatter()->column_width(40);
138  CLI11_PARSE(app, argc, argv);
139  // END parse command line using CLI ----------------------------------------------
140 
141 
142  QApplication application(argc,argv);
143 
144  string extension = inputFileName.substr(inputFileName.find_last_of(".") + 1);
145  if(extension!="vol" && extension != "p3d" && extension != "pgm3D" && extension != "pgm3d" && extension != "sdp" && extension != "pgm"
146 #ifdef WITH_ITK
147  && extension !="dcm"
148 #endif
149  ){
150  trace.info() << "File extension not recognized: "<< extension << std::endl;
151  return 0;
152  }
153 
154  if(extension=="vol" || extension=="pgm3d" || extension=="pgm3D"
155 #ifdef WITH_ITK
156  || extension =="dcm"
157 #endif
158  ){
159  trace.beginBlock( "Loading image into memory." );
160 #ifdef WITH_ITK
161  typedef DGtal::functors::Rescaling<int ,unsigned char > RescalFCT;
162  Image image = extension == "dcm" ? DicomReader< Image, RescalFCT >::importDicom( inputFileName,
163  RescalFCT(dicomMin,
164  dicomMax,
165  0, 255) ) :
166  GenericReader<Image>::import( inputFileName );
167 #else
168  Image image = GenericReader<Image>::import (inputFileName );
169 #endif
170  trace.info() << "Image loaded: "<<image<< std::endl;
171  trace.endBlock();
172 
174  trace.beginBlock( "Construct the Khalimsky space from the image domain." );
175  Domain domain = image.domain();
176  KSpace ks;
177  bool space_ok = ks.init( domain.lowerBound(), domain.upperBound(), true );
178  if (!space_ok)
179  {
180  trace.error() << "Error in the Khamisky space construction."<<std::endl;
181  return 2;
182  }
183  trace.endBlock();
185 
187  trace.beginBlock( "Wrapping a digital set around image. " );
188  typedef functors::IntervalForegroundPredicate<Image> ThresholdedImage;
189  ThresholdedImage thresholdedImage( image, thresholdMin, thresholdMax );
190  trace.endBlock();
192 
194  trace.beginBlock( "Extracting boundary by scanning the space. " );
195  typedef KSpace::SurfelSet SurfelSet;
196  typedef SetOfSurfels< KSpace, SurfelSet > MySetOfSurfels;
197  typedef DigitalSurface< MySetOfSurfels > MyDigitalSurface;
198  MySurfelAdjacency surfAdj( true ); // interior in all directions.
199  MySetOfSurfels theSetOfSurfels( ks, surfAdj );
200  Surfaces<KSpace>::sMakeBoundary( theSetOfSurfels.surfelSet(),
201  ks, thresholdedImage,
202  domain.lowerBound(),
203  domain.upperBound() );
204  MyDigitalSurface digSurf( theSetOfSurfels );
205  trace.info() << "Digital surface has " << digSurf.size() << " surfels."
206  << std::endl;
207  trace.endBlock();
209 
211  trace.beginBlock( "Displaying everything. " );
212  Viewer3D<Space,KSpace> viewer(ks);
213  viewer.setWindowTitle("Simple boundary of volume Viewer");
214  viewer.show();
215  typedef MyDigitalSurface::ConstIterator ConstIterator;
216  if ( mode == "BDRY" )
217  {
218  viewer << SetMode3D(ks.unsigns( *(digSurf.begin()) ).className(), "Basic");
219  for ( ConstIterator it = digSurf.begin(), itE = digSurf.end(); it != itE; ++it )
220  viewer << ks.unsigns( *it );
221  }else if ( mode == "INNER" )
222  for ( ConstIterator it = digSurf.begin(), itE = digSurf.end(); it != itE; ++it )
223  viewer << ks.sCoords( ks.sDirectIncident( *it, ks.sOrthDir( *it ) ) );
224  else if ( mode == "OUTER" )
225  for ( ConstIterator it = digSurf.begin(), itE = digSurf.end(); it != itE; ++it )
226  viewer << ks.sCoords( ks.sIndirectIncident( *it, ks.sOrthDir( *it ) ) );
227  else{
228  trace.error() << "Warning display mode (" << mode << ") not implemented." << std::endl;
229  trace.error() << "The display will be empty." << std::endl;
230  }
231  viewer << Viewer3D<>::updateDisplay;
232  trace.endBlock();
233  return application.exec();
234  }
235  return 0;
236 }
Definition: ATu0v1.h:57