Medical Imaging Interaction Toolkit  2018.4.99-3e3f1a6e
Medical Imaging Interaction Toolkit
mitkMRNormLinearStatisticBasedFilter.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 (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 
14 
15 #include "mitkImageToItk.h"
16 #include "mitkImageAccessByItk.h"
17 
18 #include "itkImageRegionIterator.h"
19 // MITK
20 #include <mitkITKImageImport.h>
21 #include <mitkImageCast.h>
22 #include <mitkImageAccessByItk.h>
23 // ITK
24 #include <itkLabelStatisticsImageFilter.h>
25 #include <itkMinimumMaximumImageCalculator.h>
26 
28 m_CenterMode(MRNormLinearStatisticBasedFilter::MEDIAN), m_TargetValue(0), m_TargetWidth(1)
29 {
30  this->SetNumberOfIndexedInputs(2);
31  this->SetNumberOfRequiredInputs(1);
32 }
33 
35 {
36 }
37 
39 {
40  // Process object is not const-correct so the const_cast is required here
41  auto* nonconstMask = const_cast< mitk::Image * >( mask );
42  this->SetNthInput(1, nonconstMask );
43 }
44 
46 {
47  return this->GetInput(1);
48 }
49 
51 {
52  Superclass::GenerateInputRequestedRegion();
53 
54  mitk::Image* input = this->GetInput();
55 
57 }
58 
60 {
61  mitk::Image::ConstPointer input = this->GetInput();
62  mitk::Image::Pointer output = this->GetOutput();
63 
64  itkDebugMacro(<< "GenerateOutputInformation()");
65 
66  output->Initialize(input->GetPixelType(), *input->GetTimeGeometry());
67  output->SetPropertyList(input->GetPropertyList()->Clone());
68 }
69 
70 template < typename TPixel, unsigned int VImageDimension >
71 void mitk::MRNormLinearStatisticBasedFilter::InternalComputeMask(itk::Image<TPixel, VImageDimension>* itkImage)
72 {
73  // Define all necessary Types
74  typedef itk::Image<TPixel, VImageDimension> ImageType;
75  typedef itk::Image<int, VImageDimension> MaskType;
76  typedef itk::LabelStatisticsImageFilter<ImageType, MaskType> FilterType;
77  typedef itk::MinimumMaximumImageCalculator<ImageType> MinMaxComputerType;
78 
79  typename MaskType::Pointer itkMask0 = MaskType::New();
80  mitk::CastToItkImage(this->GetMask(), itkMask0);
81 
82  typename ImageType::Pointer outImage = ImageType::New();
83  mitk::CastToItkImage(this->GetOutput(0), outImage);
84 
85  typename MinMaxComputerType::Pointer minMaxComputer = MinMaxComputerType::New();
86  minMaxComputer->SetImage(itkImage);
87  minMaxComputer->Compute();
88  double min = minMaxComputer->GetMinimum();
89  double max = minMaxComputer->GetMaximum();
90 
91  if (m_IgnoreOutlier)
92  {
93  typename FilterType::Pointer labelStatisticsImageFilter = FilterType::New();
94  labelStatisticsImageFilter->SetUseHistograms(true);
95  labelStatisticsImageFilter->SetHistogramParameters(2048, min, max);
96  labelStatisticsImageFilter->SetInput(itkImage);
97 
98  labelStatisticsImageFilter->SetLabelInput(itkMask0);
99  labelStatisticsImageFilter->Update();
100  auto histo = labelStatisticsImageFilter->GetHistogram(1);
101  min = histo->Quantile(0, 0.02);
102  max = histo->Quantile(0, 0.98);
103  }
104 
105 
106  typename FilterType::Pointer labelStatisticsImageFilter = FilterType::New();
107  labelStatisticsImageFilter->SetUseHistograms(true);
108  labelStatisticsImageFilter->SetHistogramParameters(256, min, max);
109  labelStatisticsImageFilter->SetInput( itkImage );
110 
111  labelStatisticsImageFilter->SetLabelInput(itkMask0);
112  labelStatisticsImageFilter->Update();
113  double median0 = labelStatisticsImageFilter->GetMedian(1);
114  double mean0 = labelStatisticsImageFilter->GetMean(1);
115  double stddev = labelStatisticsImageFilter->GetSigma(1);
116  double modulo0=0;
117  {
118  auto histo = labelStatisticsImageFilter->GetHistogram(1);
119  double maxFrequency=0;
120  for (auto hIter=histo->Begin();hIter!=histo->End();++hIter)
121  {
122  if (maxFrequency < hIter.GetFrequency())
123  {
124  maxFrequency = hIter.GetFrequency();
125  modulo0 = (histo->GetBinMin(0,hIter.GetInstanceIdentifier()) + histo->GetBinMax(0,hIter.GetInstanceIdentifier())) / 2.0;
126  }
127  }
128  }
129 
130  double value0=0;
131  switch (m_CenterMode)
132  {
134  value0=mean0; break;
136  value0=median0; break;
138  value0=modulo0; break;
139  }
140 
141  double offset = value0+m_TargetValue;
142  double scaling = stddev*m_TargetWidth;
143  if (scaling < 0.0001)
144  return;
145 
146  itk::ImageRegionIterator<ImageType> inIter(itkImage, itkImage->GetLargestPossibleRegion());
147  itk::ImageRegionIterator<ImageType> outIter(outImage, outImage->GetLargestPossibleRegion());
148  while (! inIter.IsAtEnd())
149  {
150  TPixel value = inIter.Value();
151  if (m_IgnoreOutlier && (value < min))
152  {
153  value = min;
154  }
155  else if (m_IgnoreOutlier && (value > max))
156  {
157  value = max;
158  }
159 
160  outIter.Set((value - offset) / scaling);
161  ++inIter;
162  ++outIter;
163  }
164 }
165 
167 {
169 }
void SetRequestedRegionToLargestPossibleRegion() override
itk::Image< unsigned char, 3 > ImageType
void InternalComputeMask(itk::Image< TPixel, VImageDimension > *itkImage)
static Vector3D offset
Image class for storing images.
Definition: mitkImage.h:72
#define AccessByItk(mitkImage, itkImageTypeFunction)
Access a MITK image by an ITK image.
static T max(T x, T y)
Definition: svm.cpp:56
static T min(T x, T y)
Definition: svm.cpp:53
InputImageType * GetInput(void)
void MITKCORE_EXPORT CastToItkImage(const mitk::Image *mitkImage, itk::SmartPointer< ItkOutputImageType > &itkOutputImage)
Cast an mitk::Image to an itk::Image with a specific type.
mitk::Image::Pointer mask
OutputType * GetOutput()
Get the output data of this image source object.