Medical Imaging Interaction Toolkit  2016.11.0
Medical Imaging Interaction Toolkit
ImageStatisticsMiniapp_v2.cpp
Go to the documentation of this file.
1 /*===================================================================
2 
3 The Medical Imaging Interaction Toolkit (MITK)
4 
5 Copyright (c) German Cancer Research Center,
6 Division of Medical and Biological Informatics.
7 All rights reserved.
8 
9 This software is distributed WITHOUT ANY WARRANTY; without
10 even the implied warranty of MERCHANTABILITY or FITNESS FOR
11 A PARTICULAR PURPOSE.
12 
13 See LICENSE.txt or http://www.mitk.org for details.
14 
15 ===================================================================*/
16 
17 #include "mitkCommandLineParser.h"
18 #include "mitkImage.h"
19 #include <mitkImageStatisticsCalculator.h>
20 #include <mitkExtendedLabelStatisticsImageFilter.h>
21 #include <mitkPlanarFigureMaskGenerator.h>
22 #include <mitkPlanarFigure.h>
23 #include "mitkIOUtil.h"
24 #include <iostream>
25 #include <usAny.h>
26 #include <fstream>
27 #include <itkImageRegionConstIterator.h>
28 #include "mitkImageAccessByItk.h"
29 #include <boost/accumulators/accumulators.hpp>
30 #include <boost/accumulators/statistics/stats.hpp>
31 #include <boost/accumulators/statistics/mean.hpp>
32 #include <boost/accumulators/statistics/variance.hpp>
33 #include <boost/accumulators/statistics/min.hpp>
34 #include <boost/accumulators/statistics/max.hpp>
35 #include <boost/accumulators/statistics/count.hpp>
36 #include <boost/accumulators/statistics/moment.hpp>
37 #include <mitkImageMaskGenerator.h>
38 #include <mitkHotspotMaskGenerator.h>
39 #include <mitkIgnorePixelMaskGenerator.h>
40 #include <chrono>
41 
42 
43 struct statistics_res{
44  double mean, variance, min, max, count, moment;
45 };
46 
47 long getTimeInMs()
48 {
49  std::chrono::milliseconds ms = std::chrono::duration_cast< std::chrono::milliseconds > (std::chrono::system_clock::now().time_since_epoch());
50  long time = ms.count();
51  return time;
52 }
53 
54 //std::string printstats(mitk::ImageStatisticsCalculator::Statistics statisticsStruct)
55 //{
56 // std::string res = "";
57 // res += std::string("Entropy ") + std::to_string(statisticsStruct.GetEntropy()) + std::string("\n");
58 // res += std::string("Kurtosis ") + std::to_string(statisticsStruct.GetKurtosis()) + std::string("\n");
59 // res += std::string("MPP ") + std::to_string(statisticsStruct.GetMPP()) + std::string("\n");
60 // res += std::string("Max ") + std::to_string(statisticsStruct.GetMax()) + std::string("\n");
61 // res += std::string("Mean ") + std::to_string(statisticsStruct.GetMean()) + std::string("\n");
62 // res += std::string("Median ") + std::to_string(statisticsStruct.GetMedian()) + std::string("\n");
63 // res += std::string("Min ") + std::to_string(statisticsStruct.GetMin()) + std::string("\n");
64 // res += std::string("N ") + std::to_string(statisticsStruct.GetN()) + std::string("\n");
65 // res += std::string("RMS ") + std::to_string(statisticsStruct.GetRMS()) + std::string("\n");
66 // res += std::string("Skewness ") + std::to_string(statisticsStruct.GetSkewness()) + std::string("\n");
67 // res += std::string("Std ") + std::to_string(statisticsStruct.GetSigma()) + std::string("\n");
68 // res += std::string("UPP ") + std::to_string(statisticsStruct.GetUPP()) + std::string("\n");
69 // res += std::string("Uniformity ") + std::to_string(statisticsStruct.GetUniformity()) + std::string("\n");
70 // res += std::string("Variance ") + std::to_string(statisticsStruct.GetVariance()) + std::string("\n");
71 // vnl_vector<int> minIndex = statisticsStruct.GetMinIndex();
72 // vnl_vector<int> maxIndex = statisticsStruct.GetMaxIndex();
73 
74 // res += "Min Index: ";
75 // for (auto it = minIndex.begin(); it != minIndex.end(); it++)
76 // {
77 // res += std::to_string(*it) + " ";
78 // }
79 // res += std::string("\n");
80 
81 // res += "Max Index: ";
82 // for (auto it = maxIndex.begin(); it != maxIndex.end(); it++)
83 // {
84 // res += std::to_string(*it) + " ";
85 // }
86 // res += std::string("\n");
87 // return res;
88 //}
89 
90 template <class T, class P>
91 void printMap(std::map<T, P> input)
92 {
93  for (auto it = input.begin(); it != input.end(); ++it)
94  {
95  std::cout << it->first<< ": " << it->second<< std::endl;
96  }
97  std::cout << std::endl;
98 }
99 
100 template < typename TPixel, unsigned int VImageDimension >
101 void get_statistics_boost(itk::Image<TPixel, VImageDimension>* itkImage, statistics_res& res){
102  typedef itk::Image<TPixel, VImageDimension> ImageType;
103 
104  itk::ImageRegionConstIterator<ImageType> it(itkImage, itkImage->GetLargestPossibleRegion());
105 
106  int ctr=0;
107  double sum=0;
108 
109  boost::accumulators::accumulator_set<double, boost::accumulators::stats<
110  boost::accumulators::tag::mean,
111  boost::accumulators::tag::variance,
112  boost::accumulators::tag::min,
113  boost::accumulators::tag::max,
114  boost::accumulators::tag::count,
115  boost::accumulators::tag::moment<2>> > acc;
116 
117  for (it.GoToBegin(); !it.IsAtEnd(); ++it)
118  {
119  acc(it.Get());
120 // currentPixel = it.Get();
121 // sum+=currentPixel;
122 // ctr+=1;
123  }
124 
125 // res.mean=(double)sum/(double)ctr;
126  res.mean = boost::accumulators::mean(acc);
127  res.variance = boost::accumulators::variance(acc);
128  res.min = boost::accumulators::min(acc);
129  res.max = boost::accumulators::max(acc);
130  res.count = boost::accumulators::count(acc);
131  res.moment = boost::accumulators::moment<2>(acc);
132 
133  std::cout << "sum: " << sum << " N: " << ctr << " mean: " << res.mean << std::endl;
134 }
135 
136 void compute_statistics(std::string inputImageFile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic3D.nrrd",
137  std::string outfname = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic3D_statistics_new.txt",
138  std::string outfname2 = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic3D_statistics_old.txt",
139  std::string maskImageFile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic3D_someSegmentation.nrrd",
140  std::string pFfile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic3D_rectangle.pf" )
141 {
142  ofstream outFile;
143  outFile.open(outfname);
144 
145  unsigned int timeStep = 0;
146 
147  mitk::Image::Pointer inputImage = mitk::IOUtil::LoadImage(inputImageFile);
148 
149  mitk::Image::Pointer maskImage;
150  maskImage = mitk::IOUtil::LoadImage(maskImageFile);
151 
152  std::vector<mitk::BaseData::Pointer> loadedObjects;
153  loadedObjects = mitk::IOUtil::Load(pFfile);
154  mitk::BaseData::Pointer pf = loadedObjects[0];
155  mitk::PlanarFigure::Pointer planarFigure = dynamic_cast<mitk::PlanarFigure*>(pf.GetPointer());
156 
157  mitk::ImageStatisticsCalculator::Pointer calculator = mitk::ImageStatisticsCalculator::New();
158  mitk::ImageStatisticsCalculator::StatisticsContainer::Pointer result;
159  calculator->SetInputImage(inputImage);
160  calculator->SetNBinsForHistogramStatistics(100);
161 
162  long start_time;
163 
164  outFile << "Calculating Statistics on Pic3D" << std::endl << std::endl;
165 
166  outFile << "New Image Statistics Calculator" << std::endl;
167 
168  start_time = getTimeInMs();
169  outFile << "1) unmasked: " << std::endl;
170  calculator->SetMask(nullptr);
171  result = calculator->GetStatistics(timeStep);
172  outFile << result->GetAsString() << std::endl;
173  outFile << "new image statistics calculator took: " << getTimeInMs()-start_time << " ms." << std::endl << std::endl;
174 
175  start_time = getTimeInMs();
176  outFile << "2) planarfigure: " << std::endl;
177  mitk::PlanarFigureMaskGenerator::Pointer planarFigMaskExtr = mitk::PlanarFigureMaskGenerator::New();
178  planarFigMaskExtr->SetPlanarFigure(planarFigure);
179  planarFigMaskExtr->SetInputImage(inputImage);
180  calculator->SetMask(planarFigMaskExtr.GetPointer());
181  result = calculator->GetStatistics(timeStep, 1);
182  outFile << result->GetAsString() << std::endl;
183  outFile << "new image statistics calculator took: " << getTimeInMs()-start_time << " ms." << std::endl << std::endl;
184 
185  start_time = getTimeInMs();
186  outFile << "3) ignore pixel value mask: " << std::endl;
187  mitk::IgnorePixelMaskGenerator::Pointer ignPixVal = mitk::IgnorePixelMaskGenerator::New();
188  ignPixVal->SetInputImage(inputImage);
189  ignPixVal->SetIgnoredPixelValue(0);
190  calculator->SetMask(ignPixVal.GetPointer());
191  result = calculator->GetStatistics(timeStep, 1);
192  outFile << result->GetAsString() << std::endl;
193  outFile << "new image statistics calculator took: " << getTimeInMs()-start_time << " ms." << std::endl << std::endl;
194 
195  start_time = getTimeInMs();
196  outFile << "4) image mask: " << std::endl;
197  mitk::ImageMaskGenerator::Pointer binaryImageMaskGen = mitk::ImageMaskGenerator::New();
198  binaryImageMaskGen->SetImageMask(maskImage);
199  calculator->SetMask(binaryImageMaskGen.GetPointer());
200  result = calculator->GetStatistics(timeStep, 1);
201  outFile << result->GetAsString() << std::endl;
202  outFile << "new image statistics calculator took: " << getTimeInMs()-start_time << " ms." << std::endl << std::endl;
203 
204  start_time = getTimeInMs();
205  outFile << "5) hotspot mask: " << std::endl;
206  mitk::HotspotMaskGenerator::Pointer hotSpotMaskGen = mitk::HotspotMaskGenerator::New();
207  hotSpotMaskGen->SetInputImage(inputImage);
208  hotSpotMaskGen->SetHotspotRadiusInMM(10);
209  hotSpotMaskGen->SetTimeStep(0);
210  calculator->SetMask(hotSpotMaskGen.GetPointer());
211  result = calculator->GetStatistics(timeStep, 1);
212  outFile << result->GetAsString() << std::endl;
213  outFile << "new image statistics calculator took: " << getTimeInMs()-start_time << " ms." << std::endl << std::endl;
214  //mitk::IOUtil::SaveImage(hotSpotMaskGen->GetMask(), "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic2DplusTHotspot.nrrd");
215 
216  outFile << "6) all time steps (image mask): " << std::endl;
217  if (inputImage->GetTimeSteps() > 1)
218  {
219  start_time = getTimeInMs();
220  for (unsigned int i=0; i < inputImage->GetTimeSteps(); i++)
221  {
222  outFile << "timestep: " << i << std::endl;
223  calculator->SetMask(binaryImageMaskGen.GetPointer());
224  result = calculator->GetStatistics(i, 1);
225  outFile << result->GetAsString() << std::endl;
226  outFile << "new image statistics calculator took: " << getTimeInMs()-start_time << " ms." << std::endl << std::endl;
227  }
228  }
229  else
230  {
231  outFile << "Input image has only 1 time step!" << std::endl;
232  }
233 
234  outFile.flush();
235  outFile.close();
236 
237  // -----------------------------------------------------------------------------------------------------------------------------------------------------------
238 // outFile.open(outfname2);
239 
240 // std::cout << std::endl << std::endl << "calculating statistics old imgstatcalc: " << std::endl;
241 // mitk::ImageStatisticsCalculator::Statistics statisticsStruct;
242 // mitk::ImageStatisticsCalculator::Pointer calculator_old;
243 
244 
245 
246 // start_time = getTimeInMs();
247 // outFile << "1) unmasked: " << std::endl;
248 // calculator_old = mitk::ImageStatisticsCalculator::New();
249 // calculator_old->SetHistogramBinSize(100);
250 // calculator_old->SetImage(inputImage);
251 // calculator_old->SetMaskingModeToNone();
252 // calculator_old->ComputeStatistics(timeStep);
253 // statisticsStruct = calculator_old->GetStatistics(timeStep);
254 // outFile << printstats(statisticsStruct) << std::endl;
255 // outFile << "old image statistics calculator took: " << getTimeInMs()-start_time << " ms." << std::endl << std::endl;
256 
257 // start_time = getTimeInMs();
258 // outFile << "2) planarFigure: " << std::endl;
259 // calculator_old = mitk::ImageStatisticsCalculator::New();
260 // calculator_old->SetHistogramBinSize(100);
261 // calculator_old->SetImage(inputImage);
262 // calculator_old->SetPlanarFigure(planarFigure);
263 // calculator_old->SetMaskingModeToPlanarFigure();
264 // calculator_old->ComputeStatistics(timeStep);
265 // statisticsStruct = calculator_old->GetStatistics(timeStep);
266 // outFile << printstats(statisticsStruct) << std::endl;
267 // outFile << "old image statistics calculator took: " << getTimeInMs()-start_time << " ms." << std::endl << std::endl;
268 
269 // start_time = getTimeInMs();
270 // outFile << "3) IgnorePixelValue: " << std::endl;
271 // calculator_old = mitk::ImageStatisticsCalculator::New();
272 // calculator_old->SetHistogramBinSize(100);
273 // calculator_old->SetImage(inputImage);
274 // calculator_old->SetMaskingModeToNone();
275 // calculator_old->SetDoIgnorePixelValue(true);
276 // calculator_old->SetIgnorePixelValue(0);
277 // calculator_old->ComputeStatistics(timeStep);
278 // statisticsStruct = calculator_old->GetStatistics(timeStep);
279 // outFile << printstats(statisticsStruct) << std::endl;
280 // outFile << "old image statistics calculator took: " << getTimeInMs()-start_time << " ms." << std::endl << std::endl;
281 // calculator_old->SetDoIgnorePixelValue(false);
282 
283 // start_time = getTimeInMs();
284 // outFile << "4) Image Mask: " << std::endl;
285 // calculator_old = mitk::ImageStatisticsCalculator::New();
286 // calculator_old->SetHistogramBinSize(100);
287 // calculator_old->SetImage(inputImage);
288 // calculator_old->SetImageMask(maskImage);
289 // calculator_old->SetMaskingModeToImage();
290 // calculator_old->ComputeStatistics(timeStep);
291 // statisticsStruct = calculator_old->GetStatistics(timeStep);
292 // outFile << printstats(statisticsStruct) << std::endl;
293 // outFile << "old image statistics calculator took: " << getTimeInMs()-start_time << " ms." << std::endl << std::endl;
294 
295 // start_time = getTimeInMs();
296 // outFile << "5) Hotspot Mask: " << std::endl;
297 // calculator_old = mitk::ImageStatisticsCalculator::New();
298 // calculator_old->SetHistogramBinSize(100);
299 // calculator_old->SetImage(inputImage);
300 // calculator_old->SetMaskingModeToNone();
301 // calculator_old->SetCalculateHotspot(true);
302 // calculator_old->SetHotspotRadiusInMM(10);
303 // calculator_old->ComputeStatistics(timeStep);
304 // statisticsStruct = calculator_old->GetStatistics(timeStep).GetHotspotStatistics();
305 // outFile << printstats(statisticsStruct) << std::endl;
306 // outFile << "old image statistics calculator took: " << getTimeInMs()-start_time << " ms." << std::endl << std::endl;
307 // calculator_old->SetCalculateHotspot(false);
308 
309 // outFile << "6) all time steps (image mask): " << std::endl;
310 // if (inputImage->GetTimeSteps() > 1)
311 // {
312 // start_time = getTimeInMs();
313 // calculator_old = mitk::ImageStatisticsCalculator::New();
314 // calculator_old->SetHistogramBinSize(100);
315 // calculator_old->SetImage(inputImage);
316 // calculator_old->SetImageMask(maskImage);
317 // calculator_old->SetMaskingModeToImage();
318 // for (unsigned int i=0; i < inputImage->GetTimeSteps(); i++)
319 // {
320 // calculator_old->ComputeStatistics(i);
321 // statisticsStruct = calculator_old->GetStatistics(i);
322 // outFile << printstats(statisticsStruct) << std::endl;
323 // outFile << "old image statistics calculator took: " << getTimeInMs()-start_time << " ms." << std::endl << std::endl;
324 // }
325 // }
326 // else
327 // {
328 // outFile << "Input image has only 1 time step!" << std::endl;
329 // }
330 
331 // outFile.flush();
332 // outFile.close();
333 }
334 
335 int main( int argc, char* argv[] )
336 {
337  std::string inputImageFile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic3D.nrrd";
338  std::string outfname = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic3D_statistics_new.txt";
339  std::string outfname2 = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic3D_statistics_old.txt";
340  std::string maskImageFile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic3D_someSegmentation.nrrd";
341  std::string pFfile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic3D_rectangle.pf";
342  compute_statistics(inputImageFile, outfname, outfname2, maskImageFile, pFfile);
343 
344  inputImageFile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic2DplusT.nrrd";
345  outfname = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic2DplusT_statistics_new.txt";
346  outfname2 = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic2DplusT_statistics_old.txt";
347  maskImageFile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic2DplusT_someSegmentation.nrrd";
348  pFfile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic2DplusT_ellipse.pf";
349  compute_statistics(inputImageFile, outfname, outfname2, maskImageFile, pFfile);
350 
351  inputImageFile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/leber_ct.pic";
352  outfname = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/leber_ct_statistics_new.txt";
353  outfname2 = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/leber_ct_statistics_old.txt";
354  maskImageFile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/leber_ct_segmentation.nrrd";
355  pFfile = "/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/leber_ct_PF.pf";
356  compute_statistics(inputImageFile, outfname, outfname2, maskImageFile, pFfile);
357 
358 
359  return EXIT_SUCCESS;
360 }
ImageType
itk::Image< unsigned char, 3 > ImageType
Definition: mitkDummyLsetReader.cpp:23
get_statistics_boost
void get_statistics_boost(itk::Image< TPixel, VImageDimension > *itkImage, statistics_res &res)
Definition: ImageStatisticsMiniapp_v2.cpp:101
mitk::ImageMaskGenerator::New
static Pointer New()
getTimeInMs
long getTimeInMs()
Definition: ImageStatisticsMiniapp_v2.cpp:47
main
int main(int argc, char *argv[])
Definition: ImageStatisticsMiniapp_v2.cpp:335
compute_statistics
void compute_statistics(std::string inputImageFile="/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic3D.nrrd", std::string outfname="/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic3D_statistics_new.txt", std::string outfname2="/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic3D_statistics_old.txt", std::string maskImageFile="/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic3D_someSegmentation.nrrd", std::string pFfile="/home/fabian/MITK/MITK_platform_project/bin/MITK-superbuild/MITK-Data/Pic3D_rectangle.pf")
Definition: ImageStatisticsMiniapp_v2.cpp:136
max
static T max(T x, T y)
Definition: svm.cpp:70
min
static T min(T x, T y)
Definition: svm.cpp:67
mitk::PlanarFigure
Base-class for geometric planar (2D) figures, such as lines, circles, rectangles, polygons...
Definition: mitkPlanarFigure.h:54
mitk::IOUtil::Load
static DataStorage::SetOfObjects::Pointer Load(const std::string &path, DataStorage &storage)
Load a file into the given DataStorage.
Definition: mitkIOUtil.cpp:483
printMap
void printMap(std::map< T, P > input)
Definition: ImageStatisticsMiniapp_v2.cpp:91
mitk::IOUtil::LoadImage
static mitk::Image::Pointer LoadImage(const std::string &path)
LoadImage Convenience method to load an arbitrary mitkImage.
Definition: mitkIOUtil.cpp:597