Medical Imaging Interaction Toolkit  2018.4.99-12ad79a3
Medical Imaging Interaction Toolkit
mitkBandpassFilterTest.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 #define _USE_MATH_DEFINES
13 #include <cmath>
14 
15 #include <mitkTestFixture.h>
16 #include <mitkTestingMacros.h>
17 #include <mitkImage.h>
18 #include <mitkImageReadAccessor.h>
20 #include <random>
21 #include <mitkIOUtil.h>
22 
23 #include "../ITKFilter/ITKUltrasound/itkFFT1DRealToComplexConjugateImageFilter.h"
24 #include "mitkImageCast.h"
25 #include "mitkITKImageImport.h"
26 #include "itkComplexToModulusImageFilter.h"
27 
28 class mitkBandpassFilterTestSuite : public mitk::TestFixture
29 {
30  CPPUNIT_TEST_SUITE(mitkBandpassFilterTestSuite);
31  MITK_TEST(testHighPass);
32  MITK_TEST(testLowPass);
33  CPPUNIT_TEST_SUITE_END();
34 
35 private:
36 
37  mitk::PhotoacousticFilterService::Pointer m_FilterService;
38  const unsigned int NUM_ITERATIONS = 10;
39  const unsigned int DATA_XY_DIM = 512;
40  const unsigned int DATA_Z_DIM = 8;
41  const float TIME_SPACING = 0.00625; // [us]
42  float FREQUENCY_RESOLUTION = 1 / (TIME_SPACING / 1e6 * DATA_XY_DIM); // [Hz]
43  const float MAX_FREQUENCY = FREQUENCY_RESOLUTION * DATA_XY_DIM / 2.f; // [Hz]
44  const float HIGHPASS_FREQENCY = MAX_FREQUENCY * 0.8f; // [Hz]
45  const float LOWPASS_FREQENCY = MAX_FREQUENCY * 0.1f; // [Hz]
46  const float ALPHA = 0.01; // 0 = box, 1 = von Hann; changing this may make the test invalid
47  const float EPSILON_FFT = 0.0001f;
48 
49 public:
50 
51  void setUp() override
52  {
53  m_FilterService = mitk::PhotoacousticFilterService::New();
54  }
55 
56  void test(float HighPass, float LowPass, float HighPassAlpha, float LowPassAlpha, bool useLow, bool useHigh)
57  {
58  std::random_device r;
59  std::default_random_engine randGen(r());
60  std::uniform_real_distribution<float> randDistrHighPass(HighPass * 0.01f, HighPass * 0.2f);
61  std::uniform_real_distribution<float> randDistrLowPass(LowPass * 1.5f, LowPass * 2.f);
62  float* data = new float[DATA_XY_DIM*DATA_XY_DIM*DATA_Z_DIM];
63 
65  unsigned int dimension[3]{ DATA_XY_DIM, DATA_XY_DIM, DATA_Z_DIM };
66  inputImage->Initialize(mitk::MakeScalarPixelType<float>(), 3, dimension);
67  mitk::Vector3D spacing;
68  spacing[0] = 1;
69  spacing[1] = TIME_SPACING;
70  spacing[2] = 1;
71  inputImage->SetSpacing(spacing);
72 
73  for (unsigned int iteration = 0; iteration < NUM_ITERATIONS; ++iteration)
74  {
75  // fill image with zeros
76  for (unsigned int i = 0; i < DATA_XY_DIM*DATA_XY_DIM*DATA_Z_DIM; ++i)
77  {
78  data[i] = 0;
79  }
80 
81  // write specific frequencies to the image
82  if (useHigh)
83  addFrequency(randDistrHighPass(randGen), TIME_SPACING, data, dimension);
84  if (useLow)
85  addFrequency(randDistrLowPass(randGen), TIME_SPACING, data, dimension);
86 
87  inputImage->SetImportVolume(data, 0, 0, mitk::Image::ImportMemoryManagementType::CopyMemory);
88 
89  if (!useHigh)
90  HighPass = 0;
91 
92  if (!useLow)
93  LowPass = MAX_FREQUENCY;
94 
95  mitk::Image::Pointer outputImage = m_FilterService->ApplyBandpassFilter(inputImage, HighPass, LowPass, HighPassAlpha, LowPassAlpha, TIME_SPACING / 1e6, 0, false);
96 
97  // do a fourier transform, and check whether the part of the image that has been filtered is zero
98  typedef itk::Image< float, 3 > RealImageType;
99  RealImageType::Pointer image;
100 
101  mitk::CastToItkImage(outputImage, image);
102 
104  // typedef ForwardFFTFilterType::OutputImageType ComplexImageType;
105  ForwardFFTFilterType::Pointer forwardFFTFilter = ForwardFFTFilterType::New();
106  forwardFFTFilter->SetInput(image);
107  forwardFFTFilter->SetDirection(1);
108  forwardFFTFilter->UpdateOutputInformation();
109  forwardFFTFilter->Update();
110 
111  auto fftResult = forwardFFTFilter->GetOutput();
112 
113  /*
114  std::string img = "D:/img" + std::to_string(iteration) + ".nrrd";
115  mitk::IOUtil::Save(inputImage, img);
116  std::string res = "D:/res" + std::to_string(iteration) + ".nrrd";
117  mitk::IOUtil::Save(outputImage, res);*/
118 
119  // the resulting image should consist only of zeros, as we filtered the frequencies out
120  for (unsigned int z = 0; z < DATA_Z_DIM; ++z)
121  {
122  for (unsigned int y = 0; y < DATA_XY_DIM / 2; ++y)
123  {
124  if (y < (unsigned int)std::floor(0.95 * HighPass / FREQUENCY_RESOLUTION) || y >(unsigned int)std::ceil(1.05 * LowPass / FREQUENCY_RESOLUTION))
125  {
126  for (unsigned int x = 0; x < DATA_XY_DIM; ++x)
127  {
128  // unsigned int outPos = x + y * DATA_XY_DIM + z * DATA_XY_DIM * DATA_XY_DIM;
129  std::complex<float> value = fftResult->GetPixel({ x,y,z });
130  CPPUNIT_ASSERT_MESSAGE(std::string("Expected 0, got (" + std::to_string(value.real()) + " + " + std::to_string(value.imag()) + "i) at " + std::to_string(x) + "-" + std::to_string(y) + "-" + std::to_string(z)),
131  (value.real() < EPSILON_FFT) && (value.imag() < EPSILON_FFT));
132  }
133  }
134  }
135 
136  for (unsigned int y = DATA_XY_DIM / 2; y < DATA_XY_DIM; ++y)
137  {
138  if (y > DATA_XY_DIM - (unsigned int)std::floor(HighPass / FREQUENCY_RESOLUTION) || y < DATA_XY_DIM - (unsigned int)std::ceil(LowPass / FREQUENCY_RESOLUTION))
139  {
140  for (unsigned int x = 0; x < DATA_XY_DIM; ++x)
141  {
142  // unsigned int outPos = x + y * DATA_XY_DIM + z * DATA_XY_DIM * DATA_XY_DIM;
143  std::complex<float> value = fftResult->GetPixel({ x,y,z });
144  CPPUNIT_ASSERT_MESSAGE(std::string("Expected 0, got (" + std::to_string(value.real()) + " + " + std::to_string(value.imag()) + "i) at " + std::to_string(x) + "-" + std::to_string(y) + "-" + std::to_string(z)),
145  (value.real() < EPSILON_FFT) && (value.imag() < EPSILON_FFT));
146  }
147  }
148  }
149  }
150  }
151 
152  delete[] data;
153  }
154 
155  // write a fixed-frequency signal to the image
156  void addFrequency(float freq, float timeSpacing, float* data, unsigned int* dim)
157  {
158  for (unsigned int z = 0; z < dim[2]; ++z)
159  {
160  for (unsigned int y = 0; y < dim[1]; ++y)
161  {
162  for (unsigned int x = 0; x < dim[0]; ++x)
163  {
164  data[x + y * dim[0] + z * dim[0] * dim[1]] += std::sin(freq * timeSpacing * y);
165  }
166  }
167  }
168  }
169 
170  void testHighPass()
171  {
172  MITK_INFO << "Performing HighPass test";
173  test(HIGHPASS_FREQENCY, 0, ALPHA, ALPHA, false, true);
174  }
175 
176  void testLowPass()
177  {
178  MITK_INFO << "Performing LowPass test";
179  test(0, LOWPASS_FREQENCY, ALPHA, ALPHA, true, false);
180  }
181 
182  void tearDown() override
183  {
184  m_FilterService = nullptr;
185  }
186 };
187 
188 MITK_TEST_SUITE_REGISTRATION(mitkBandpassFilter)
MITK_TEST_SUITE_REGISTRATION(mitkImageToItk)
#define MITK_INFO
Definition: mitkLogMacros.h:18
Perform the Fast Fourier Transform, in the forward direction, with real inputs, but only along one di...
Follow Up Storage - Class to facilitate loading/accessing structured follow-up data.
Definition: testcase.h:28
#define MITK_TEST(TESTMETHOD)
Adds a test to the current test suite.
Test fixture for parameterized tests.
mitk::Image::Pointer image
static Pointer New()
void MITKCORE_EXPORT CastToItkImage(const mitk::Image *mitkImage, itk::SmartPointer< ItkOutputImageType > &itkOutputImage)
Cast an mitk::Image to an itk::Image with a specific type.