Medical Imaging Interaction Toolkit  2018.4.99-b585543d
Medical Imaging Interaction Toolkit
mitkCESTImageNormalizationFilter.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 <mitkCustomTagParser.h>
16 #include <mitkImage.h>
17 #include <mitkImageAccessByItk.h>
18 #include <mitkImageCast.h>
19 
20 #include <boost/algorithm/string.hpp>
21 
23 {
24 }
25 
27 {
28 }
29 
31 {
32  mitk::Image::ConstPointer inputImage = this->GetInput(0);
33  if ((inputImage->GetDimension() != 4))
34  {
35  mitkThrow() << "mitk::CESTImageNormalizationFilter:GenerateData works only with 4D images, sorry.";
36  return;
37  }
38 
39  auto resultMitkImage = this->GetOutput();
41 
42  auto originalTimeGeometry = this->GetInput()->GetTimeGeometry();
43  auto resultTimeGeometry = mitk::ProportionalTimeGeometry::New();
44 
45  unsigned int numberOfNonM0s = m_NonM0Indices.size();
46  resultTimeGeometry->Expand(numberOfNonM0s);
47 
48  for (unsigned int index = 0; index < numberOfNonM0s; ++index)
49  {
50  resultTimeGeometry->SetTimeStepGeometry(originalTimeGeometry->GetGeometryCloneForTimeStep(m_NonM0Indices.at(index)), index);
51  }
52  resultMitkImage->SetTimeGeometry(resultTimeGeometry);
53 
54  resultMitkImage->SetPropertyList(this->GetInput()->GetPropertyList()->Clone());
55  resultMitkImage->GetPropertyList()->SetStringProperty(mitk::CustomTagParser::m_OffsetsPropertyName.c_str(), m_RealOffsets.c_str());
56  // remove uids
57  resultMitkImage->GetPropertyList()->DeleteProperty("DICOM.0008.0018");
58  resultMitkImage->GetPropertyList()->DeleteProperty("DICOM.0020.000D");
59  resultMitkImage->GetPropertyList()->DeleteProperty("DICOM.0020.000E");
60 
61 }
62 
63 std::vector<double> ExtractOffsets(const mitk::Image* image)
64 {
65  std::vector<double> result;
66 
67  if (image)
68  {
69  std::string offsets = "";
70  std::vector<std::string> parts;
71  if (image->GetPropertyList()->GetStringProperty(mitk::CustomTagParser::m_OffsetsPropertyName.c_str(), offsets) && !offsets.empty())
72  {
73  boost::algorithm::trim(offsets);
74  boost::split(parts, offsets, boost::is_any_of(" "));
75 
76  for (auto part : parts)
77  {
78  std::istringstream iss(part);
79  iss.imbue(std::locale("C"));
80  double d;
81  iss >> d;
82  result.push_back(d);
83  }
84  }
85  }
86 
87  return result;
88 }
89 
90 
91 template <typename TPixel, unsigned int VImageDimension>
92 void mitk::CESTImageNormalizationFilter::NormalizeTimeSteps(const itk::Image<TPixel, VImageDimension>* image)
93 {
94  typedef itk::Image<TPixel, VImageDimension> ImageType;
95  typedef itk::Image<double, VImageDimension> OutputImageType;
96 
97  auto offsets = ExtractOffsets(this->GetInput());
98 
99  // determine normalization images
100  std::vector<unsigned int> mZeroIndices;
101  std::stringstream offsetsWithoutM0;
102  offsetsWithoutM0.imbue(std::locale("C"));
103  m_NonM0Indices.clear();
104  for (unsigned int index = 0; index < offsets.size(); ++index)
105  {
106  if ((offsets.at(index) < -299) || (offsets.at(index) > 299))
107  {
108  mZeroIndices.push_back(index);
109  }
110  else
111  {
112  offsetsWithoutM0 << offsets.at(index) << " ";
113  m_NonM0Indices.push_back(index);
114  }
115  }
116 
117  auto resultImage = OutputImageType::New();
118  typename ImageType::RegionType targetEntireRegion = image->GetLargestPossibleRegion();
119  targetEntireRegion.SetSize(3, m_NonM0Indices.size());
120  resultImage->SetRegions(targetEntireRegion);
121  resultImage->Allocate();
122  resultImage->FillBuffer(0);
123 
124  unsigned int numberOfTimesteps = image->GetLargestPossibleRegion().GetSize(3);
125 
126  typename ImageType::RegionType lowerMZeroRegion = image->GetLargestPossibleRegion();
127  lowerMZeroRegion.SetSize(3, 1);
128  typename ImageType::RegionType upperMZeroRegion = image->GetLargestPossibleRegion();
129  upperMZeroRegion.SetSize(3, 1);
130  typename ImageType::RegionType sourceRegion = image->GetLargestPossibleRegion();
131  sourceRegion.SetSize(3, 1);
132  typename OutputImageType::RegionType targetRegion = resultImage->GetLargestPossibleRegion();
133  targetRegion.SetSize(3, 1);
134  unsigned int targetTimestep = 0;
135  for (unsigned int sourceTimestep = 0; sourceTimestep < numberOfTimesteps; ++sourceTimestep)
136  {
137  unsigned int lowerMZeroIndex = mZeroIndices[0];
138  unsigned int upperMZeroIndex = mZeroIndices[0];
139  for (unsigned int loop = 0; loop < mZeroIndices.size(); ++loop)
140  {
141  if (mZeroIndices[loop] <= sourceTimestep)
142  {
143  lowerMZeroIndex = mZeroIndices[loop];
144  }
145  if (mZeroIndices[loop] > sourceTimestep)
146  {
147  upperMZeroIndex = mZeroIndices[loop];
148  break;
149  }
150  }
151  bool isMZero = (lowerMZeroIndex == sourceTimestep);
152 
153  double weight = 0.0;
154  if (lowerMZeroIndex == upperMZeroIndex)
155  {
156  weight = 1.0;
157  }
158  else
159  {
160  weight = 1.0 - double(sourceTimestep - lowerMZeroIndex) / double(upperMZeroIndex - lowerMZeroIndex);
161  }
162 
163 
164  if (isMZero)
165  {
166  //do nothing
167  }
168  else
169  {
170  lowerMZeroRegion.SetIndex(3, lowerMZeroIndex);
171  upperMZeroRegion.SetIndex(3, upperMZeroIndex);
172  sourceRegion.SetIndex(3, sourceTimestep);
173  targetRegion.SetIndex(3, targetTimestep);
174 
175  itk::ImageRegionConstIterator<ImageType> lowerMZeroIterator(image, lowerMZeroRegion);
176  itk::ImageRegionConstIterator<ImageType> upperMZeroIterator(image, upperMZeroRegion);
177  itk::ImageRegionConstIterator<ImageType> sourceIterator(image, sourceRegion);
178  itk::ImageRegionIterator<OutputImageType> targetIterator(resultImage.GetPointer(), targetRegion);
179 
180  while (!sourceIterator.IsAtEnd())
181  {
182  double normalizationFactor = weight * lowerMZeroIterator.Get() + (1.0 - weight) * upperMZeroIterator.Get();
183  if (mitk::Equal(normalizationFactor, 0))
184  {
185  targetIterator.Set(0);
186  }
187  else
188  {
189  targetIterator.Set(double(sourceIterator.Get()) / normalizationFactor);
190  }
191 
192  ++lowerMZeroIterator;
193  ++upperMZeroIterator;
194  ++sourceIterator;
195  ++targetIterator;
196  }
197  ++targetTimestep;
198  }
199  }
200 
201  // get Pointer to output image
202  mitk::Image::Pointer resultMitkImage = this->GetOutput();
203  // write into output image
204  mitk::CastToMitkImage<OutputImageType>(resultImage, resultMitkImage);
205 
206  m_RealOffsets = offsetsWithoutM0.str();
207 }
208 
210 {
211  mitk::Image::ConstPointer input = this->GetInput();
212  mitk::Image::Pointer output = this->GetOutput();
213 
214  itkDebugMacro(<< "GenerateOutputInformation()");
215 }
216 
217 bool mitk::IsNotNormalizedCESTImage(const Image* cestImage)
218 {
219  auto offsets = ExtractOffsets(cestImage);
220 
221  for (auto offset : offsets)
222  {
223  if (offset < -299 || offset > 299)
224  {
225  return true;
226  }
227  }
228  return false;
229 };
#define AccessFixedDimensionByItk(mitkImage, itkImageTypeFunction, dimension)
Access a mitk-image with known dimension by an itk-image.
itk::Image< unsigned char, 3 > ImageType
void GenerateData() override
Method generating the output of this filter. Called in the updated process of the pipeline...
~CESTImageNormalizationFilter() override
standard destructor
static Vector3D offset
const mitk::TimeGeometry * GetTimeGeometry() const
Return the TimeGeometry of the data as const pointer.
Definition: mitkBaseData.h:61
void NormalizeTimeSteps(const itk::Image< TPixel, VImageDimension > *image)
itk::ImageRegion< RegionDimension > RegionType
std::vector< unsigned int > m_NonM0Indices
non M0 indices
#define mitkThrow()
Image class for storing images.
Definition: mitkImage.h:72
mitk::Image OutputImageType
Some convenient typedefs.
std::string m_RealOffsets
Offsets without M0s.
mitk::Image::Pointer image
MITKCEST_EXPORT bool IsNotNormalizedCESTImage(const Image *cestImage)
static Pointer New()
mitk::PropertyList::Pointer GetPropertyList() const
Get the data&#39;s property list.
void GenerateOutputInformation() override
Method generating the output information of this filter (e.g. image dimension, image type...
InputImageType * GetInput(void)
MITKNEWMODULE_EXPORT bool Equal(mitk::ExampleDataStructure *leftHandSide, mitk::ExampleDataStructure *rightHandSide, mitk::ScalarType eps, bool verbose)
Returns true if the example data structures are considered equal.
std::vector< double > ExtractOffsets(const mitk::Image *image)
OutputType * GetOutput()
Get the output data of this image source object.
static const std::string m_OffsetsPropertyName
name of the property for the offsets, including normalization offsets