33 #include "DGtal/base/Common.h" 34 #include "DGtal/helpers/StdDefs.h" 35 #include "DGtal/topology/helpers/Surfaces.h" 37 #include "DGtal/images/ImageContainerBySTLVector.h" 38 #include "DGtal/io/readers/GenericReader.h" 39 #include "DGtal/io/writers/GenericWriter.h" 40 #include "DGtal/images/IntervalForegroundPredicate.h" 41 #include <DGtal/topology/SetOfSurfels.h> 42 #include "DGtal/topology/DigitalSurface.h" 43 #include "DGtal/topology/helpers/Surfaces.h" 46 #include <boost/program_options/options_description.hpp> 47 #include <boost/program_options/parsers.hpp> 48 #include <boost/program_options/variables_map.hpp> 51 using namespace DGtal;
55 namespace po = boost::program_options;
110 typedef ImageContainerBySTLVector < Z3i::Domain, unsigned char > Image3D;
112 std::vector<unsigned int> getHistoFromImage(
const Image3D &image){
113 const Image3D::Domain &imgDom = image.domain();
114 std::vector<unsigned int> vectHisto(UCHAR_MAX);
115 for(Image3D::Domain::ConstIterator it=imgDom.begin(); it!= imgDom.end(); ++it){
116 vectHisto[image(*it)]++;
124 getOtsuThreshold(
const Image3D &image){
125 std::vector<unsigned int> histo = getHistoFromImage(image);
126 unsigned int imageSize = image.domain().size();
127 unsigned int sumA = 0;
128 unsigned int sumB = imageSize;
131 unsigned int sumMuAll= 0;
132 for(
unsigned int t=0; t< histo.size();t++){
133 sumMuAll+=histo[t]*t;
136 unsigned int thresholdRes=0;
138 for(
unsigned int t=0; t< histo.size(); t++){
149 double muAr=muA/(double)sumA;
150 double muBr=muB/(double)sumB;
151 double sigma= (double)sumA*(
double)sumB*(muAr-muBr)*(muAr-muBr);
162 int main(
int argc,
char** argv )
165 typedef Z3i::KSpace::SurfelSet SurfelSet;
166 typedef SetOfSurfels< Z3i::KSpace, SurfelSet > MySetOfSurfels;
167 typedef DigitalSurface< MySetOfSurfels > MyDigitalSurface;
170 po::options_description general_opt(
" \n Allowed options are ");
171 general_opt.add_options()
172 (
"help,h",
"display this message")
173 (
"input,i", po::value<std::string>(),
"volumetric input file (.vol, .pgm, .pgm3d, .longvol) " )
174 (
"output,o", po::value<std::string>(),
"volumetric output file (.vol, .pgm, .pgm3d, .longvol) " )
175 (
"labelBackground",
"option to define a label to regions associated to object background. ")
176 (
"thresholdMin,m", po::value<int>()->default_value(0),
"min threshold (if not given the max threshold is computed with Otsu algorithm)." )
177 (
"thresholdMax,M", po::value<int>(),
"max threshold (default 255)" );
181 po::variables_map vm;
183 po::store(po::parse_command_line(argc, argv, general_opt), vm);
184 }
catch(
const std::exception& ex){
186 trace.info()<<
"Error checking program options: "<< ex.what()<< endl;
189 if( !parseOK || vm.count(
"help")||argc<=1)
191 std::cout <<
"Usage: " << argv[0] <<
" [input] [output]\n" 192 <<
"Segment volumetric file from a simple threshold which can be set automatically from the otsu estimation.\n" 193 <<
"The segmentation result is given by an integer label given in the resulting image." 194 << general_opt <<
"\n";
195 std::cout <<
"Example:\n" 196 <<
"volSegment -i ${DGtal}/examples/samples/lobster.vol -o segmentation.vol \n";
201 if(! vm.count(
"input") ||! vm.count(
"output"))
203 trace.error() <<
" Input and output filename are needed to be defined" << endl;
208 string inputFilename = vm[
"input"].as<std::string>();
209 string outputFilename = vm[
"output"].as<std::string>();
211 trace.info() <<
"Reading input file " << inputFilename ;
212 Image3D inputImage = DGtal::GenericReader<Image3D>::import(inputFilename);
213 Image3D imageResuSegmentation(inputImage.domain());
215 trace.info() <<
" [done] " << std::endl ;
216 std::ofstream outStream;
217 outStream.open(outputFilename.c_str());
218 int minTh = vm[
"thresholdMin"].as<
int>();
220 if(!vm.count(
"thresholdMax")){
221 maxTh = getOtsuThreshold(inputImage);
222 trace.info() <<
"maximal threshold value not specified, using Otsu value: " << maxTh << std::endl;
224 maxTh = vm[
"thresholdMax"].as<
int>();
226 trace.info() <<
"Processing image to output file " << outputFilename << std::endl;
228 functors::IntervalForegroundPredicate<Image3D> simplePredicate ( inputImage, minTh, maxTh );
229 SurfelAdjacency< Z3i::KSpace::dimension > SAdj (
true );
231 bool space_ok = K.init( inputImage.domain().lowerBound(), inputImage.domain().upperBound(), false );
233 trace.error() <<
"problem initializing 3d space" << endl;
237 std::vector< std::vector<Z3i::SCell > > vectConnectedSCell;
238 Surfaces<Z3i::KSpace>::extractAllConnectedSCell(vectConnectedSCell,K, SAdj, simplePredicate,
false);
239 trace.progressBar(0, vectConnectedSCell.size());
240 for(
unsigned int i = 0; i<vectConnectedSCell.size(); i++)
242 trace.progressBar(i, vectConnectedSCell.size());
243 MySetOfSurfels aSet(K, SAdj);
244 Z3i::Point lowerPoint, upperPoint;
247 for(std::vector<Z3i::SCell>::const_iterator it= vectConnectedSCell.at(i).begin(); it != vectConnectedSCell.at(i).end(); ++it)
249 aSet.surfelSet().insert(aSet.surfelSet().begin(), *it);
250 unsigned int orth_dir = K.sOrthDir( *it );
251 p1 = K.sCoords( K.sIncident( *it, orth_dir,
true ) );
252 p2 = K.sCoords( K.sIncident( *it, orth_dir,
false ) );
253 if(p1[0] < lowerPoint[0]) lowerPoint[0]= p1[0];
254 if(p1[1] < lowerPoint[1]) lowerPoint[1]= p1[1];
255 if(p1[2] < lowerPoint[2]) lowerPoint[2]= p1[2];
257 if(p1[0] > upperPoint[0]) upperPoint[0]= p1[0];
258 if(p1[1] > upperPoint[1]) upperPoint[1]= p1[1];
259 if(p1[2] > upperPoint[2]) upperPoint[2]= p1[2];
261 if(p2[0] < lowerPoint[0]) lowerPoint[0]= p2[0];
262 if(p2[1] < lowerPoint[1]) lowerPoint[1]= p2[1];
263 if(p2[2] < lowerPoint[2]) lowerPoint[2]= p2[2];
265 if(p2[0] > upperPoint[0]) upperPoint[0]= p2[0];
266 if(p2[1] > upperPoint[1]) upperPoint[1]= p2[1];
267 if(p2[2] > upperPoint[2]) upperPoint[2]= p2[2];
272 kRestr.init( lowerPoint, upperPoint,
false );
273 if(simplePredicate(p2)){
274 DGtal::Surfaces<Z3i::KSpace>::uFillInterior( kRestr, aSet.surfelPredicate(),
275 imageResuSegmentation,
277 }
else if (vm.count(
"labelBackground")){
278 DGtal::Surfaces<Z3i::KSpace>::uFillExterior( kRestr, aSet.surfelPredicate(),
279 imageResuSegmentation,
283 trace.progressBar(vectConnectedSCell.size(), vectConnectedSCell.size());
284 trace.info() << std::endl;
285 GenericWriter<Image3D>::exportFile(outputFilename, imageResuSegmentation);