31 #include <boost/program_options/options_description.hpp> 32 #include <boost/program_options/parsers.hpp> 33 #include <boost/program_options/variables_map.hpp> 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" 48 #include "DGtal/io/readers/DicomReader.h" 50 #include "DGtal/io/Color.h" 51 #include "DGtal/io/colormaps/GradientColorMap.h" 55 using namespace DGtal;
106 namespace po = boost::program_options;
109 int main(
int argc,
char** argv )
111 typedef SpaceND<3,int> Space;
112 typedef KhalimskySpaceND<3,int> KSpace;
113 typedef HyperRectDomain<Space> Domain;
114 typedef ImageSelector<Domain, unsigned char>::Type Image;
115 typedef DigitalSetSelector< Domain, BIG_DS+HIGH_BEL_DS >::Type DigitalSet;
116 typedef SurfelAdjacency<KSpace::dimension> MySurfelAdjacency;
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" )
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")
129 (
"mode", po::value<std::string>()->default_value(
"INNER"),
"set mode for display: INNER: inner voxels, OUTER: outer voxels, BDRY: surfels") ;
132 po::variables_map vm;
134 po::store(po::parse_command_line(argc, argv, general_opt), vm);
135 }
catch(
const std::exception& ex){
137 trace.info()<<
"Error checking program options: "<< ex.what()<< endl;
140 if( !parseOK || vm.count(
"help")||argc<=1)
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";
148 if(! vm.count(
"input"))
150 trace.error() <<
" The file name was defined" << endl;
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>();
158 QApplication application(argc,argv);
160 string extension = inputFilename.substr(inputFilename.find_last_of(
".") + 1);
161 if(extension!=
"vol" && extension !=
"p3d" && extension !=
"pgm3D" && extension !=
"pgm3d" && extension !=
"sdp" && extension !=
"pgm" 166 trace.info() <<
"File extension not recognized: "<< extension << std::endl;
170 if(extension==
"vol" || extension==
"pgm3d" || extension==
"pgm3D" 175 trace.beginBlock(
"Loading image into memory." );
177 int dicomMin = vm[
"dicomMin"].as<
int>();
178 int dicomMax = vm[
"dicomMax"].as<
int>();
179 typedef DGtal::functors::Rescaling<int ,unsigned char > RescalFCT;
180 Image image = extension ==
"dcm" ? DicomReader< Image, RescalFCT >::importDicom( inputFilename,
184 GenericReader<Image>::import( inputFilename );
186 Image image = GenericReader<Image>::import (inputFilename );
188 trace.info() <<
"Image loaded: "<<image<< std::endl;
192 trace.beginBlock(
"Construct the Khalimsky space from the image domain." );
193 Domain domain = image.domain();
195 bool space_ok = ks.init( domain.lowerBound(), domain.upperBound(), true );
198 trace.error() <<
"Error in the Khamisky space construction."<<std::endl;
205 trace.beginBlock(
"Wrapping a digital set around image. " );
206 typedef functors::IntervalForegroundPredicate<Image> ThresholdedImage;
207 ThresholdedImage thresholdedImage( image, thresholdMin, thresholdMax );
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 );
217 MySetOfSurfels theSetOfSurfels( ks, surfAdj );
218 Surfaces<KSpace>::sMakeBoundary( theSetOfSurfels.surfelSet(),
219 ks, thresholdedImage,
221 domain.upperBound() );
222 MyDigitalSurface digSurf( theSetOfSurfels );
223 trace.info() <<
"Digital surface has " << digSurf.size() <<
" surfels." 229 trace.beginBlock(
"Displaying everything. " );
230 Viewer3D<Space,KSpace> viewer(ks);
231 viewer.setWindowTitle(
"Simple boundary of volume Viewer");
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 ) ) );
245 trace.error() <<
"Warning display mode (" << mode <<
") not implemented." << std::endl;
246 trace.error() <<
"The display will be empty." << std::endl;
248 viewer << Viewer3D<>::updateDisplay;
250 return application.exec();