Thursday, April 02 2020 15:38:49
MITK
| Site: E130-PC03 |
| Build Name: Linux-g++ |
| Coverage File: ./Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp |
1 | /*============================================================================
2 |
3 | The Medical Imaging Interaction Toolkit (MITK)
4 |
5 | Copyright (c) German Cancer Research Center (DKFZ)
6 | All rights reserved.
7 |
8 | Use of this source code is governed by a 3-clause BSD license that can be
9 | found in the LICENSE file.
10 |
11 | ============================================================================*/
12 |
13 | #include "mitkImageStatisticsCalculator.h"
14 | #include <mitkExtendedLabelStatisticsImageFilter.h>
15 | #include <mitkExtendedStatisticsImageFilter.h>
16 | #include <mitkImage.h>
17 | #include <mitkImageAccessByItk.h>
18 | #include <mitkImageCast.h>
19 | #include <mitkImageStatisticsConstants.h>
20 | #include <mitkImageTimeSelector.h>
21 | #include <mitkImageToItk.h>
22 | #include <mitkMaskUtilities.h>
23 | #include <mitkMinMaxImageFilterWithIndex.h>
24 | #include <mitkMinMaxLabelmageFilterWithIndex.h>
25 | #include <mitkitkMaskImageFilter.h>
26 |
27 | namespace mitk
28 | {
29 46 | void ImageStatisticsCalculator::SetInputImage(const mitk::Image *image)
30 | {
31 46 | if (image != m_Image)
32 | {
33 46 | m_Image = image;
34 46 | this->Modified();
35 | }
36 46 | }
37 |
38 30 | void ImageStatisticsCalculator::SetMask(mitk::MaskGenerator *mask)
39 | {
40 30 | if (mask != m_MaskGenerator)
41 | {
42 30 | m_MaskGenerator = mask;
43 30 | this->Modified();
44 | }
45 30 | }
46 |
47 24 | void ImageStatisticsCalculator::SetSecondaryMask(mitk::MaskGenerator *mask)
48 | {
49 24 | if (mask != m_SecondaryMaskGenerator)
50 | {
51 2 | m_SecondaryMaskGenerator = mask;
52 2 | this->Modified();
53 | }
54 24 | }
55 |
56 24 | void ImageStatisticsCalculator::SetNBinsForHistogramStatistics(unsigned int nBins)
57 | {
58 24 | if (nBins != m_nBinsForHistogramStatistics)
59 | {
60 4 | m_nBinsForHistogramStatistics = nBins;
61 4 | this->Modified();
62 4 | this->m_UseBinSizeOverNBins = false;
63 | }
64 24 | if (m_UseBinSizeOverNBins)
65 | {
66 0 | this->Modified();
67 0 | this->m_UseBinSizeOverNBins = false;
68 | }
69 24 | }
70 |
71 0 | unsigned int ImageStatisticsCalculator::GetNBinsForHistogramStatistics() const
72 | {
73 0 | return m_nBinsForHistogramStatistics;
74 | }
75 |
76 0 | void ImageStatisticsCalculator::SetBinSizeForHistogramStatistics(double binSize)
77 | {
78 0 | if (binSize != m_binSizeForHistogramStatistics)
79 | {
80 0 | m_binSizeForHistogramStatistics = binSize;
81 0 | this->Modified();
82 0 | this->m_UseBinSizeOverNBins = true;
83 | }
84 0 | if (!m_UseBinSizeOverNBins)
85 | {
86 0 | this->Modified();
87 0 | this->m_UseBinSizeOverNBins = true;
88 | }
89 0 | }
90 |
91 0 | double ImageStatisticsCalculator::GetBinSizeForHistogramStatistics() const { return m_binSizeForHistogramStatistics; }
92 |
93 70 | mitk::ImageStatisticsContainer* ImageStatisticsCalculator::GetStatistics(LabelIndex label)
94 | {
95 70 | if (m_Image.IsNull())
96 | {
97 0 | mitkThrow() << "no image";
98 | }
99 |
100 70 | if (!m_Image->IsInitialized())
101 | {
102 0 | mitkThrow() << "Image not initialized!";
103 | }
104 |
105 70 | if (IsUpdateRequired(label))
106 | {
107 46 | auto timeGeometry = m_Image->GetTimeGeometry();
108 | // always compute statistics on all timesteps
109 92 | for (unsigned int timeStep = 0; timeStep < m_Image->GetTimeSteps(); timeStep++)
110 | {
111 46 | if (m_MaskGenerator.IsNotNull())
112 | {
113 30 | m_MaskGenerator->SetTimeStep(timeStep);
114 | //See T25625: otherwise, the mask is not computed again after setting a different time step
115 30 | m_MaskGenerator->Modified();
116 30 | m_InternalMask = m_MaskGenerator->GetMask();
117 30 | if (m_MaskGenerator->GetReferenceImage().IsNotNull())
118 | {
119 25 | m_InternalImageForStatistics = m_MaskGenerator->GetReferenceImage();
120 | }
121 | else
122 | {
123 5 | m_InternalImageForStatistics = m_Image;
124 | }
125 | }
126 | else
127 | {
128 16 | m_InternalImageForStatistics = m_Image;
129 | }
130 |
131 46 | if (m_SecondaryMaskGenerator.IsNotNull())
132 | {
133 2 | m_SecondaryMaskGenerator->SetTimeStep(timeStep);
134 2 | m_SecondaryMask = m_SecondaryMaskGenerator->GetMask();
135 | }
136 |
137 92 | ImageTimeSelector::Pointer imgTimeSel = ImageTimeSelector::New();
138 46 | imgTimeSel->SetInput(m_InternalImageForStatistics);
139 46 | imgTimeSel->SetTimeNr(timeStep);
140 46 | imgTimeSel->UpdateLargestPossibleRegion();
141 46 | imgTimeSel->Update();
142 46 | m_ImageTimeSlice = imgTimeSel->GetOutput();
143 |
144 | // Calculate statistics with/without mask
145 46 | if (m_MaskGenerator.IsNull() && m_SecondaryMaskGenerator.IsNull())
146 | {
147 | // 1) calculate statistics unmasked:
148 14 | AccessByItk_2(m_ImageTimeSlice, InternalCalculateStatisticsUnmasked, timeGeometry, timeStep)
149 | }
150 | else
151 | {
152 | // 2) calculate statistics masked
153 32 | AccessByItk_2(m_ImageTimeSlice, InternalCalculateStatisticsMasked, timeGeometry, timeStep)
154 | }
155 | }
156 | }
157 |
158 70 | auto it = m_StatisticContainers.find(label);
159 70 | if (it != m_StatisticContainers.end())
160 | {
161 70 | return (it->second).GetPointer();
162 | }
163 | else
164 | {
165 0 | mitkThrow() << "unknown label";
166 | return nullptr;
167 | }
168 | }
169 |
170 | template <typename TPixel, unsigned int VImageDimension>
171 28 | void ImageStatisticsCalculator::InternalCalculateStatisticsUnmasked(
172 | typename itk::Image<TPixel, VImageDimension> *image, const TimeGeometry *timeGeometry, TimeStepType timeStep)
173 | {
174 | typedef typename itk::Image<TPixel, VImageDimension> ImageType;
175 | typedef typename itk::ExtendedStatisticsImageFilter<ImageType> ImageStatisticsFilterType;
176 | typedef typename itk::MinMaxImageFilterWithIndex<ImageType> MinMaxFilterType;
177 |
178 | // reset statistics container if exists
179 56 | ImageStatisticsContainer::Pointer statisticContainerForImage;
180 28 | LabelIndex labelNoMask = 1;
181 28 | auto it = m_StatisticContainers.find(labelNoMask);
182 28 | if (it != m_StatisticContainers.end())
183 | {
184 0 | statisticContainerForImage = it->second;
185 | }
186 | else
187 | {
188 28 | statisticContainerForImage = ImageStatisticsContainer::New();
189 28 | statisticContainerForImage->SetTimeGeometry(const_cast<mitk::TimeGeometry*>(timeGeometry));
190 28 | m_StatisticContainers.emplace(labelNoMask, statisticContainerForImage);
191 | }
192 |
193 56 | auto statObj = ImageStatisticsContainer::ImageStatisticsObject();
194 |
195 56 | typename ImageStatisticsFilterType::Pointer statisticsFilter = ImageStatisticsFilterType::New();
196 28 | statisticsFilter->SetInput(image);
197 28 | statisticsFilter->SetCoordinateTolerance(0.001);
198 28 | statisticsFilter->SetDirectionTolerance(0.001);
199 |
200 | // TODO: this is single threaded. Implement our own image filter that does this multi threaded
201 | // typename itk::MinimumMaximumImageCalculator<ImageType>::Pointer imgMinMaxFilter =
202 | // itk::MinimumMaximumImageCalculator<ImageType>::New(); imgMinMaxFilter->SetImage(image);
203 | // imgMinMaxFilter->Compute();
204 56 | vnl_vector<int> minIndex, maxIndex;
205 |
206 28 | typename MinMaxFilterType::Pointer minMaxFilter = MinMaxFilterType::New();
207 28 | minMaxFilter->SetInput(image);
208 28 | minMaxFilter->UpdateLargestPossibleRegion();
209 28 | typename ImageType::PixelType minval = minMaxFilter->GetMin();
210 28 | typename ImageType::PixelType maxval = minMaxFilter->GetMax();
211 |
212 28 | typename ImageType::IndexType tmpMinIndex = minMaxFilter->GetMinIndex();
213 28 | typename ImageType::IndexType tmpMaxIndex = minMaxFilter->GetMaxIndex();
214 |
215 | // typename ImageType::IndexType tmpMinIndex = imgMinMaxFilter->GetIndexOfMinimum();
216 | // typename ImageType::IndexType tmpMaxIndex = imgMinMaxFilter->GetIndexOfMaximum();
217 |
218 28 | minIndex.set_size(tmpMaxIndex.GetIndexDimension());
219 28 | maxIndex.set_size(tmpMaxIndex.GetIndexDimension());
220 |
221 112 | for (unsigned int i = 0; i < tmpMaxIndex.GetIndexDimension(); i++)
222 | {
223 84 | minIndex[i] = tmpMinIndex[i];
224 84 | maxIndex[i] = tmpMaxIndex[i];
225 | }
226 |
227 28 | statObj.AddStatistic(mitk::ImageStatisticsConstants::MINIMUMPOSITION(), minIndex);
228 28 | statObj.AddStatistic(mitk::ImageStatisticsConstants::MAXIMUMPOSITION(), maxIndex);
229 |
230 | // convert m_binSize in m_nBins if necessary
231 | unsigned int nBinsForHistogram;
232 28 | if (m_UseBinSizeOverNBins)
233 | {
234 0 | nBinsForHistogram = std::max(static_cast<double>(std::ceil(maxval - minval)) / m_binSizeForHistogramStatistics,
235 0 | 10.); // do not allow less than 10 bins
236 | }
237 | else
238 | {
239 28 | nBinsForHistogram = m_nBinsForHistogramStatistics;
240 | }
241 |
242 28 | statisticsFilter->SetHistogramParameters(nBinsForHistogram, minval, maxval);
243 |
244 | try
245 | {
246 28 | statisticsFilter->Update();
247 | }
248 | catch (const itk::ExceptionObject &e)
249 | {
250 | mitkThrow() << "Image statistics calculation failed due to following ITK Exception: \n " << e.what();
251 | }
252 |
253 28 | auto voxelVolume = GetVoxelVolume<TPixel, VImageDimension>(image);
254 |
255 28 | auto numberOfPixels = image->GetLargestPossibleRegion().GetNumberOfPixels();
256 28 | auto volume = static_cast<double>(numberOfPixels) * voxelVolume;
257 28 | auto variance = statisticsFilter->GetSigma() * statisticsFilter->GetSigma();
258 28 | auto rms =
259 28 | std::sqrt(std::pow(statisticsFilter->GetMean(), 2.) + statisticsFilter->GetVariance()); // variance = sigma^2
260 |
261 28 | statObj.AddStatistic(mitk::ImageStatisticsConstants::NUMBEROFVOXELS(),
262 | static_cast<ImageStatisticsContainer::VoxelCountType>(numberOfPixels));
263 28 | statObj.AddStatistic(mitk::ImageStatisticsConstants::VOLUME(), volume);
264 28 | statObj.AddStatistic(mitk::ImageStatisticsConstants::MEAN(), statisticsFilter->GetMean());
265 28 | statObj.AddStatistic(mitk::ImageStatisticsConstants::MINIMUM(),
266 28 | static_cast<ImageStatisticsContainer::RealType>(statisticsFilter->GetMinimum()));
267 28 | statObj.AddStatistic(mitk::ImageStatisticsConstants::MAXIMUM(),
268 28 | static_cast<ImageStatisticsContainer::RealType>(statisticsFilter->GetMaximum()));
269 28 | statObj.AddStatistic(mitk::ImageStatisticsConstants::STANDARDDEVIATION(), statisticsFilter->GetSigma());
270 28 | statObj.AddStatistic(mitk::ImageStatisticsConstants::VARIANCE(), variance);
271 28 | statObj.AddStatistic(mitk::ImageStatisticsConstants::SKEWNESS(), statisticsFilter->GetSkewness());
272 28 | statObj.AddStatistic(mitk::ImageStatisticsConstants::KURTOSIS(), statisticsFilter->GetKurtosis());
273 28 | statObj.AddStatistic(mitk::ImageStatisticsConstants::RMS(), rms);
274 28 | statObj.AddStatistic(mitk::ImageStatisticsConstants::MPP(), statisticsFilter->GetMPP());
275 28 | statObj.AddStatistic(mitk::ImageStatisticsConstants::ENTROPY(), statisticsFilter->GetEntropy());
276 28 | statObj.AddStatistic(mitk::ImageStatisticsConstants::MEDIAN(), statisticsFilter->GetMedian());
277 28 | statObj.AddStatistic(mitk::ImageStatisticsConstants::UNIFORMITY(), statisticsFilter->GetUniformity());
278 28 | statObj.AddStatistic(mitk::ImageStatisticsConstants::UPP(), statisticsFilter->GetUPP());
279 28 | statObj.m_Histogram = statisticsFilter->GetHistogram().GetPointer();
280 28 | statisticContainerForImage->SetStatisticsForTimeStep(timeStep, statObj);
281 28 | }
282 |
283 | template <typename TPixel, unsigned int VImageDimension>
284 152 | double ImageStatisticsCalculator::GetVoxelVolume(typename itk::Image<TPixel, VImageDimension> *image) const
285 | {
286 152 | auto spacing = image->GetSpacing();
287 152 | double voxelVolume = 1.;
288 580 | for (unsigned int i = 0; i < image->GetImageDimension(); i++)
289 | {
290 428 | voxelVolume *= spacing[i];
291 | }
292 152 | return voxelVolume;
293 | }
294 |
295 | template <typename TPixel, unsigned int VImageDimension>
296 64 | void ImageStatisticsCalculator::InternalCalculateStatisticsMasked(typename itk::Image<TPixel, VImageDimension> *image,
297 | const TimeGeometry *timeGeometry,
298 | unsigned int timeStep)
299 | {
300 | typedef itk::Image<TPixel, VImageDimension> ImageType;
301 | typedef itk::Image<MaskPixelType, VImageDimension> MaskType;
302 | typedef typename MaskType::PixelType LabelPixelType;
303 | typedef itk::ExtendedLabelStatisticsImageFilter<ImageType, MaskType> ImageStatisticsFilterType;
304 | typedef MaskUtilities<TPixel, VImageDimension> MaskUtilType;
305 | typedef typename itk::MinMaxLabelImageFilterWithIndex<ImageType, MaskType> MinMaxLabelFilterType;
306 | typedef typename ImageType::PixelType InputImgPixelType;
307 |
308 | // workaround: if m_SecondaryMaskGenerator ist not null but m_MaskGenerator is! (this is the case if we request a
309 | // 'ignore zuero valued pixels' mask in the gui but do not define a primary mask)
310 64 | bool swapMasks = false;
311 64 | if (m_SecondaryMask.IsNotNull() && m_InternalMask.IsNull())
312 | {
313 4 | m_InternalMask = m_SecondaryMask;
314 4 | m_SecondaryMask = nullptr;
315 4 | swapMasks = true;
316 | }
317 |
318 | // maskImage has to have the same dimension as image
319 128 | typename MaskType::Pointer maskImage = MaskType::New();
320 | try
321 | {
322 | // try to access the pixel values directly (no copying or casting). Only works if mask pixels are of pixelType
323 | // unsigned short
324 64 | maskImage = ImageToItkImage<MaskPixelType, VImageDimension>(m_InternalMask);
325 | }
326 | catch (const itk::ExceptionObject &)
327 |
328 | {
329 | // if the pixel type of the mask is not short, then we have to make a copy of m_InternalMask (and cast the values)
330 | CastToItkImage(m_InternalMask, maskImage);
331 | }
332 |
333 | // if we have a secondary mask (say a ignoreZeroPixelMask) we need to combine the masks (corresponds to AND)
334 64 | if (m_SecondaryMask.IsNotNull())
335 | {
336 | // dirty workaround for a bug when pf mask + any other mask is used in conjunction. We need a proper fix for this
337 | // (Fabian Isensee is responsible and probably working on it!)
338 0 | if (m_InternalMask->GetDimension() == 2 &&
339 0 | (m_SecondaryMask->GetDimension() == 3 || m_SecondaryMask->GetDimension() == 4))
340 | {
341 0 | mitk::Image::ConstPointer old_img = m_SecondaryMaskGenerator->GetReferenceImage();
342 0 | m_SecondaryMaskGenerator->SetInputImage(m_MaskGenerator->GetReferenceImage());
343 0 | m_SecondaryMask = m_SecondaryMaskGenerator->GetMask();
344 0 | m_SecondaryMaskGenerator->SetInputImage(old_img);
345 | }
346 0 | typename MaskType::Pointer secondaryMaskImage = MaskType::New();
347 0 | secondaryMaskImage = ImageToItkImage<MaskPixelType, VImageDimension>(m_SecondaryMask);
348 |
349 | // secondary mask should be a ignore zero value pixel mask derived from image. it has to be cropped to the mask
350 | // region (which may be planar or simply smaller)
351 0 | typename MaskUtilities<MaskPixelType, VImageDimension>::Pointer secondaryMaskMaskUtil =
352 | MaskUtilities<MaskPixelType, VImageDimension>::New();
353 0 | secondaryMaskMaskUtil->SetImage(secondaryMaskImage.GetPointer());
354 0 | secondaryMaskMaskUtil->SetMask(maskImage.GetPointer());
355 0 | typename MaskType::Pointer adaptedSecondaryMaskImage = secondaryMaskMaskUtil->ExtractMaskImageRegion();
356 |
357 0 | typename itk::MaskImageFilter2<MaskType, MaskType, MaskType>::Pointer maskFilter =
358 | itk::MaskImageFilter2<MaskType, MaskType, MaskType>::New();
359 0 | maskFilter->SetInput1(maskImage);
360 0 | maskFilter->SetInput2(adaptedSecondaryMaskImage);
361 0 | maskFilter->SetMaskingValue(
362 0 | 1); // all pixels of maskImage where secondaryMaskImage==1 will be kept, all the others are set to 0
363 0 | maskFilter->UpdateLargestPossibleRegion();
364 0 | maskImage = maskFilter->GetOutput();
365 | }
366 |
367 128 | typename MaskUtilType::Pointer maskUtil = MaskUtilType::New();
368 64 | maskUtil->SetImage(image);
369 64 | maskUtil->SetMask(maskImage.GetPointer());
370 |
371 | // if mask is smaller than image, extract the image region where the mask is
372 128 | typename ImageType::Pointer adaptedImage = ImageType::New();
373 |
374 64 | adaptedImage = maskUtil->ExtractMaskImageRegion(); // this also checks mask sanity
375 |
376 | // find min, max, minindex and maxindex
377 128 | typename MinMaxLabelFilterType::Pointer minMaxFilter = MinMaxLabelFilterType::New();
378 64 | minMaxFilter->SetInput(adaptedImage);
379 64 | minMaxFilter->SetLabelInput(maskImage);
380 64 | minMaxFilter->UpdateLargestPossibleRegion();
381 |
382 | // set histogram parameters for each label individually (min/max may be different for each label)
383 | typedef typename std::map<LabelPixelType, InputImgPixelType> MapType;
384 | typedef typename std::pair<LabelPixelType, InputImgPixelType> PairType;
385 |
386 128 | std::vector<LabelPixelType> relevantLabels = minMaxFilter->GetRelevantLabels();
387 128 | MapType minVals;
388 128 | MapType maxVals;
389 128 | std::map<LabelPixelType, unsigned int> nBins;
390 |
391 188 | for (LabelPixelType label : relevantLabels)
392 | {
393 124 | minVals.insert(PairType(label, minMaxFilter->GetMin(label)));
394 124 | maxVals.insert(PairType(label, minMaxFilter->GetMax(label)));
395 |
396 | unsigned int nBinsForHistogram;
397 124 | if (m_UseBinSizeOverNBins)
398 | {
399 0 | nBinsForHistogram =
400 0 | std::max(static_cast<double>(std::ceil(minMaxFilter->GetMax(label) - minMaxFilter->GetMin(label))) /
401 0 | m_binSizeForHistogramStatistics,
402 0 | 10.); // do not allow less than 10 bins
403 | }
404 | else
405 | {
406 124 | nBinsForHistogram = m_nBinsForHistogramStatistics;
407 | }
408 |
409 124 | nBins.insert(typename std::pair<LabelPixelType, unsigned int>(label, nBinsForHistogram));
410 | }
411 |
412 128 | typename ImageStatisticsFilterType::Pointer imageStatisticsFilter = ImageStatisticsFilterType::New();
413 64 | imageStatisticsFilter->SetDirectionTolerance(0.001);
414 64 | imageStatisticsFilter->SetCoordinateTolerance(0.001);
415 64 | imageStatisticsFilter->SetInput(adaptedImage);
416 64 | imageStatisticsFilter->SetLabelInput(maskImage);
417 64 | imageStatisticsFilter->SetHistogramParametersForLabels(nBins, minVals, maxVals);
418 64 | imageStatisticsFilter->Update();
419 |
420 128 | std::list<int> labels = imageStatisticsFilter->GetRelevantLabels();
421 64 | auto it = labels.begin();
422 |
423 188 | while (it != labels.end())
424 | {
425 248 | ImageStatisticsContainer::Pointer statisticContainerForLabelImage;
426 124 | auto labelIt = m_StatisticContainers.find(*it);
427 | // reset if statisticContainer already exist
428 124 | if (labelIt != m_StatisticContainers.end())
429 | {
430 0 | statisticContainerForLabelImage = labelIt->second;
431 | }
432 | // create new statisticContainer
433 | else
434 | {
435 124 | statisticContainerForLabelImage = ImageStatisticsContainer::New();
436 124 | statisticContainerForLabelImage->SetTimeGeometry(const_cast<mitk::TimeGeometry*>(timeGeometry));
437 | // link label (*it) to statisticContainer
438 124 | m_StatisticContainers.emplace(*it, statisticContainerForLabelImage);
439 | }
440 |
441 248 | ImageStatisticsContainer::ImageStatisticsObject statObj;
442 |
443 | // find min, max, minindex and maxindex
444 | // make sure to only look in the masked region, use a masker for this
445 |
446 248 | vnl_vector<int> minIndex, maxIndex;
447 124 | mitk::Point3D worldCoordinateMin;
448 124 | mitk::Point3D worldCoordinateMax;
449 124 | mitk::Point3D indexCoordinateMin;
450 124 | mitk::Point3D indexCoordinateMax;
451 124 | m_InternalImageForStatistics->GetGeometry()->IndexToWorld(minMaxFilter->GetMinIndex(*it), worldCoordinateMin);
452 124 | m_InternalImageForStatistics->GetGeometry()->IndexToWorld(minMaxFilter->GetMaxIndex(*it), worldCoordinateMax);
453 124 | m_Image->GetGeometry()->WorldToIndex(worldCoordinateMin, indexCoordinateMin);
454 124 | m_Image->GetGeometry()->WorldToIndex(worldCoordinateMax, indexCoordinateMax);
455 |
456 124 | minIndex.set_size(3);
457 124 | maxIndex.set_size(3);
458 |
459 | // for (unsigned int i=0; i < tmpMaxIndex.GetIndexDimension(); i++)
460 496 | for (unsigned int i = 0; i < 3; i++)
461 | {
462 372 | minIndex[i] = indexCoordinateMin[i];
463 372 | maxIndex[i] = indexCoordinateMax[i];
464 | }
465 |
466 124 | statObj.AddStatistic(mitk::ImageStatisticsConstants::MINIMUMPOSITION(), minIndex);
467 124 | statObj.AddStatistic(mitk::ImageStatisticsConstants::MAXIMUMPOSITION(), maxIndex);
468 |
469 124 | assert(std::abs(minMaxFilter->GetMax(*it) - imageStatisticsFilter->GetMaximum(*it)) < mitk::eps);
470 124 | assert(std::abs(minMaxFilter->GetMin(*it) - imageStatisticsFilter->GetMinimum(*it)) < mitk::eps);
471 |
472 124 | auto voxelVolume = GetVoxelVolume<TPixel, VImageDimension>(image);
473 124 | auto numberOfVoxels =
474 124 | static_cast<unsigned long>(imageStatisticsFilter->GetSum(*it) / (double)imageStatisticsFilter->GetMean(*it));
475 124 | auto volume = static_cast<double>(numberOfVoxels) * voxelVolume;
476 248 | auto rms = std::sqrt(std::pow(imageStatisticsFilter->GetMean(*it), 2.) +
477 124 | imageStatisticsFilter->GetVariance(*it)); // variance = sigma^2
478 124 | auto variance = imageStatisticsFilter->GetSigma(*it) * imageStatisticsFilter->GetSigma(*it);
479 |
480 124 | statObj.AddStatistic(mitk::ImageStatisticsConstants::NUMBEROFVOXELS(), numberOfVoxels);
481 124 | statObj.AddStatistic(mitk::ImageStatisticsConstants::VOLUME(), volume);
482 124 | statObj.AddStatistic(mitk::ImageStatisticsConstants::MEAN(), imageStatisticsFilter->GetMean(*it));
483 124 | statObj.AddStatistic(mitk::ImageStatisticsConstants::MINIMUM(), imageStatisticsFilter->GetMinimum(*it));
484 124 | statObj.AddStatistic(mitk::ImageStatisticsConstants::MAXIMUM(), imageStatisticsFilter->GetMaximum(*it));
485 124 | statObj.AddStatistic(mitk::ImageStatisticsConstants::STANDARDDEVIATION(), imageStatisticsFilter->GetSigma(*it));
486 124 | statObj.AddStatistic(mitk::ImageStatisticsConstants::VARIANCE(), variance);
487 124 | statObj.AddStatistic(mitk::ImageStatisticsConstants::SKEWNESS(), imageStatisticsFilter->GetSkewness(*it));
488 124 | statObj.AddStatistic(mitk::ImageStatisticsConstants::KURTOSIS(), imageStatisticsFilter->GetKurtosis(*it));
489 124 | statObj.AddStatistic(mitk::ImageStatisticsConstants::RMS(), rms);
490 124 | statObj.AddStatistic(mitk::ImageStatisticsConstants::MPP(), imageStatisticsFilter->GetMPP(*it));
491 124 | statObj.AddStatistic(mitk::ImageStatisticsConstants::ENTROPY(), imageStatisticsFilter->GetEntropy(*it));
492 124 | statObj.AddStatistic(mitk::ImageStatisticsConstants::MEDIAN(), imageStatisticsFilter->GetMedian(*it));
493 124 | statObj.AddStatistic(mitk::ImageStatisticsConstants::UNIFORMITY(), imageStatisticsFilter->GetUniformity(*it));
494 124 | statObj.AddStatistic(mitk::ImageStatisticsConstants::UPP(), imageStatisticsFilter->GetUPP(*it));
495 124 | statObj.m_Histogram = imageStatisticsFilter->GetHistogram(*it).GetPointer();
496 |
497 124 | statisticContainerForLabelImage->SetStatisticsForTimeStep(timeStep, statObj);
498 124 | ++it;
499 | }
500 |
501 | // swap maskGenerators back
502 64 | if (swapMasks)
503 | {
504 4 | m_SecondaryMask = m_InternalMask;
505 4 | m_InternalMask = nullptr;
506 | }
507 64 | }
508 |
509 70 | bool ImageStatisticsCalculator::IsUpdateRequired(LabelIndex label) const
510 | {
511 70 | unsigned long thisClassTimeStamp = this->GetMTime();
512 70 | unsigned long inputImageTimeStamp = m_Image->GetMTime();
513 |
514 70 | auto it = m_StatisticContainers.find(label);
515 70 | if (it == m_StatisticContainers.end())
516 | {
517 46 | return true;
518 | }
519 |
520 24 | unsigned long statisticsTimeStamp = it->second->GetMTime();
521 |
522 24 | if (thisClassTimeStamp > statisticsTimeStamp) // inputs have changed
523 | {
524 0 | return true;
525 | }
526 |
527 24 | if (inputImageTimeStamp > statisticsTimeStamp) // image has changed
528 | {
529 0 | return true;
530 | }
531 |
532 24 | if (m_MaskGenerator.IsNotNull())
533 | {
534 12 | unsigned long maskGeneratorTimeStamp = m_MaskGenerator->GetMTime();
535 12 | if (maskGeneratorTimeStamp > statisticsTimeStamp) // there is a mask generator and it has changed
536 | {
537 0 | return true;
538 | }
539 | }
540 |
541 24 | if (m_SecondaryMaskGenerator.IsNotNull())
542 | {
543 2 | unsigned long maskGeneratorTimeStamp = m_SecondaryMaskGenerator->GetMTime();
544 2 | if (maskGeneratorTimeStamp > statisticsTimeStamp) // there is a secondary mask generator and it has changed
545 | {
546 0 | return true;
547 | }
548 | }
549 |
550 24 | return false;
551 | }
552 | } // namespace mitk
553 |
