32 #include <boost/program_options/options_description.hpp> 33 #include <boost/program_options/parsers.hpp> 34 #include <boost/program_options/variables_map.hpp> 36 #include "DGtal/base/Common.h" 37 #include "DGtal/base/BasicFunctors.h" 38 #include "DGtal/kernel/SpaceND.h" 39 #include "DGtal/kernel/domains/HyperRectDomain.h" 40 #include "DGtal/images/ImageSelector.h" 41 #include "DGtal/images/IntervalForegroundPredicate.h" 42 #include "DGtal/topology/KhalimskySpaceND.h" 43 #include "DGtal/topology/DigitalSurface.h" 44 #include "DGtal/topology/SetOfSurfels.h" 45 #include "DGtal/io/boards/Board3D.h" 46 #include "DGtal/io/readers/PointListReader.h" 47 #include "DGtal/io/readers/GenericReader.h" 49 #include "DGtal/io/readers/DicomReader.h" 51 #include "DGtal/io/Color.h" 52 #include "DGtal/io/colormaps/GradientColorMap.h" 56 using namespace DGtal;
60 namespace po = boost::program_options;
104 int main(
int argc,
char** argv )
106 typedef SpaceND<3,int> Space;
107 typedef KhalimskySpaceND<3,int> KSpace;
108 typedef HyperRectDomain<Space> Domain;
109 typedef ImageSelector<Domain, unsigned char>::Type Image;
110 typedef DigitalSetSelector< Domain, BIG_DS+HIGH_BEL_DS >::Type DigitalSet;
111 typedef SurfelAdjacency<KSpace::dimension> MySurfelAdjacency;
114 po::options_description general_opt(
"Allowed options are: ");
115 general_opt.add_options()
116 (
"help,h",
"display this message")
117 (
"input,i", po::value<std::string>(),
"vol file (.vol) , pgm3d (.p3d or .pgm3d, pgm (with 3 dims)) file or sdp (sequence of discrete points)" )
118 (
"output,o", po::value<std::string>(),
"output obj file (.obj)" )
119 (
"thresholdMin,m", po::value<int>()->default_value(0),
"threshold min (excluded) to define binary shape" )
120 (
"thresholdMax,M", po::value<int>()->default_value(255),
"threshold max (included) to define binary shape" )
122 (
"dicomMin", po::value<int>()->default_value(-1000),
"set minimum density threshold on Hounsfield scale")
123 (
"dicomMax", po::value<int>()->default_value(3000),
"set maximum density threshold on Hounsfield scale")
125 (
"mode", po::value<std::string>()->default_value(
"INNER"),
"set mode for display: INNER: inner voxels, OUTER: outer voxels, BDRY: surfels")
126 (
"normalization,n",
"Normalization so that the geometry fits in [-1/2,1/2]^3") ;
129 po::variables_map vm;
131 po::store(po::parse_command_line(argc, argv, general_opt), vm);
132 }
catch(
const std::exception& ex){
134 trace.info()<<
"Error checking program options: "<< ex.what()<< endl;
137 if( !parseOK || vm.count(
"help")||argc<=1)
139 std::cout <<
"Usage: " << argv[0] <<
" -i [input] -o [output]\n" 140 <<
"Export the boundary of a volume file to OBJ format. 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
141 << general_opt <<
"\n";
145 if(! vm.count(
"input"))
147 trace.error() <<
" The file name was defined" << endl;
151 if(! vm.count(
"output"))
153 trace.error() <<
" The output filename was defined" << endl;
157 string inputFilename = vm[
"input"].as<std::string>();
158 int thresholdMin = vm[
"thresholdMin"].as<
int>();
159 int thresholdMax = vm[
"thresholdMax"].as<
int>();
160 string mode = vm[
"mode"].as<
string>();
161 bool normalization =
false;
162 if (vm.count(
"normalization"))
163 normalization =
true;
165 string extension = inputFilename.substr(inputFilename.find_last_of(
".") + 1);
166 if(extension!=
"vol" && extension !=
"p3d" && extension !=
"pgm3D" && extension !=
"pgm3d" && extension !=
"sdp" && extension !=
"pgm" 171 trace.info() <<
"File extension not recognized: "<< extension << std::endl;
175 if(extension==
"vol" || extension==
"pgm3d" || extension==
"pgm3D" 180 trace.beginBlock(
"Loading image into memory." );
182 int dicomMin = vm[
"dicomMin"].as<
int>();
183 int dicomMax = vm[
"dicomMax"].as<
int>();
184 typedef DGtal::functors::Rescaling<int ,unsigned char > RescalFCT;
185 Image image = extension ==
"dcm" ? DicomReader< Image, RescalFCT >::importDicom( inputFilename,
189 GenericReader<Image>::import( inputFilename );
191 Image image = GenericReader<Image>::import (inputFilename );
193 trace.info() <<
"Image loaded: "<<image<< std::endl;
196 trace.beginBlock(
"Construct the Khalimsky space from the image domain." );
197 Domain domain = image.domain();
199 bool space_ok = ks.init( domain.lowerBound(), domain.upperBound(), true );
202 trace.error() <<
"Error in the Khamisky space construction."<<std::endl;
207 trace.beginBlock(
"Wrapping a digital set around image. " );
208 typedef functors::IntervalForegroundPredicate<Image> ThresholdedImage;
209 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." 227 trace.beginBlock(
"Displaying everything. " );
228 Board3D<Space,KSpace> board(ks);
230 board << SetMode3D( ks.unsigns( *digSurf.begin() ).className(),
"Basic" );
232 typedef MyDigitalSurface::ConstIterator ConstIterator;
233 if ( mode ==
"BDRY" )
234 for ( ConstIterator it = digSurf.begin(), itE = digSurf.end(); it != itE; ++it )
235 board << ks.unsigns( *it );
236 else if ( mode ==
"INNER" )
237 for ( ConstIterator it = digSurf.begin(), itE = digSurf.end(); it != itE; ++it )
238 board << ks.sCoords( ks.sDirectIncident( *it, ks.sOrthDir( *it ) ) );
239 else if ( mode ==
"OUTER" )
240 for ( ConstIterator it = digSurf.begin(), itE = digSurf.end(); it != itE; ++it )
241 board << ks.sCoords( ks.sIndirectIncident( *it, ks.sOrthDir( *it ) ) );
243 string outputFilename = vm[
"output"].as<std::string>();
245 board.saveOBJ(outputFilename, normalization);