DGtalTools  1.2.0
volShapeMetrics.cpp
1 
29 #include <iostream>
30 #include <DGtal/base/Common.h>
31 #include <DGtal/io/readers/GenericReader.h>
32 #include <DGtal/io/writers/GenericWriter.h>
33 #include <DGtal/helpers/StdDefs.h>
34 #include <DGtal/images/Image.h>
35 #include <DGtal/images/ImageContainerBySTLVector.h>
36 #include <DGtal/images/imagesSetsUtils/SetFromImage.h>
37 #include <DGtal/geometry/volumes/distance/DistanceTransformation.h>
38 #include <DGtal/math/Statistic.h>
39 
40 #include "CLI11.hpp"
41 
42 #include <limits>
43 
44 using namespace std;
45 using namespace DGtal;
46 using namespace Z3i;
47 
48 
111 typedef ImageContainerBySTLVector < Z3i::Domain, int > Image3D;
112 typedef ImageContainerBySTLVector < Z2i::Domain, int > Image2D;
113 
114 
115 bool
116 isDiff(const Image3D &imageA, int aMin, int aMax, const Image3D &imageB,
117  int bMin, int bMax, const Point &pt){
118  bool isRefOn = (imageA(pt)<= aMax) && (imageA(pt) >= aMin);
119  bool isCompOn = (imageB(pt)<= bMax) && (imageB(pt) >= bMin);
120  return ((!isRefOn && isCompOn) || (isRefOn && !isCompOn));
121 }
122 
123 
124 
125 void
126 getStatsFromDistanceMap(Statistic<double> & stats, const Image3D &imageA, int aMin, int aMax,
127  const Image3D & imageB, int bMin, int bMax,
128  bool statOnFalsePositiveOnly=false, Point *ptMax=0){
129 
130  // Get the digital set from ref image by computing the surface (use -1 and +1 since the interval of append function are open)
131  Z3i::DigitalSet set3dRef (imageA.domain());
132  SetFromImage<Z3i::DigitalSet>::append<Image3D>(set3dRef, imageA, aMin-1,aMax);
133  typedef functors::NotPointPredicate<Z3i::DigitalSet> NegPredicate;
134 
135 
136  // Applying the distance transform on the digital surface of the set:
137  typedef DistanceTransformation<Z3i::Space, NegPredicate, Z3i::L2Metric> DTL2;
138  const NegPredicate aPredicate( set3dRef );
139  DTL2 dtL2( imageA.domain(), aPredicate, Z3i::l2Metric );
140 
141  // Get the set of point of imageB: (use -1 and +1 since the interval of append function are open)
142  Z3i::DigitalSet set3dComp (imageB.domain());
143  SetFromImage<Z3i::DigitalSet>::append<Image3D>(set3dComp, imageB, bMin-1, bMax);
144 
145  unsigned int nbAdded=0;
146  double maxDist=0;
147  //Applying stats from the set to be compared (from imageB)
148  for(Z3i::DigitalSet::ConstIterator it= set3dComp.begin(); it!= set3dComp.end(); ++it){
149  if((!statOnFalsePositiveOnly) || (isDiff(imageA, aMin, aMax, imageB, bMin, bMax, *it))){
150  DTL2::Value distance = dtL2(*it);
151  stats.addValue(distance);
152  nbAdded++;
153  if(maxDist<distance){
154  maxDist=distance;
155  if(ptMax!=0){
156  (*ptMax)[0]=(*it)[0];
157  (*ptMax)[1]=(*it)[1];
158  (*ptMax)[2]=(*it)[2];
159  }
160  }
161  }
162  }
163 
164  if(nbAdded==0)
165  trace.error() << "No point added to statistics, will failed..." << endl;
166 }
167 
168 
169 void
170 getVoxelsStats(const Image3D &imageA, int aMin, int aMax, const Image3D &imageB,
171  int bMin, int bMax, bool exportStatVoxels, std::vector<Point> &vectPtBinA,
172  std::vector<Point> &vectPtCompBinCompA, std::vector<Point> &vectPtBnotInA,
173  std::vector<Point> &vectPtnotBInA, bool precisionRecallFMean ){
174  int numBinA = 0; // true positif with A as ref shape.
175  int numCompBinCompA = 0; // true neg with A as ref shape.
176  int numBnotInA = 0; // false pos with A as ref shape
177  int numNotBinA = 0; // false neg with A as ref shape
178  int numTotalInA = 0; // total pos in reference shape
179  int numTotalInB = 0; // total pos in compared shape
180  int numTotalinCompA = 0; //total in complement A
181  int numTotalinCompB = 0; //total in complement B
182 
183  for(Image3D::Domain::ConstIterator it = imageA.domain().begin(); it!=imageA.domain().end(); it++){
184  // voxels in A
185  if(imageA(*it) <= aMax && imageA(*it) >= aMin){
186  numTotalInA++;
187  //voxels in B
188  if(imageB(*it) <= bMax && imageB(*it) >= bMin){
189  numTotalInB++;
190  numBinA++;
191  if(exportStatVoxels) vectPtBinA.push_back(*it);
192  }else{
193  numTotalinCompB++;
194  numNotBinA++;
195  if(exportStatVoxels) vectPtnotBInA.push_back(*it);
196  }
197  // voxels outside A
198  }else{
199  numTotalinCompA++;
200  // voxels in B
201  if(imageB(*it) <= bMax && imageB(*it) >= bMin){
202  numBnotInA++;
203  numTotalInB++;
204  if(exportStatVoxels) vectPtBnotInA.push_back(*it);
205  }else{
206  numTotalinCompB++;
207  numCompBinCompA++;
208  if(exportStatVoxels) vectPtCompBinCompA.push_back(*it);
209  }
210  }
211  }
212 
213  std::cout << numBinA << " " << numCompBinCompA << " " << numBnotInA
214  << " " << numNotBinA << " " << numTotalInA << " "
215  << numTotalInB << " " << numTotalinCompA << " " << numTotalinCompB;
216  if(precisionRecallFMean){
217  double precision = (double)numBinA/(numBinA + numBnotInA);
218  double recall = (double)numBinA/(numBinA+numNotBinA);
219  double fmean = (2.0*precision*recall)/(precision+recall);
220  std::cout << " " << precision<< " " << recall << " " << fmean ;
221  }
222 }
223 
224 
225 // total ref: True Positive, True Negative, False Positive, False Negative
226 void
227 getVoxelsStats(const Image3D &imageA, int aMin, int aMax, const Image3D &imageB,
228  int bMin, int bMax, bool precisionRecallFMean ){
229  std::vector<Point> v1, v2, v3, v4;
230  return getVoxelsStats(imageA, aMin, aMax, imageB, bMin, bMax, false, v1, v2, v3, v4, precisionRecallFMean);
231 }
232 
233 
234 
235 
236 
237 void
238 exportSetofPoints(string filename, std::vector<Point> aVectPoint){
239  std::ofstream ofs;
240  ofs.open(filename.c_str(), std::ofstream::out );
241  ofs<< "# Set of 3d points with format: x y z" << std::endl;
242  for (unsigned int i =0; i< aVectPoint.size(); i++){
243  ofs << aVectPoint.at(i)[0] << " " << aVectPoint.at(i)[1] << " "<< aVectPoint.at(i)[2] << std::endl;
244  }
245  ofs.close();
246 }
247 
248 
249 
250 int main(int argc, char**argv)
251 {
252 
253 
254  // parse command line ----------------------------------------------
255  // A = ref
256  // B= comp
257  // parse command line using CLI ----------------------------------------------
258  CLI::App app;
259  std::string volAFilename;
260  std::string volBFilename;
261 
262  int aMin {0};
263  int aMax {128} ;
264  int bMin {0};
265  int bMax {128};
266  bool noDistanceComparisons {false};
267  bool distancesFromBnotInAOnly {false};
268  bool displayTFstats {false};
269  bool exportSDP {false};
270 
271  app.description("Apply shape measures for comparing two volumetric images A and B (shape defined from thresholds).\n It can compute: \n - voxel count from voxel partition (number of voxel from (B-A), (A-B) ...): usefull to determine classical statistics like false positive related stats.\n - euclidean distance between two volumetric images A and B\n Basic usage: \t volShapeMetrics --volA <volAFilename> --volB <volBFilename>\nTypical use :\n volShapeMetrics -a imageA.vol --aMin 128 --aMax 255 -b imageB.vol --bMin 128 --bMax 255 --distancesFromBnotInAOnly \n");
272 
273 
274 
275  app.add_option("-a,--volA,1", volAFilename, "Input filename of volume A (vol format, and other pgm3d can also be used)." )
276  ->required()
277  ->check(CLI::ExistingFile);
278 
279  app.add_option("-b,--volB,2", volBFilename, "Input filename of volume B (vol format, and other pgm3d can also be used).")
280  ->required()
281  ->check(CLI::ExistingFile);
282  app.add_option("--aMin",aMin, "min threshold for a voxel to be considered as belonging to the object of volume A. (default 0)",true);
283  app.add_option("--aMax",aMax, "max threshold for a voxel to be considered as belonging to the object of volume A. (default 128)",true);
284 
285  app.add_option("--bMin",bMin, "min threshold for a voxel to be considered as belonging to the object of volume B. (default 0)",true);
286  app.add_option("--bMax",bMax, "max threshold for a voxel to be considered as belonging to the object of volume B. (default 128)",true);
287 
288  app.add_flag("--noDistanceComparisons", noDistanceComparisons, "to avoid to apply distance map computation if the distance comparaison are not needed.");
289 
290  app.add_flag("--distancesFromBnotInAOnly",distancesFromBnotInAOnly, "apply distance map measures only for voxels of B which are not in A (else the measure are given from all distances of the object B).");
291 
292  app.add_flag("--displayTFstats", displayTFstats, "Change the comparison diplay by using the true/false/positive/negative notation and considering the shape A as reference. It also display precision/recall/f-mean statistics.");
293  app.add_flag("--exportSDP", exportSDP, "Export voxels belonging to each categorie (voxels of ( B in A) , (NOT in B and NOT in A), (B and NOT in A) and (Voxels of NOT in B and in A)).");
294 
295 
296  app.get_formatter()->column_width(40);
297  CLI11_PARSE(app, argc, argv);
298  // END parse command line using CLI ----------------------------------------------
299 
300 
301  Image3D imageA = GenericReader<Image3D>::import(volAFilename);
302  Image3D imageB = GenericReader<Image3D>::import(volBFilename);
303 
304  std::cout << "# Shape comparisons (generated with volShapeMetrics) given with the reference shape A: "<< volAFilename
305  << " (defined with threshold min: " << aMin << " and max: " << aMax << " )"<< endl
306  << "# and with the compared shape B: "<< volBFilename << " (defined with threshold min: " << bMin << " and max: " << bMax << " )"<< endl;
307  if(displayTFstats){
308  std::cout << "# #True_Positive #TrueNegative #FalsePositive #FalseNegative #TotalinA #TotalInB #TotalComplementOfRef #TotalComplementOfComp Precision Recall F-Mean ";
309  }else{
310  std::cout << "# #(Voxels of B in A) #(Voxels of NOT in B and NOT in A) #(Voxels of B and NOT in A) #(Voxels of NOT in B and in A) #(Voxels in A) #(Voxels in B) #(Voxels not in A) #(Voxels not in B) ";
311  }
312 
313  if(!noDistanceComparisons){
314  std::cout << " Max(MinDistance(shape B to shape A) Mean(MinDistance(shape B to shape A) Variance(MinDistance(shape B to shape A)) "
315  << " Mediane(MinDistance(shape B to shape A) Farthest point of B to A ";
316  if(distancesFromBnotInAOnly){
317  std::cout << "*** for parts of B which are not in A only ***" ;
318  }
319  }
320  std::cout << std::endl;
321 
322  if(exportSDP){
323  std::vector<Point> voxelsBinA, voxelsNotInBNotInA, voxelsBNotInA, voxelsNotInBInA;
324  if(displayTFstats){
325  getVoxelsStats(imageA, aMin, aMax, imageB, bMin, bMax, true, voxelsBinA, voxelsNotInBNotInA, voxelsBNotInA, voxelsNotInBInA, true);
326  exportSetofPoints("truePos.sdp", voxelsBinA);
327  exportSetofPoints("trueNeg.sdp", voxelsNotInBNotInA);
328  exportSetofPoints("falsePos.sdp", voxelsBNotInA);
329  exportSetofPoints("falseNeg.sdp", voxelsNotInBInA);
330  }else{
331  std::vector<Point> voxelsBinA, voxelsNotInBNotInA, voxelsBNotInA, voxelsNotInBInA;
332  getVoxelsStats(imageA, aMin, aMax, imageB, bMin, bMax, true, voxelsBinA, voxelsNotInBNotInA, voxelsBNotInA, voxelsNotInBInA, false);
333  exportSetofPoints("inBinA.sdp", voxelsBinA);
334  exportSetofPoints("notinBnotinA.sdp", voxelsNotInBNotInA);
335  exportSetofPoints("inBnotinA.sdp", voxelsBNotInA);
336  exportSetofPoints("notinBinA.sdp", voxelsNotInBInA);
337  }
338  }else{
339  getVoxelsStats(imageA, aMin, aMax, imageB, bMin, bMax, displayTFstats);
340  }
341 
342  if(!noDistanceComparisons){
343  trace.info() << "Computing Distance Map stats ...";
344  Statistic<double> statDistances(true);
345  Point ptMax;
346  getStatsFromDistanceMap(statDistances, imageA, aMin, aMax, imageB, bMin, bMax, false, &ptMax );
347  std::cout << " " << statDistances.max() << " " << statDistances.mean() << " " << statDistances.variance() << " "<< statDistances.median() << " " << ptMax[0] << " " << ptMax[1] << " " << ptMax[2] << std::endl;
348 
349  trace.info() << " [done] " << std::endl;
350  }
351 
352  return 1;
353 }
Definition: ATu0v1.h:57