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