DGtalTools  1.2.0
volSegment.cpp
1 
29 #include <iostream>
30 #include <fstream>
31 
32 #include "DGtal/base/Common.h"
33 #include "DGtal/helpers/StdDefs.h"
34 #include "DGtal/topology/helpers/Surfaces.h"
35 
36 #include "DGtal/images/ImageContainerBySTLVector.h"
37 #include "DGtal/io/readers/GenericReader.h"
38 #include "DGtal/io/writers/GenericWriter.h"
39 #include "DGtal/images/IntervalForegroundPredicate.h"
40 #include <DGtal/topology/SetOfSurfels.h>
41 #include "DGtal/topology/DigitalSurface.h"
42 #include "DGtal/topology/helpers/Surfaces.h"
43 
44 #include "CLI11.hpp"
45 
46 
47 using namespace std;
48 using namespace DGtal;
49 
50 
103 typedef ImageContainerBySTLVector < Z3i::Domain, unsigned char > Image3D;
104 
105 std::vector<unsigned int> getHistoFromImage(const Image3D &image){
106  const Image3D::Domain &imgDom = image.domain();
107  std::vector<unsigned int> vectHisto(UCHAR_MAX);
108  for(Image3D::Domain::ConstIterator it=imgDom.begin(); it!= imgDom.end(); ++it){
109  vectHisto[image(*it)]++;
110  }
111  return vectHisto;
112 }
113 
114 
115 unsigned int
116 getOtsuThreshold(const Image3D &image){
117  std::vector<unsigned int> histo = getHistoFromImage(image);
118  unsigned int imageSize = image.domain().size();
119  unsigned int sumA = 0;
120  unsigned int sumB = imageSize;
121  unsigned int muA=0;
122  unsigned int muB=0;
123  unsigned int sumMuAll= 0;
124  for( unsigned int t=0; t< histo.size();t++){
125  sumMuAll+=histo[t]*t;
126  }
127 
128  unsigned int thresholdRes=0;
129  double valMax=0.0;
130  for( unsigned int t=0; t< histo.size(); t++){
131  sumA+=histo[t];
132  if(sumA==0)
133  continue;
134  sumB=imageSize-sumA;
135  if(sumB==0){
136  break;
137  }
138 
139  muA+=histo[t]*t;
140  muB=sumMuAll-muA;
141  double muAr=muA/(double)sumA;
142  double muBr=muB/(double)sumB;
143  double sigma= (double)sumA*(double)sumB*(muAr-muBr)*(muAr-muBr);
144  if(valMax<=sigma){
145  valMax=sigma;
146  thresholdRes=t;
147  }
148  }
149  return thresholdRes;
150 }
151 
152 
153 
154 int main( int argc, char** argv )
155 {
156  typedef Z3i::KSpace::SurfelSet SurfelSet;
157  typedef SetOfSurfels< Z3i::KSpace, SurfelSet > MySetOfSurfels;
158 
159  // parse command line using CLI ----------------------------------------------
160  CLI::App app;
161  std::string inputFileName;
162  std::string outputFileName {"result.vol"};
163  bool labelBackground {false};
164  int minTh {0};
165  int maxTh {255};
166 
167  app.description("Segment volumetric file from a simple threshold which can be set automatically from the otsu estimation.\n The segmentation result is given by an integer label given in the resulting image. Example:\n volSegment ${DGtal}/examples/samples/lobster.vol segmentation.vol \n");
168  app.add_option("-i,--input,1", inputFileName, "volumetric input file (.vol, .pgm, .pgm3d, .longvol)." )
169  ->required()
170  ->check(CLI::ExistingFile);
171  app.add_option("-o,--output,2", outputFileName, "volumetric output file (.vol, .pgm, .pgm3d, .longvol)", true);
172  auto labelOpt = app.add_flag("--labelBackground",labelBackground, "option to define a label to regions associated to object background.");
173  app.add_option("-m,--thresholdMin",minTh, "min threshold (if not given the max threshold is computed with Otsu algorithm).", true );
174  auto maxThOpt = app.add_option("-M,--thresholdMax", maxTh, "max threshold", true );
175 
176 
177  app.get_formatter()->column_width(40);
178  CLI11_PARSE(app, argc, argv);
179  // END parse command line using CLI ----------------------------------------------
180 
181  trace.info() << "Reading input file " << inputFileName ;
182  Image3D inputImage = DGtal::GenericReader<Image3D>::import(inputFileName);
183  Image3D imageResuSegmentation(inputImage.domain());
184 
185  trace.info() << " [done] " << std::endl ;
186  std::ofstream outStream;
187  outStream.open(outputFileName.c_str());
188  if(maxThOpt->count()==0){
189  maxTh = getOtsuThreshold(inputImage);
190  trace.info() << "maximal threshold value not specified, using Otsu value: " << maxTh << std::endl;
191  }
192  trace.info() << "Processing image to output file " << outputFileName << std::endl;
193 
194  functors::IntervalForegroundPredicate<Image3D> simplePredicate ( inputImage, minTh, maxTh );
195  SurfelAdjacency< Z3i::KSpace::dimension > SAdj ( true );
196  Z3i::KSpace K;
197  bool space_ok = K.init( inputImage.domain().lowerBound(), inputImage.domain().upperBound(), false );
198  if(!space_ok){
199  trace.error() << "problem initializing 3d space" << endl;
200  }
201 
202  std::vector< std::vector<Z3i::SCell > > vectConnectedSCell;
203  Surfaces<Z3i::KSpace>::extractAllConnectedSCell(vectConnectedSCell,K, SAdj, simplePredicate, false);
204  trace.progressBar(0, vectConnectedSCell.size());
205  for(unsigned int i = 0; i<vectConnectedSCell.size(); i++)
206  {
207  trace.progressBar(i, vectConnectedSCell.size());
208  MySetOfSurfels aSet(K, SAdj);
209  Z3i::Point lowerPoint, upperPoint;
210  Z3i::Point p1;
211  Z3i::Point p2;
212  for(std::vector<Z3i::SCell>::const_iterator it= vectConnectedSCell.at(i).begin(); it != vectConnectedSCell.at(i).end(); ++it)
213  {
214  aSet.surfelSet().insert(aSet.surfelSet().begin(), *it);
215  unsigned int orth_dir = K.sOrthDir( *it );
216  p1 = K.sCoords( K.sIncident( *it, orth_dir, true ) );
217  p2 = K.sCoords( K.sIncident( *it, orth_dir, false ) );
218  if(p1[0] < lowerPoint[0]) lowerPoint[0]= p1[0];
219  if(p1[1] < lowerPoint[1]) lowerPoint[1]= p1[1];
220  if(p1[2] < lowerPoint[2]) lowerPoint[2]= p1[2];
221 
222  if(p1[0] > upperPoint[0]) upperPoint[0]= p1[0];
223  if(p1[1] > upperPoint[1]) upperPoint[1]= p1[1];
224  if(p1[2] > upperPoint[2]) upperPoint[2]= p1[2];
225 
226  if(p2[0] < lowerPoint[0]) lowerPoint[0]= p2[0];
227  if(p2[1] < lowerPoint[1]) lowerPoint[1]= p2[1];
228  if(p2[2] < lowerPoint[2]) lowerPoint[2]= p2[2];
229 
230  if(p2[0] > upperPoint[0]) upperPoint[0]= p2[0];
231  if(p2[1] > upperPoint[1]) upperPoint[1]= p2[1];
232  if(p2[2] > upperPoint[2]) upperPoint[2]= p2[2];
233 
234  }
235 
236  Z3i::KSpace kRestr ;
237  kRestr.init( lowerPoint, upperPoint, false );
238  if(simplePredicate(p2)){
239  DGtal::Surfaces<Z3i::KSpace>::uFillInterior( kRestr, aSet.surfelPredicate(),
240  imageResuSegmentation,
241  i, false, false);
242  }else if (labelOpt->count() > 0){
243  DGtal::Surfaces<Z3i::KSpace>::uFillExterior( kRestr, aSet.surfelPredicate(),
244  imageResuSegmentation,
245  i+1, false, false);
246  }
247  }
248  trace.progressBar(vectConnectedSCell.size(), vectConnectedSCell.size());
249  trace.info() << std::endl;
250  GenericWriter<Image3D>::exportFile(outputFileName, imageResuSegmentation);
251  return 0;
252 }
253 
Definition: ATu0v1.h:57