DGtalTools  0.9.4
3dVolBoundaryViewer.cpp
1 
29 #include <iostream>
31 #include <boost/program_options/options_description.hpp>
32 #include <boost/program_options/parsers.hpp>
33 #include <boost/program_options/variables_map.hpp>
34 
35 #include "DGtal/base/Common.h"
36 #include "DGtal/base/BasicFunctors.h"
37 #include "DGtal/kernel/SpaceND.h"
38 #include "DGtal/kernel/domains/HyperRectDomain.h"
39 #include "DGtal/images/ImageSelector.h"
40 #include "DGtal/images/IntervalForegroundPredicate.h"
41 #include "DGtal/topology/KhalimskySpaceND.h"
42 #include "DGtal/topology/DigitalSurface.h"
43 #include "DGtal/topology/SetOfSurfels.h"
44 #include "DGtal/io/viewers/Viewer3D.h"
45 #include "DGtal/io/readers/PointListReader.h"
46 #include "DGtal/io/readers/GenericReader.h"
47 #ifdef WITH_ITK
48 #include "DGtal/io/readers/DicomReader.h"
49 #endif
50 #include "DGtal/io/Color.h"
51 #include "DGtal/io/colormaps/GradientColorMap.h"
52 
53 
54 using namespace std;
55 using namespace DGtal;
56 
57 
106 namespace po = boost::program_options;
108 
109 int main( int argc, char** argv )
110 {
111  typedef SpaceND<3,int> Space;
116  typedef SurfelAdjacency<KSpace::dimension> MySurfelAdjacency;
117 
118  // parse command line ----------------------------------------------
119  po::options_description general_opt("Allowed options are: ");
120  general_opt.add_options()
121  ("help,h", "display this message")
122  ("input,i", po::value<std::string>(), "vol file (.vol) , pgm3d (.p3d or .pgm3d, pgm (with 3 dims)) file or sdp (sequence of discrete points)" )
123  ("thresholdMin,m", po::value<int>()->default_value(0), "threshold min (excluded) to define binary shape" )
124  ("thresholdMax,M", po::value<int>()->default_value(255), "threshold max (included) to define binary shape" )
125 #ifdef WITH_ITK
126  ("dicomMin", po::value<int>()->default_value(-1000), "set minimum density threshold on Hounsfield scale")
127  ("dicomMax", po::value<int>()->default_value(3000), "set maximum density threshold on Hounsfield scale")
128 #endif
129  ("mode", po::value<std::string>()->default_value("INNER"), "set mode for display: INNER: inner voxels, OUTER: outer voxels, BDRY: surfels") ;
130 
131  bool parseOK=true;
132  po::variables_map vm;
133  try{
134  po::store(po::parse_command_line(argc, argv, general_opt), vm);
135  }catch(const std::exception& ex){
136  parseOK=false;
137  trace.info()<< "Error checking program options: "<< ex.what()<< endl;
138  }
139  po::notify(vm);
140  if( !parseOK || vm.count("help")||argc<=1)
141  {
142  std::cout << "Usage: " << argv[0] << " -i [input]\n"
143  << "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."<< endl
144  << general_opt << "\n";
145  return 0;
146  }
147 
148  if(! vm.count("input"))
149  {
150  trace.error() << " The file name was defined" << endl;
151  return 0;
152  }
153  string inputFilename = vm["input"].as<std::string>();
154  int thresholdMin = vm["thresholdMin"].as<int>();
155  int thresholdMax = vm["thresholdMax"].as<int>();
156  string mode = vm["mode"].as<string>();
157 
158  QApplication application(argc,argv);
159 
160  string extension = inputFilename.substr(inputFilename.find_last_of(".") + 1);
161  if(extension!="vol" && extension != "p3d" && extension != "pgm3D" && extension != "pgm3d" && extension != "sdp" && extension != "pgm"
162 #ifdef WITH_ITK
163  && extension !="dcm"
164 #endif
165  ){
166  trace.info() << "File extension not recognized: "<< extension << std::endl;
167  return 0;
168  }
169 
170  if(extension=="vol" || extension=="pgm3d" || extension=="pgm3D"
171 #ifdef WITH_ITK
172  || extension =="dcm"
173 #endif
174  ){
175  trace.beginBlock( "Loading image into memory." );
176 #ifdef WITH_ITK
177  int dicomMin = vm["dicomMin"].as<int>();
178  int dicomMax = vm["dicomMax"].as<int>();
180  Image image = extension == "dcm" ? DicomReader< Image, RescalFCT >::importDicom( inputFilename,
181  RescalFCT(dicomMin,
182  dicomMax,
183  0, 255) ) :
184  GenericReader<Image>::import( inputFilename );
185 #else
186  Image image = GenericReader<Image>::import (inputFilename );
187 #endif
188  trace.info() << "Image loaded: "<<image<< std::endl;
189  trace.endBlock();
190 
192  trace.beginBlock( "Construct the Khalimsky space from the image domain." );
193  Domain domain = image.domain();
194  KSpace ks;
195  bool space_ok = ks.init( domain.lowerBound(), domain.upperBound(), true );
196  if (!space_ok)
197  {
198  trace.error() << "Error in the Khamisky space construction."<<std::endl;
199  return 2;
200  }
201  trace.endBlock();
203 
205  trace.beginBlock( "Wrapping a digital set around image. " );
206  typedef functors::IntervalForegroundPredicate<Image> ThresholdedImage;
207  ThresholdedImage thresholdedImage( image, thresholdMin, thresholdMax );
208  trace.endBlock();
210 
212  trace.beginBlock( "Extracting boundary by scanning the space. " );
213  typedef KSpace::SurfelSet SurfelSet;
214  typedef SetOfSurfels< KSpace, SurfelSet > MySetOfSurfels;
215  typedef DigitalSurface< MySetOfSurfels > MyDigitalSurface;
216  MySurfelAdjacency surfAdj( true ); // interior in all directions.
217  MySetOfSurfels theSetOfSurfels( ks, surfAdj );
218  Surfaces<KSpace>::sMakeBoundary( theSetOfSurfels.surfelSet(),
219  ks, thresholdedImage,
220  domain.lowerBound(),
221  domain.upperBound() );
222  MyDigitalSurface digSurf( theSetOfSurfels );
223  trace.info() << "Digital surface has " << digSurf.size() << " surfels."
224  << std::endl;
225  trace.endBlock();
227 
229  trace.beginBlock( "Displaying everything. " );
230  Viewer3D<Space,KSpace> viewer(ks);
231  viewer.setWindowTitle("Simple boundary of volume Viewer");
232  viewer.show();
233  typedef MyDigitalSurface::ConstIterator ConstIterator;
234  if ( mode == "BDRY" ){
235  viewer << SetMode3D(ks.unsigns( *(digSurf.begin()) ).className(), "Basic");
236  for ( ConstIterator it = digSurf.begin(), itE = digSurf.end(); it != itE; ++it )
237  viewer << ks.unsigns( *it );
238  }else if ( mode == "INNER" )
239  for ( ConstIterator it = digSurf.begin(), itE = digSurf.end(); it != itE; ++it )
240  viewer << ks.sCoords( ks.sDirectIncident( *it, ks.sOrthDir( *it ) ) );
241  else if ( mode == "OUTER" )
242  for ( ConstIterator it = digSurf.begin(), itE = digSurf.end(); it != itE; ++it )
243  viewer << ks.sCoords( ks.sIndirectIncident( *it, ks.sOrthDir( *it ) ) );
244  else{
245  trace.error() << "Warning display mode (" << mode << ") not implemented." << std::endl;
246  trace.error() << "The display will be empty." << std::endl;
247  }
248  viewer << Viewer3D<>::updateDisplay;
249  trace.endBlock();
250  return application.exec();
251  }
252  return 0;
253 }
void beginBlock(const std::string &keyword="")
DigitalSetSelector< Domain, BIG_DS+HIGH_BEL_DS >::Type DigitalSet
SpaceND< 2, Integer > Space
std::set< SCell > SurfelSet
STL namespace.
double endBlock()
std::string className() const
SCell sDirectIncident(const SCell &p, Dimension k) const
const Point & lowerBound() const
Dimension sOrthDir(const SCell &s) const
bool init(const Point &lower, const Point &upper, bool isClosed)
Point sCoords(const SCell &c) const
Trace trace(traceWriterTerm)
std::ostream & info()
SCell sIndirectIncident(const SCell &p, Dimension k) const
const Domain & domain() const
Cell unsigns(const SCell &p) const
std::ostream & error()
typename Self::Domain Domain