Medical Imaging Interaction Toolkit  2018.4.99-f51274ea
Medical Imaging Interaction Toolkit
mitkPAVolumeManipulator.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 #include "mitkPAExceptions.h"
16 #include "itkResampleImageFilter.h"
17 #include "itkGaussianInterpolateImageFunction.h"
18 
19 // Includes for image casting between ITK and MITK
20 #include <mitkImageCast.h>
21 #include <mitkITKImageImport.h>
22 
23 #include <mitkIOUtil.h>
24 
25 mitk::pa::VolumeManipulator::VolumeManipulator() {}
26 
27 mitk::pa::VolumeManipulator::~VolumeManipulator() {}
28 
29 void mitk::pa::VolumeManipulator::ThresholdImage(Volume::Pointer image, double threshold)
30 {
31  for (unsigned int z = 0; z < image->GetZDim(); z++)
32  for (unsigned int y = 0; y < image->GetYDim(); y++)
33  for (unsigned int x = 0; x < image->GetXDim(); x++)
34  if (image->GetData(x, y, z) > threshold)
35  image->SetData(1, x, y, z);
36  else
37  image->SetData(0, x, y, z);
38 }
39 
40 void mitk::pa::VolumeManipulator::MultiplyImage(Volume::Pointer image, double factor)
41 {
42  for (unsigned int z = 0; z < image->GetZDim(); z++)
43  for (unsigned int y = 0; y < image->GetYDim(); y++)
44  for (unsigned int x = 0; x < image->GetXDim(); x++)
45  image->SetData(image->GetData(x, y, z)*factor, x, y, z);
46 }
47 
49 {
50  for (unsigned int z = 0; z < image->GetZDim(); z++)
51  for (unsigned int y = 0; y < image->GetYDim(); y++)
52  for (unsigned int x = 0; x < image->GetXDim(); x++)
53  {
54  if (image->GetData(x, y, z) < 0)
55  {
56  MITK_ERROR << "Signal was smaller than 0. There is an error in the input image file.";
57  throw InvalidValueException("Signal was smaller than 0. There is an error in the input image file.");
58  }
59  image->SetData(log10(image->GetData(x, y, z)), x, y, z);
60  }
61 }
62 
63 void mitk::pa::VolumeManipulator::RescaleImage(InSilicoTissueVolume::Pointer image, double ratio)
64 {
65  MITK_INFO << "Rescaling images..";
66  double sigma = image->GetSpacing() / (ratio * 2);
67  image->SetAbsorptionVolume(RescaleImage(image->GetAbsorptionVolume(), ratio, sigma));
68  image->SetScatteringVolume(RescaleImage(image->GetScatteringVolume(), ratio, sigma));
69  image->SetAnisotropyVolume(RescaleImage(image->GetAnisotropyVolume(), ratio, sigma));
70  image->SetSegmentationVolume(RescaleImage(image->GetSegmentationVolume(), ratio, 0));
71  MITK_INFO << "Rescaling images..[Done]";
72 }
73 
74 mitk::pa::Volume::Pointer mitk::pa::VolumeManipulator::RescaleImage(Volume::Pointer image, double ratio, double sigma)
75 {
76  MITK_INFO << "Rescaling image..";
77  typedef itk::Image<double, 3> ImageType;
78  typedef itk::ResampleImageFilter<ImageType, ImageType> FilterType;
79  typedef itk::GaussianInterpolateImageFunction<ImageType, double> InterpolatorType;
80 
81  auto input = image->AsMitkImage();
82  ImageType::Pointer itkInput = ImageType::New();
83  mitk::CastToItkImage(input, itkInput);
84 
85  ImageType::SizeType outputSize;
86  outputSize[0] = input->GetDimensions()[0] * ratio;
87  outputSize[1] = input->GetDimensions()[1] * ratio;
88  outputSize[2] = input->GetDimensions()[2] * ratio;
89 
90  FilterType::Pointer resampleImageFilter = FilterType::New();
91  resampleImageFilter->SetInput(itkInput);
92  resampleImageFilter->SetSize(outputSize);
93  if (sigma > mitk::eps)
94  {
95  auto interpolator = InterpolatorType::New();
96  interpolator->SetSigma(sigma);
97  resampleImageFilter->SetInterpolator(interpolator);
98  }
99  resampleImageFilter->SetOutputSpacing(input->GetGeometry()->GetSpacing()[0] / ratio);
100 
101  MITK_INFO << "Update..";
102  resampleImageFilter->UpdateLargestPossibleRegion();
103  MITK_INFO << "Update..[Done]";
104 
105  ImageType::Pointer output = resampleImageFilter->GetOutput();
106  mitk::Image::Pointer mitkOutput = mitk::Image::New();
107 
108  GrabItkImageMemory(output, mitkOutput);
109  MITK_INFO << "Rescaling image..[Done]";
110  return Volume::New(mitkOutput);
111 }
112 
128 void mitk::pa::VolumeManipulator::GaussianBlur3D(mitk::pa::Volume::Pointer paVolume, double sigma)
129 {
130  double* volume = paVolume->GetData();
131  long width = paVolume->GetYDim();
132  long height = paVolume->GetXDim();
133  long depth = paVolume->GetZDim();
134  const long plane = width*height;
135  const long numel = plane*depth;
136  double lambda, dnu;
137  double nu, boundaryscale, postscale;
138  double *ptr;
139  long i, x, y, z;
140  int step;
141 
142  if (sigma <= 0)
143  return;
144 
145  lambda = (sigma*sigma) / (8.0);
146  dnu = (1.0 + 2.0*lambda - sqrt(1.0 + 4.0*lambda)) / (2.0*lambda);
147  nu = dnu;
148  boundaryscale = 1.0 / (1.0 - dnu);
149  postscale = pow(dnu / lambda, 12);
150 
151  /* Filter horizontally along each row */
152  for (z = 0; z < depth; z++)
153  {
154  for (y = 0; y < height; y++)
155  {
156  for (step = 0; step < 4; step++)
157  {
158  ptr = volume + width*(y + height*z);
159  ptr[0] *= boundaryscale;
160 
161  /* Filter rightwards */
162  for (x = 1; x < width; x++)
163  {
164  ptr[x] += nu*ptr[x - 1];
165  }
166 
167  ptr[x = width - 1] *= boundaryscale;
168  /* Filter leftwards */
169  for (; x > 0; x--)
170  {
171  ptr[x - 1] += nu*ptr[x];
172  }
173  }
174  }
175  }
176  /* Filter vertically along each column */
177  for (z = 0; z < depth; z++)
178  {
179  for (x = 0; x < width; x++)
180  {
181  for (step = 0; step < 4; step++)
182  {
183  ptr = volume + x + plane*z;
184  ptr[0] *= boundaryscale;
185 
186  /* Filter downwards */
187  for (i = width; i < plane; i += width)
188  {
189  ptr[i] += nu*ptr[i - width];
190  }
191 
192  ptr[i = plane - width] *= boundaryscale;
193 
194  /* Filter upwards */
195  for (; i > 0; i -= width)
196  {
197  ptr[i - width] += nu*ptr[i];
198  }
199  }
200  }
201  }
202 
203  /* Filter along z-dimension */
204  for (y = 0; y < height; y++)
205  {
206  for (x = 0; x < width; x++)
207  {
208  for (step = 0; step < 4; step++)
209  {
210  ptr = volume + x + width*y;
211  ptr[0] *= boundaryscale;
212 
213  for (i = plane; i < numel; i += plane)
214  {
215  ptr[i] += nu*ptr[i - plane];
216  }
217 
218  ptr[i = numel - plane] *= boundaryscale;
219 
220  for (; i > 0; i -= plane)
221  {
222  ptr[i - plane] += nu*ptr[i];
223  }
224  }
225  }
226  }
227 
228  for (i = 0; i < numel; i++)
229  {
230  volume[i] *= postscale;
231  }
232 }
static void ThresholdImage(Volume::Pointer image, double threshold)
ThresholdImage applies a binary threshold filter to this image.
#define MITK_INFO
Definition: mitkLogMacros.h:18
itk::Image< unsigned char, 3 > ImageType
#define MITK_ERROR
Definition: mitkLogMacros.h:20
static void Log10Image(Volume::Pointer image)
Image::Pointer GrabItkImageMemory(itk::SmartPointer< ItkOutputImageType > &itkimage, mitk::Image *mitkImage=nullptr, const BaseGeometry *geometry=nullptr, bool update=true)
Grabs the memory of an itk::Image (with a specific type) and puts it into an mitk::Image.The memory is managed by the mitk::Image after calling this function. The itk::Image remains valid until the mitk::Image decides to free the memory.
static void RescaleImage(InSilicoTissueVolume::Pointer image, double ratio)
mitk::Image::Pointer image
static Pointer New()
static void MultiplyImage(Volume::Pointer image, double factor)
Multiplies the image with a given factor.
void MITKCORE_EXPORT CastToItkImage(const mitk::Image *mitkImage, itk::SmartPointer< ItkOutputImageType > &itkOutputImage)
Cast an mitk::Image to an itk::Image with a specific type.
MITKCORE_EXPORT const ScalarType eps
static void GaussianBlur3D(Volume::Pointer paVolume, double sigma)
applies a Gaussian blur to an image
static Volume::Pointer New(double *data, unsigned int xDim, unsigned int yDim, unsigned int zDim, double spacing)
returns smartpointer reference to a new instance of this objects. The given data array will be freed ...
this exception is thrown if an invalid value is supposed to be processed