Medical Imaging Interaction Toolkit  2018.4.99-1bab67a2
Medical Imaging Interaction Toolkit
mitkBandpassFilter.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 
13 #define _USE_MATH_DEFINES
14 #include <cmath>
15 
16 #include "mitkBandpassFilter.h"
17 
18 #include "../ITKFilter/ITKUltrasound/itkFFT1DComplexConjugateToRealImageFilter.h"
19 #include "../ITKFilter/ITKUltrasound/itkFFT1DRealToComplexConjugateImageFilter.h"
20 #include "itkComplexToModulusImageFilter.h"
21 #include "itkMultiplyImageFilter.h"
22 #include "mitkIOUtil.h"
23 #include "mitkITKImageImport.h"
24 #include "mitkImageCast.h"
25 
27  : m_HighPass(0),
28  m_LowPass(50),
29  m_HighPassAlpha(1),
30  m_LowPassAlpha(1),
31  m_FilterService(mitk::PhotoacousticFilterService::New())
32 {
33  MITK_INFO << "Instantiating BandpassFilter...";
34  SetNumberOfIndexedInputs(1);
35  SetNumberOfIndexedOutputs(1);
36  MITK_INFO << "Instantiating BandpassFilter...[Done]";
37 }
38 
40 {
41  MITK_INFO << "Destructed BandpassFilter.";
42 }
43 
45 {
46  auto input = GetInput();
47 
48  std::string type = input->GetPixelType().GetTypeAsString();
49  if (!(type == "scalar (float)" || type == " (float)"))
50  {
51  MITK_ERROR << "This filter can currently only handle float type images.";
52  mitkThrow() << "This filter can currently only handle float type images.";
53  }
54  if (m_HighPass > m_LowPass)
55  {
56  MITK_ERROR << "High Pass is higher than Low Pass; aborting.";
57  mitkThrow() << "High Pass is higher than Low Pass; aborting.";
58  }
59 }
60 
61 itk::Image<float, 3U>::Pointer BPFunction(mitk::Image::Pointer reference,
62  float cutoffFrequencyPixelHighPass,
63  float cutoffFrequencyPixelLowPass,
64  float alphaHighPass,
65  float alphaLowPass)
66 {
67  float *imageData = new float[reference->GetDimension(0) * reference->GetDimension(1)];
68 
69  float width = cutoffFrequencyPixelLowPass - cutoffFrequencyPixelHighPass;
70  float center = cutoffFrequencyPixelHighPass + width / 2.f;
71 
72  for (unsigned int n = 0; n < reference->GetDimension(1); ++n)
73  {
74  imageData[reference->GetDimension(0) * n] = 0;
75  }
76  for (int n = 0; n < width; ++n)
77  {
78  imageData[reference->GetDimension(0) * (int)(n + center - (width / 2.f))] = 1;
79  if (n <= (alphaHighPass * (width - 1)) / 2.f)
80  {
81  if (alphaHighPass > 0.00001f)
82  {
83  imageData[reference->GetDimension(0) * (int)(n + center - (width / 2.f))] =
84  (1 + cos(itk::Math::pi * (2 * n / (alphaHighPass * (width - 1)) - 1))) / 2;
85  }
86  else
87  {
88  imageData[reference->GetDimension(0) * (int)(n + center - (width / 2.f))] = 1;
89  }
90  }
91  else if (n >= (width - 1) * (1 - alphaLowPass / 2.f))
92  {
93  if (alphaLowPass > 0.00001f)
94  {
95  imageData[reference->GetDimension(0) * (int)(n + center - (width / 2.f))] =
96  (1 + cos(itk::Math::pi * (2 * n / (alphaLowPass * (width - 1)) + 1 - 2 / alphaLowPass))) / 2;
97  }
98  else
99  {
100  imageData[reference->GetDimension(0) * (int)(n + center - (width / 2.f))] = 1;
101  }
102  }
103  }
104 
105  for (unsigned int n = reference->GetDimension(1) / 2; n < reference->GetDimension(1); ++n)
106  {
107  imageData[reference->GetDimension(0) * n] =
108  imageData[(reference->GetDimension(1) - (n + 1)) * reference->GetDimension(0)];
109  }
110 
111  for (unsigned int line = 1; line < reference->GetDimension(0); ++line)
112  {
113  for (unsigned int sample = 0; sample < reference->GetDimension(1); ++sample)
114  {
115  imageData[reference->GetDimension(0) * sample + line] = imageData[reference->GetDimension(0) * sample];
116  }
117  }
118 
119  typedef itk::Image<float, 3U> ImageType;
120  ImageType::RegionType region;
121  ImageType::IndexType start;
122  start.Fill(0);
123 
124  region.SetIndex(start);
125 
126  ImageType::SizeType size;
127  size[0] = reference->GetDimension(0);
128  size[1] = reference->GetDimension(1);
129  size[2] = reference->GetDimension(2);
130 
131  region.SetSize(size);
132 
133  ImageType::SpacingType SpacingItk;
134  SpacingItk[0] = reference->GetGeometry()->GetSpacing()[0];
135  SpacingItk[1] = reference->GetGeometry()->GetSpacing()[1];
136  SpacingItk[2] = reference->GetGeometry()->GetSpacing()[2];
137 
138  ImageType::Pointer image = ImageType::New();
139  image->SetRegions(region);
140  image->Allocate();
141  image->FillBuffer(itk::NumericTraits<float>::Zero);
142  image->SetSpacing(SpacingItk);
143 
144  ImageType::IndexType pixelIndex;
145 
146  for (unsigned int slice = 0; slice < reference->GetDimension(2); ++slice)
147  {
148  for (unsigned int line = 0; line < reference->GetDimension(0); ++line)
149  {
150  for (unsigned int sample = 0; sample < reference->GetDimension(1); ++sample)
151  {
152  pixelIndex[0] = line;
153  pixelIndex[1] = sample;
154  pixelIndex[2] = slice;
155 
156  image->SetPixel(pixelIndex, imageData[line + sample * reference->GetDimension(0)]);
157  }
158  }
159  }
160 
161  delete[] imageData;
162 
163  return image;
164 }
165 
167 {
169  auto input = GetInput();
170  auto output = GetOutput();
171  mitk::Image::Pointer resampledInput = input;
172 
173  double powerOfTwo = std::log2(input->GetDimension(1));
174  int finalSize = 0;
175  double spacingResize = 1;
176 
177  // check if this is a power of two by checking that log2 is int
178  if (std::fmod(powerOfTwo, 1.0) >= std::numeric_limits<double>::epsilon())
179  {
180  finalSize = (int)pow(2, std::ceil(powerOfTwo));
181  double dim[2] = {(double)input->GetDimension(0), (double)finalSize };
182  resampledInput = m_FilterService->ApplyResamplingToDim(input, dim);
183  spacingResize = (double)input->GetDimension(1) / finalSize;
184  }
185 
186  // do a fourier transform, multiply with an appropriate window for the filter, and transform back
187  typedef itk::Image<float, 3> RealImageType;
188  RealImageType::Pointer image;
189  mitk::CastToItkImage(resampledInput, image);
191  typedef ForwardFFTFilterType::OutputImageType ComplexImageType;
192  ForwardFFTFilterType::Pointer forwardFFTFilter = ForwardFFTFilterType::New();
193  forwardFFTFilter->SetInput(image);
194  forwardFFTFilter->SetDirection(1);
195 
196  try
197  {
198  forwardFFTFilter->UpdateOutputInformation();
199  }
200  catch (itk::ExceptionObject &error)
201  {
202  std::cerr << "Error: " << error << std::endl;
203  MITK_ERROR << "Bandpass could not be applied";
204  mitkThrow() << "bandpass error";
205  }
206 
207  if (m_HighPass > m_LowPass)
208  mitkThrow() << "High pass frequency higher than low pass frequency, abort";
209 
210  float singleVoxel = spacingResize / (m_TimeSpacing * resampledInput->GetDimension(1)); // [Hz]
211  if(m_IsBFImage)
212  singleVoxel = spacingResize / (resampledInput->GetGeometry()->GetSpacing()[1] / 1e3 / m_SpeedOfSound * resampledInput->GetDimension(1)); // [Hz]
213  float cutoffPixelHighPass = std::min((m_HighPass / singleVoxel), (float)resampledInput->GetDimension(1) / 2.0f);
214  float cutoffPixelLowPass = std::min((m_LowPass / singleVoxel), (float)resampledInput->GetDimension(1) / 2.0f);
215 
216  MITK_INFO << "SingleVoxel: " << singleVoxel;
217  MITK_INFO << "HighPass: " << m_HighPass;
218  MITK_INFO << "LowPass: " << m_LowPass;
219 
220  MITK_INFO << "cutoffPixelHighPass: " << cutoffPixelHighPass;
221  MITK_INFO << "cutoffPixelLowPass: " << cutoffPixelLowPass;
222 
223  RealImageType::Pointer fftMultiplicator =
224  BPFunction(resampledInput, cutoffPixelHighPass, cutoffPixelLowPass, m_HighPassAlpha, m_LowPassAlpha);
225 
226  typedef itk::MultiplyImageFilter<ComplexImageType, RealImageType, ComplexImageType> MultiplyFilterType;
227  MultiplyFilterType::Pointer multiplyFilter = MultiplyFilterType::New();
228  multiplyFilter->SetInput1(forwardFFTFilter->GetOutput());
229  multiplyFilter->SetInput2(fftMultiplicator);
230 
231  /*
232  GrabItkImageMemory(fftMultiplicator, output);
233  mitk::IOUtil::Save(output, "D:/fftImage.nrrd");
234 
235  typedef itk::ComplexToModulusImageFilter< ComplexImageType, RealImageType > modulusType;
236  modulusType::Pointer modul = modulusType::New();
237 
238  modul->SetInput(multiplyFilter->GetOutput());
239  GrabItkImageMemory(modul->GetOutput(), output);
240  mitk::IOUtil::Save(output, "D:/fftout.nrrd");
241 
242  modul->SetInput(forwardFFTFilter->GetOutput());
243  GrabItkImageMemory(modul->GetOutput(), output);
244  mitk::IOUtil::Save(output, "D:/fftin.nrrd");*/
245 
247  InverseFilterType::Pointer inverseFFTFilter = InverseFilterType::New();
248  inverseFFTFilter->SetInput(multiplyFilter->GetOutput());
249  inverseFFTFilter->SetDirection(1);
250 
251  GrabItkImageMemory(inverseFFTFilter->GetOutput(), output);
252 
253  double dim[2] = { (double)input->GetDimension(0), (double)input->GetDimension(1) };
254  auto resampledOutput = m_FilterService->ApplyResamplingToDim(output, dim);
255 
256  output->Initialize(mitk::MakeScalarPixelType<float>(), 3, input->GetDimensions());
257  output->SetSpacing(resampledOutput->GetGeometry()->GetSpacing());
258  ImageReadAccessor copy(resampledOutput);
259  output->SetImportVolume(copy.GetData());
260 }
static char * line
Definition: svm.cpp:2870
itk::Image< float, 3U >::Pointer BPFunction(mitk::Image::Pointer reference, float cutoffFrequencyPixelHighPass, float cutoffFrequencyPixelLowPass, float alphaHighPass, float alphaLowPass)
#define MITK_INFO
Definition: mitkLogMacros.h:18
Perform the Fast Fourier Transform, in the forward direction, with real inputs, but only along one di...
Class holding methods to apply all Filters within the Photoacoustics Algorithms Module.
itk::Image< unsigned char, 3 > ImageType
#define MITK_ERROR
Definition: mitkLogMacros.h:20
DataCollection - Class to facilitate loading/accessing structured data.
Perform the Fast Fourier Transform, in the reverse direction, with real output, but only along one di...
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.
#define mitkThrow()
mitk::Image::Pointer image
static T min(T x, T y)
Definition: svm.cpp:53
InputImageType * GetInput(void)
mitk::PhotoacousticFilterService::Pointer m_FilterService
void MITKCORE_EXPORT CastToItkImage(const mitk::Image *mitkImage, itk::SmartPointer< ItkOutputImageType > &itkOutputImage)
Cast an mitk::Image to an itk::Image with a specific type.
void GenerateData() override
OutputType * GetOutput()
Get the output data of this image source object.
ImageReadAccessor class to get locked read access for a particular image part.
const void * GetData() const
Gives const access to the data.