MITK-IGT
IGT Extension of MITK
Loading...
Searching...
No Matches
mitkHummelProtocolEvaluation.cpp
Go to the documentation of this file.
1/*============================================================================
2
3The Medical Imaging Interaction Toolkit (MITK)
4
5Copyright (c) German Cancer Research Center (DKFZ)
6All rights reserved.
7
8Use of this source code is governed by a 3-clause BSD license that can be
9found in the LICENSE file.
10
11============================================================================*/
12#define _USE_MATH_DEFINES
13#include "mitkHummelProtocolEvaluation.h"
14
15#include <boost/serialization/array_wrapper.hpp>
16#include <boost/accumulators/accumulators.hpp>
17#include <boost/accumulators/statistics.hpp>
18#include <boost/accumulators/statistics/stats.hpp>
19#include <boost/accumulators/statistics/mean.hpp>
20#include <boost/accumulators/statistics/moment.hpp>
21#include <boost/accumulators/statistics/tail_quantile.hpp>
22#include <boost/math/distributions/normal.hpp>
23#include <algorithm>
24
25bool mitk::HummelProtocolEvaluation::Evaluate15cmDistances(mitk::PointSet::Pointer p, HummelProtocolMeasurementVolume m, std::vector<HummelProtocolDistanceError> &Results)
26{
27 if (m != mitk::HummelProtocolEvaluation::standard) { MITK_WARN << "15 cm distances are only evaluated for standard volumes, aborting!"; return false; }
28
29 MITK_INFO << "########### 15 cm distance errors #############";
30
31 //convert measurements to matrix
32 std::array<std::array<mitk::Point3D, 10> ,9> matrix = ParseMatrixStandardVolume(p);
33
34 //these are the variables for the results:
35 std::vector<double> distances;
36 std::vector<std::string> descriptions;
37
38 //evaluation of rows
39 int distanceCounter = 0;
40 for (int row = 0; row < 9; row++) //rows
41 for (int distance = 0; distance < 7; distance++)
42 {
43 distanceCounter++;
44 mitk::Point3D point1 = p->GetPoint(row * 10 + distance);
45 mitk::Point3D point2 = p->GetPoint(row * 10 + distance + 3);
46 distances.push_back(point1.EuclideanDistanceTo(point2));
47 std::stringstream description;
48 description << "Distance(" << distanceCounter << ") " << (row + 1) << "/" << (distance + 1) << " to " << (row + 1) << "/" << (distance + 4);
49 descriptions.push_back(description.str());
50 }
51
52 //evaluation of columns
53 for (int column = 0; column < 10; column++)
54 for (int row = 0; row < 6; row++)
55 {
56 distanceCounter++;
57 mitk::Point3D point1 = matrix[row][column];
58 mitk::Point3D point2 = matrix[row + 3][column];
59 distances.push_back(point1.EuclideanDistanceTo(point2));
60 std::stringstream description;
61 description << "Distance(" << distanceCounter << ") " << (row + 1) << "/" << (column + 1) << " to " << (row + 4) << "/" << (column + 1);
62 descriptions.push_back(description.str());
63 }
64
65 //compute all errors
66 for (std::size_t i = 0; i < distances.size(); ++i)
67 {
68 HummelProtocolDistanceError currentError;
69 currentError.distanceError = fabs(distances.at(i) - (double)150.0);
70 currentError.description = descriptions.at(i);
71 Results.push_back(currentError);
72 MITK_INFO << "Error " << currentError.description << " : " << currentError.distanceError;
73 }
74
75 //compute statistics
76 std::vector<HummelProtocolDistanceError> statistics = mitk::HummelProtocolEvaluation::ComputeStatistics(Results);
77 for (auto currentError : statistics)
78 {
79 Results.push_back(currentError);
80 MITK_INFO << currentError.description << " : " << currentError.distanceError;
81 }
82
83 return true;
84}
85
86bool mitk::HummelProtocolEvaluation::Evaluate30cmDistances(mitk::PointSet::Pointer p, HummelProtocolMeasurementVolume m, std::vector<HummelProtocolDistanceError> &Results)
87{
88 if (m != mitk::HummelProtocolEvaluation::standard) { MITK_WARN << "30 cm distances are only evaluated for standard volumes, aborting!"; return false; }
89 MITK_INFO << "########### 30 cm distance errors #############";
90
91 //convert measurements to matrix
92 std::array<std::array<mitk::Point3D, 10> ,9> matrix = ParseMatrixStandardVolume(p);
93
94 //these are the variables for the results:
95 std::vector<double> distances;
96 std::vector<std::string> descriptions;
97
98 //evaluation of rows
99 int distanceCounter = 0;
100 for (int row = 0; row < 9; row++) //rows
101 for (int distance = 0; distance < 4; distance++)
102 {
103 distanceCounter++;
104 mitk::Point3D point1 = p->GetPoint(row * 10 + distance);
105 mitk::Point3D point2 = p->GetPoint(row * 10 + distance + 6);
106 distances.push_back(point1.EuclideanDistanceTo(point2));
107 std::stringstream description;
108 description << "Distance(" << distanceCounter << ") " << (row + 1) << "/" << (distance + 1) << " to " << (row + 1) << "/" << (distance + 7);
109 descriptions.push_back(description.str());
110 }
111
112 //evaluation of columns
113 for (int column = 0; column < 10; column++)
114 for (int row = 0; row < 3; row++)
115 {
116 distanceCounter++;
117 mitk::Point3D point1 = matrix[row][column];
118 mitk::Point3D point2 = matrix[row + 6][column];
119 distances.push_back(point1.EuclideanDistanceTo(point2));
120 std::stringstream description;
121 description << "Distance(" << distanceCounter << ") " << (row + 1) << "/" << (column + 1) << " to " << (row + 7) << "/" << (column + 1);
122 descriptions.push_back(description.str());
123 }
124
125 //compute all errors
126 for (std::size_t i = 0; i < distances.size(); ++i)
127 {
128 HummelProtocolDistanceError currentError;
129 currentError.distanceError = fabs(distances.at(i) - (double)300.0);
130 currentError.description = descriptions.at(i);
131 Results.push_back(currentError);
132 MITK_INFO << "Error " << currentError.description << " : " << currentError.distanceError;
133 }
134
135 //compute statistics
136 std::vector<HummelProtocolDistanceError> statistics = mitk::HummelProtocolEvaluation::ComputeStatistics(Results);
137 for (auto currentError : statistics)
138 {
139 Results.push_back(currentError);
140 MITK_INFO << currentError.description << " : " << currentError.distanceError;
141 }
142
143 return true;
144}
145
146bool mitk::HummelProtocolEvaluation::EvaluateAccumulatedDistances(mitk::PointSet::Pointer p, HummelProtocolMeasurementVolume m, std::vector<HummelProtocolDistanceError> &Results)
147{
148 if (m != mitk::HummelProtocolEvaluation::standard) { MITK_WARN << "Accumulated distances are only evaluated for standard volumes, aborting!"; return false; }
149
150 MITK_INFO << "########### accumulated distance errors #############";
151
152 int distanceCounter = 0;
153
154 //evaluation of rows
155 for (int row = 0; row < 9; row++) //rows
156 for (int distance = 0; distance < 9; distance++)
157 {
158 distanceCounter++;
159 mitk::Point3D point1 = p->GetPoint(row * 10);
160 mitk::Point3D point2 = p->GetPoint(row * 10 + distance + 1);
161 std::stringstream description;
162 description << "Distance(" << distanceCounter << ") " << (row + 1) << "/1 to " << (row + 1) << "/" << (distance + 2);
163 //compute error
164 HummelProtocolDistanceError currentError;
165 currentError.distanceError = fabs(point1.EuclideanDistanceTo(point2) - (double)(50.0*(distance+1)));
166 currentError.description = description.str();
167 Results.push_back(currentError);
168 MITK_INFO << "Error " << currentError.description << " : " << currentError.distanceError;
169 }
170
171
172
173 //compute statistics
174 std::vector<HummelProtocolDistanceError> statistics = mitk::HummelProtocolEvaluation::ComputeStatistics(Results);
175 for (auto currentError : statistics)
176 {
177 Results.push_back(currentError);
178 MITK_INFO << currentError.description << " : " << currentError.distanceError;
179 }
180
181return true;
182}
183
184
185
186bool mitk::HummelProtocolEvaluation::Evaluate5cmDistances(mitk::PointSet::Pointer p, HummelProtocolMeasurementVolume m, std::vector<HummelProtocolDistanceError> &Results)
187{
188MITK_INFO << "########### 5 cm distance errors #############";
189std::vector<double> distances;
190std::vector<std::string> descriptions;
191switch (m)
192{
193case small:
194 if (p->GetSize() != 12) {
195 MITK_WARN << "Wrong number of points: " << p->GetSize() << " (expected 12), aborting";
196 return false;
197 }
198 MITK_INFO << "Computing Hummel protocol distance errors for small measurement volumes (12 points)...";
199
200 //row 1
201 distances.push_back(p->GetPoint(0).EuclideanDistanceTo(p->GetPoint(1))); //0
202 descriptions.push_back("Distance 4/4 to 4/5");
203 distances.push_back(p->GetPoint(1).EuclideanDistanceTo(p->GetPoint(2))); //1
204 descriptions.push_back("Distance 4/5 to 4/6");
205 distances.push_back(p->GetPoint(2).EuclideanDistanceTo(p->GetPoint(3))); //2
206 descriptions.push_back("Distance 4/6 to 4/7");
207 //row 2
208 distances.push_back(p->GetPoint(4).EuclideanDistanceTo(p->GetPoint(5))); //3
209 descriptions.push_back("Distance 5/4 to 5/5");
210 distances.push_back(p->GetPoint(5).EuclideanDistanceTo(p->GetPoint(6))); //4
211 descriptions.push_back("Distance 5/5 to 5/6");
212 distances.push_back(p->GetPoint(6).EuclideanDistanceTo(p->GetPoint(7))); //5
213 descriptions.push_back("Distance 5/6 to 5/7");
214 //row 3
215 distances.push_back(p->GetPoint(8).EuclideanDistanceTo(p->GetPoint(9))); //6
216 descriptions.push_back("Distance 6/4 to 6/5");
217 distances.push_back(p->GetPoint(9).EuclideanDistanceTo(p->GetPoint(10))); //7
218 descriptions.push_back("Distance 6/5 to 6/6");
219 distances.push_back(p->GetPoint(10).EuclideanDistanceTo(p->GetPoint(11))); //8
220 descriptions.push_back("Distance 6/6 to 6/7");
221 //column 1
222 distances.push_back(p->GetPoint(0).EuclideanDistanceTo(p->GetPoint(4))); //9
223 descriptions.push_back("Distance 4/4 to 5/4");
224 distances.push_back(p->GetPoint(4).EuclideanDistanceTo(p->GetPoint(8))); //10
225 descriptions.push_back("Distance 5/4 to 6/4");
226 //column 2
227 distances.push_back(p->GetPoint(1).EuclideanDistanceTo(p->GetPoint(5))); //11
228 descriptions.push_back("Distance 4/5 to 5/5");
229 distances.push_back(p->GetPoint(5).EuclideanDistanceTo(p->GetPoint(9))); //12
230 descriptions.push_back("Distance 5/5 to 6/5");
231 //column 3
232 distances.push_back(p->GetPoint(2).EuclideanDistanceTo(p->GetPoint(6))); //13
233 descriptions.push_back("Distance 4/6 to 5/6");
234 distances.push_back(p->GetPoint(6).EuclideanDistanceTo(p->GetPoint(10))); //14
235 descriptions.push_back("Distance 5/6 to 6/6");
236 //column 4
237 distances.push_back(p->GetPoint(3).EuclideanDistanceTo(p->GetPoint(7))); //15
238 descriptions.push_back("Distance 4/7 to 5/7");
239 distances.push_back(p->GetPoint(7).EuclideanDistanceTo(p->GetPoint(11))); //16
240 descriptions.push_back("Distance 5/7 to 6/7");
241
242break;
243
244case medium:
245{
246 if (p->GetSize() != 25) {
247 MITK_WARN << "Wrong number of points (expected 25), aborting";
248 return false;
249 }
250 MITK_INFO << "Computing Hummel protocol distance errors for medium measurement volumes (25 points)...";
251
252 int distanceCounter = 0;
253
254 //convert measurements to matrix
255 std::array<std::array<mitk::Point3D, 5>, 5> matrix = ParseMatrixMediumVolume(p);
256
257 //evaluation of rows
258 for (int row = 0; row < 5; row++) //rows
259 for (int distance = 0; distance < 4; distance++)
260 {
261 distanceCounter++;
262 mitk::Point3D point1 = p->GetPoint(row * 5 + distance);
263 mitk::Point3D point2 = p->GetPoint(row * 5 + distance + 1);
264 distances.push_back(point1.EuclideanDistanceTo(point2));
265 std::stringstream description;
266 description << "Distance(" << distanceCounter << ") " << (row + 1) << "/" << (distance + 1) << " to " << (row + 1) << "/" << (distance + 2);
267 descriptions.push_back(description.str());
268 }
269
270 //evaluation of columns
271 for (int column = 0; column < 5; column++)
272 for (int row = 0; row < 4; row++)
273 {
274 distanceCounter++;
275 mitk::Point3D point1 = matrix[row][column];
276 mitk::Point3D point2 = matrix[row + 1][column];
277 MITK_INFO << "Point 1:" << point1 << "/Point 2:" << point2 << "/Distance:" << point1.EuclideanDistanceTo(point2);
278 distances.push_back(point1.EuclideanDistanceTo(point2));
279 std::stringstream description;
280 description << "Distance(" << distanceCounter << ") " << (row + 1) << "/" << (column + 1) << " to " << (row + 2) << "/" << (column + 1);
281 descriptions.push_back(description.str());
282 }
283 }
284 break;
285
286case standard:
287{
288 if (p->GetSize() != 90) {
289 MITK_WARN << "Wrong number of points (expected 90), aborting";
290 return false;
291 }
292 MITK_INFO << "Computing Hummel protocol distance errors for standard measurement volumes (90 points)...";
293
294 int distanceCounter = 0;
295
296 //convert measurements to matrix
297 std::array<std::array<mitk::Point3D, 10>, 9> matrix = ParseMatrixStandardVolume(p);
298
299 //evaluation of rows
300 for (int row = 0; row < 9; row++) //rows
301 for (int distance = 0; distance < 9; distance++)
302 {
303 distanceCounter++;
304 mitk::Point3D point1 = p->GetPoint(row * 10 + distance);
305 mitk::Point3D point2 = p->GetPoint(row * 10 + distance + 1);
306 distances.push_back(point1.EuclideanDistanceTo(point2));
307 std::stringstream description;
308 description << "Distance(" << distanceCounter << ") " << (row + 1) << "/" << (distance + 1) << " to " << (row + 1) << "/" << (distance + 2);
309 descriptions.push_back(description.str());
310 }
311
312 //evaluation of columns
313 for (int column = 0; column < 10; column++)
314 for (int row = 0; row < 8; row++)
315 {
316 distanceCounter++;
317 mitk::Point3D point1 = matrix[row][column];
318 mitk::Point3D point2 = matrix[row + 1][column];
319 distances.push_back(point1.EuclideanDistanceTo(point2));
320 std::stringstream description;
321 description << "Distance(" << distanceCounter << ") " << (row + 1) << "/" << (column + 1) << " to " << (row + 2) << "/" << (column + 1);
322 descriptions.push_back(description.str());
323 }
324 }
325 break;
326
327
328}
329
330//compute all errors
331for (std::size_t i = 0; i < distances.size(); ++i)
332{
333HummelProtocolDistanceError currentError;
334currentError.distanceError = fabs(distances.at(i) - (double)50.0);
335currentError.description = descriptions.at(i);
336Results.push_back(currentError);
337MITK_INFO << "Error " << currentError.description << " : " << currentError.distanceError;
338}
339
340//compute statistics
341std::vector<HummelProtocolDistanceError> statistics = mitk::HummelProtocolEvaluation::ComputeStatistics(Results);
342for (auto currentError : statistics)
343{
344 Results.push_back(currentError);
345 MITK_INFO << currentError.description << " : " << currentError.distanceError;
346}
347
348return true;
349}
350
351std::array<std::array<mitk::Point3D, 10>, 9> mitk::HummelProtocolEvaluation::ParseMatrixStandardVolume(mitk::PointSet::Pointer p)
352{
353
354 std::array<std::array<mitk::Point3D, 10> ,9> returnValue;
355
356 if (p->GetSize() != 90)
357 {
358 MITK_WARN << "PointSet does not have the right size. Expected 90 got " << p->GetSize() << " ... aborting!";
359 return returnValue;
360 }
361 for (int row = 0; row < 9; row++)
362 for (int column = 0; column < 10; column++)
363 returnValue[row][column] = p->GetPoint(row * 10 + column);
364 return returnValue;
365}
366
367std::array<std::array<mitk::Point3D, 5>, 5> mitk::HummelProtocolEvaluation::ParseMatrixMediumVolume(mitk::PointSet::Pointer p)
368{
369
370 std::array<std::array<mitk::Point3D, 5>, 5> returnValue;
371
372 if (p->GetSize() != 25)
373 {
374 MITK_WARN << "PointSet does not have the right size. Expected 25 got " << p->GetSize() << " ... aborting!";
375 return returnValue;
376 }
377 for (int row = 0; row < 5; row++)
378 for (int column = 0; column < 5; column++)
379 {
380 returnValue[row][column] = p->GetPoint(row * 5 + column);
381 //MITK_INFO << "m " << row << "/" << column << ":" << p->GetPoint(row * 5 + column);
382 }
383
384 return returnValue;
385
386}
387
388
389std::vector<mitk::HummelProtocolEvaluation::HummelProtocolDistanceError> mitk::HummelProtocolEvaluation::ComputeStatistics(std::vector<mitk::HummelProtocolEvaluation::HummelProtocolDistanceError> values)
390{
391 std::vector<mitk::HummelProtocolEvaluation::HummelProtocolDistanceError> returnValue;
392
393 //convert input values to boost / using boost accumulators for statistics
394 boost::accumulators::accumulator_set<double, boost::accumulators::features<boost::accumulators::tag::mean,
395 //boost::accumulators::tag::median,
396 boost::accumulators::tag::variance,
397 boost::accumulators::tag::max,
398 boost::accumulators::tag::min
399 //boost::accumulators::tag::tail<boost::accumulators::left>
400 > > acc;
401 for (mitk::HummelProtocolEvaluation::HummelProtocolDistanceError each : values)
402 {
403 acc(each.distanceError);
404 }
405
406 returnValue.push_back({ static_cast<double>(values.size()), "N" });
407 returnValue.push_back({ boost::accumulators::mean(acc), "Mean" });
408 //double quantile25th = boost::accumulators::quantile(acc, boost::accumulators::quantile_probability = 0.25);
409 //returnValue.push_back({ boost::accumulators::median(acc), "Median" });
410 //returnValue.push_back({ boost::accumulators::variance(acc), "Variance" });
411 returnValue.push_back({ boost::accumulators::min(acc), "Min" });
412 returnValue.push_back({ boost::accumulators::max(acc), "Max" });
413
414 //don't get the boost stuff working correctly, so computing the quantiles, median and standard deviation by myself:
415 std::vector<double> quantile;
416 for (mitk::HummelProtocolEvaluation::HummelProtocolDistanceError each : values)
417 {quantile.push_back(each.distanceError);}
418
419 auto const Q1 = quantile.size() / 4;
420 auto const Q2 = quantile.size() / 2;
421 auto const Q3 = Q1 + Q2;
422
423 std::sort(quantile.begin(),quantile.end());
424
425 returnValue.push_back({ quantile[Q1], "Quartile 1" });
426 returnValue.push_back({ quantile[Q2], "Median" });
427 returnValue.push_back({ quantile[Q3], "Quartile 3" });
428
429 double mean = boost::accumulators::mean(acc);
430 double errorSum = 0;
431 for (mitk::HummelProtocolEvaluation::HummelProtocolDistanceError each : values)
432 {
433 double error = pow((each.distanceError - mean),2);
434 errorSum += error;
435 }
436 double variance = errorSum / values.size();
437 double sampleVariance = errorSum / (values.size()-1);
438 double standardDev = sqrt(variance);
439 double sampleStandardDev = sqrt(sampleVariance);
440
441 returnValue.push_back({ variance, "Variance" });
442 returnValue.push_back({ sampleVariance, "Sample Variance" });
443 returnValue.push_back({ standardDev, "Standard Deviation" });
444 returnValue.push_back({ sampleStandardDev, "Sample Standard Deviation" });
445
446
447 return returnValue;
448
449}
mitk::HummelProtocolEvaluation::Evaluate15cmDistances
static bool Evaluate15cmDistances(mitk::PointSet::Pointer p, HummelProtocolMeasurementVolume m, std::vector< HummelProtocolDistanceError > &Results)
Definition mitkHummelProtocolEvaluation.cpp:25
mitk::HummelProtocolEvaluation::Evaluate30cmDistances
static bool Evaluate30cmDistances(mitk::PointSet::Pointer p, HummelProtocolMeasurementVolume m, std::vector< HummelProtocolDistanceError > &Results)
Definition mitkHummelProtocolEvaluation.cpp:86
mitk::HummelProtocolEvaluation::Evaluate5cmDistances
static bool Evaluate5cmDistances(mitk::PointSet::Pointer p, HummelProtocolMeasurementVolume m, std::vector< HummelProtocolDistanceError > &Results)
Definition mitkHummelProtocolEvaluation.cpp:186
mitk::HummelProtocolEvaluation::ComputeStatistics
static std::vector< HummelProtocolDistanceError > ComputeStatistics(std::vector< HummelProtocolDistanceError > values)
Definition mitkHummelProtocolEvaluation.cpp:389
mitk::HummelProtocolEvaluation::ParseMatrixStandardVolume
static std::array< std::array< mitk::Point3D, 10 >, 9 > ParseMatrixStandardVolume(mitk::PointSet::Pointer p)
Definition mitkHummelProtocolEvaluation.cpp:351
mitk::HummelProtocolEvaluation::ParseMatrixMediumVolume
static std::array< std::array< mitk::Point3D, 5 >, 5 > ParseMatrixMediumVolume(mitk::PointSet::Pointer p)
Definition mitkHummelProtocolEvaluation.cpp:367
mitk::HummelProtocolEvaluation::EvaluateAccumulatedDistances
static bool EvaluateAccumulatedDistances(mitk::PointSet::Pointer p, HummelProtocolMeasurementVolume m, std::vector< HummelProtocolDistanceError > &Results)
Definition mitkHummelProtocolEvaluation.cpp:146