DGtalTools  1.2.0
img2freeman.cpp
1 #include "DGtal/io/colormaps/GrayscaleColorMap.h"
2 #include "DGtal/io/readers/GenericReader.h"
3 #include "DGtal/images/ImageContainerBySTLVector.h"
4 #include "DGtal/images/ImageSelector.h"
5 #include "DGtal/geometry/curves/FreemanChain.h"
6 #include "DGtal/geometry/helpers/ContourHelper.h"
7 #include "DGtal/topology/helpers/Surfaces.h"
8 
9 #include "CLI11.hpp"
10 
11 #include <vector>
12 #include <string>
13 #include <climits>
14 
15 
16 using namespace DGtal;
77 typedef ImageSelector < Z2i::Domain, unsigned char>::Type Image;
78 
79 
80 std::vector<unsigned int> getHistoFromImage(const Image &image){
81  const Image::Domain &imgDom = image.domain();
82  std::vector<unsigned int> vectHisto(UCHAR_MAX);
83  for(Image::Domain::ConstIterator it=imgDom.begin(); it!= imgDom.end(); ++it){
84  vectHisto[image(*it)]++;
85  }
86  return vectHisto;
87 }
88 
89 unsigned int
90 getOtsuThreshold(const Image &image){
91  std::vector<unsigned int> histo = getHistoFromImage(image);
92  unsigned int imageSize = image.domain().size();
93  unsigned int sumA = 0;
94  unsigned int sumB = imageSize;
95  unsigned int muA=0;
96  unsigned int muB=0;
97  unsigned int sumMuAll= 0;
98  for( unsigned int t=0; t< histo.size();t++){
99  sumMuAll+=histo[t]*t;
100  }
101 
102  unsigned int thresholdRes=0;
103  double valMax=0.0;
104  for( unsigned int t=0; t< histo.size(); t++){
105  sumA+=histo[t];
106  if(sumA==0)
107  continue;
108  sumB=imageSize-sumA;
109  if(sumB==0){
110  break;
111  }
112 
113  muA+=histo[t]*t;
114  muB=sumMuAll-muA;
115  double muAr=muA/(double)sumA;
116  double muBr=muB/(double)sumB;
117  double sigma= (double)sumA*(double)sumB*(muAr-muBr)*(muAr-muBr);
118  if(valMax<=sigma){
119  valMax=sigma;
120  thresholdRes=t;
121  }
122  }
123  return thresholdRes;
124 }
125 
126 struct CompContours{
127  bool operator()(std::vector<Z2i::Point> a, std::vector<Z2i::Point> b ){
128  return a.size() > b.size();
129  }
130 };
131 
132 
133 void saveAllContoursAsFc( std::vector< std::vector< Z2i::Point > > vectContoursBdryPointels,
134  unsigned int minSize, bool sort=false){
135  CompContours comp;
136  if(sort){
137  std::sort(vectContoursBdryPointels.begin(), vectContoursBdryPointels.end(), comp);
138  }
139  for(unsigned int k=0; k<vectContoursBdryPointels.size(); k++){
140  if(vectContoursBdryPointels.at(k).size()>minSize){
141  FreemanChain<Z2i::Integer> fc (vectContoursBdryPointels.at(k));
142  std::cout << fc.x0 << " " << fc.y0 << " " << fc.chain << std::endl;
143 
144  }
145  }
146 }
147 
148 
149 void saveSelContoursAsFC(std::vector< std::vector< Z2i::Point > > vectContoursBdryPointels,
150  unsigned int minSize, Z2i::Point refPoint, double selectDistanceMax,
151  bool sort=false){
152  CompContours comp;
153  if(sort){
154  std::sort(vectContoursBdryPointels.begin(), vectContoursBdryPointels.end(), comp);
155  }
156 
157  for(unsigned int k=0; k<vectContoursBdryPointels.size(); k++){
158  if(vectContoursBdryPointels.at(k).size()>minSize){
159  Z2i::RealPoint ptMean = ContourHelper::getBarycenter(vectContoursBdryPointels.at(k));
160  unsigned int distance = (unsigned int)ceil(sqrt((ptMean[0]-refPoint[0])*(ptMean[0]-refPoint[0])+
161  (ptMean[1]-refPoint[1])*(ptMean[1]-refPoint[1])));
162  if(distance<=selectDistanceMax){
163  FreemanChain<Z2i::Integer> fc (vectContoursBdryPointels.at(k));
164  std::cout << fc.x0 << " " << fc.y0 << " " << fc.chain << std::endl;
165  }
166  }
167  }
168 }
169 
170 
171 
172 
173 int main( int argc, char** argv )
174 {
175  // parse command line using CLI ----------------------------------------------
176  CLI::App app;
177  std::string inputFileName;
178  std::string outputFileName {"result.fc"};
179 
180  double minThreshold {128};
181  double maxThreshold {255};
182  unsigned int minSize {0};
183  bool select {false};
184  bool sortCnt {false};
185 
186  Z2i::Point selectCenter;
187  unsigned int selectDistanceMax = 0;
188  std::vector<int> cntConstraints;
189  std::vector<int> vectRangeMin, vectRangeMax, vectRange;
190 
191  app.description("Extract FreemanChains from thresholded image.\n Basic example: \t img2freeman [options] --input <imageName> -min 128 -max 255 > contours.fc \n Note that if you don't specify any threshold a threshold threshold max is automatically defined from the Otsu algorithm with min=0. ");
192  app.add_option("-i,--input,1", inputFileName, "input image file name (any 2D image format accepted by DGtal::GenericReader)." )
193  ->required()
194  ->check(CLI::ExistingFile);
195  app.add_option("-m,--min", minThreshold, "min image threshold value (default 128)");
196  app.add_option("-M,--max", maxThreshold, "max image threshold value (default 255)");
197  app.add_flag("--sort", sortCnt,"to sort the resulting freemanchain by decreasing size." );
198  app.add_option("-s,--minSize", minSize,"minSize of the extracted freeman chain (default 0)" );
199  app.add_option("contourSelect",cntConstraints,"Select contour according reference point and maximal distance: ex. --contourSelect X Y distanceMax" )
200  -> expected(3);
201  app.add_option("-r,--thresholdRangeMin",vectRangeMin, "use a range interval as threshold (from min) : --thresholdRangeMin min increment max : for each possible i, it define a digital sets [min, min+((i+1)*increment)] such that min+((i+1)*increment)< max and extract their boundary." )
202  -> expected(3);
203  app.add_option("-R,--thresholdRangeMax",vectRangeMax, "use a range interval as threshold (from max) : --thresholdRangeMax min increment max : for each possible i, it define a digital sets [ max-((i)*increment), max] such that max-((i)*increment)>min and extract their boundary." )
204  -> expected(3);
205 
206 
207  app.get_formatter()->column_width(40);
208  CLI11_PARSE(app, argc, argv);
209  // END parse command line using CLI ----------------------------------------------
210 
211 
212  bool thresholdRange=vectRangeMax.size()==3 || vectRangeMin.size()==3;
213  typedef functors::IntervalThresholder<Image::Value> Binarizer;
214  Image image = GenericReader<Image>::import( inputFileName );
215 
216 
217  if(cntConstraints.size()==3){
218  select=true;
219  selectCenter[0]= cntConstraints.at(0);
220  selectCenter[1]= cntConstraints.at(1);
221  selectDistanceMax= (unsigned int) cntConstraints.at(2);
222  }
223 
224  int min, max, increment;
225  if(! thresholdRange){
226  min=(int)minThreshold;
227  max= (int)maxThreshold;
228  increment = (int)(maxThreshold - minThreshold);
229  if(minThreshold == 128 && maxThreshold == 255) {
230  min=0;
231  trace.info() << "Min/Max threshold values not specified, set min to 0 and computing max with the otsu algorithm...";
232  max = getOtsuThreshold(image);
233  trace.info() << "[done] (max= " << max << ") "<< std::endl;
234  }
235 
236  }else{
237  vectRange = (vectRangeMin.size()==3) ? vectRangeMin : vectRangeMax;
238  min=vectRange.at(0);
239  increment=vectRange.at(1);
240  max = vectRange.at(2);
241  minThreshold=min;
242  maxThreshold=max;
243  }
244 
245 
246  Z2i::KSpace ks;
247  if(! ks.init( image.domain().lowerBound(),
248  image.domain().upperBound(), true )){
249  trace.error() << "Problem in KSpace initialisation"<< std::endl;
250  }
251 
252 
253  if (!thresholdRange){
254  Binarizer b(min, max);
255  functors::PointFunctorPredicate<Image,Binarizer> predicate(image, b);
256  trace.info() << "DGtal contour extraction from thresholds ["<< min << "," << max << "]" ;
257  SurfelAdjacency<2> sAdj( true );
258  std::vector< std::vector< Z2i::Point > > vectContoursBdryPointels;
259  Surfaces<Z2i::KSpace>::extractAllPointContours4C( vectContoursBdryPointels,
260  ks, predicate, sAdj );
261  if(select){
262  saveSelContoursAsFC(vectContoursBdryPointels, minSize, selectCenter, selectDistanceMax, sortCnt);
263  }else{
264  saveAllContoursAsFc(vectContoursBdryPointels, minSize, sortCnt);
265  }
266  }else{
267  for(int i=0; minThreshold+i*increment< maxThreshold; i++){
268  if(vectRangeMin.size()==3){
269  min = (int)(minThreshold+(i)*increment);
270  }
271  if(vectRangeMax.size()==3){
272  max = (int)(maxThreshold-(i)*increment);
273  }
274  Binarizer b(min, max);
275  functors::PointFunctorPredicate<Image,Binarizer> predicate(image, b);
276 
277  trace.info() << "DGtal contour extraction from thresholds ["<< min << "," << max << "]" ;
278  SurfelAdjacency<2> sAdj( true );
279  std::vector< std::vector< Z2i::Point > > vectContoursBdryPointels;
280  Surfaces<Z2i::KSpace>::extractAllPointContours4C( vectContoursBdryPointels,
281  ks, predicate, sAdj );
282  if(select){
283  saveSelContoursAsFC(vectContoursBdryPointels, minSize,
284  selectCenter, selectDistanceMax, sortCnt);
285  }else{
286  saveAllContoursAsFc(vectContoursBdryPointels,
287  minSize, sortCnt);
288  }
289  trace.info() << " [done]" << std::endl;
290  }
291  }
292 }
293 
Definition: ATu0v1.h:57