diff --git a/Modules/PhotoacousticsAlgorithms/source/filters/mitkBandpassFilter.cpp b/Modules/PhotoacousticsAlgorithms/source/filters/mitkBandpassFilter.cpp index a7f0232883..255f9db70c 100644 --- a/Modules/PhotoacousticsAlgorithms/source/filters/mitkBandpassFilter.cpp +++ b/Modules/PhotoacousticsAlgorithms/source/filters/mitkBandpassFilter.cpp @@ -1,235 +1,256 @@ /*=================================================================== mitkCastToFloatImageFilter The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ +#define _USE_MATH_DEFINES +#include + #include "mitkBandpassFilter.h" #include "../ITKFilter/ITKUltrasound/itkFFT1DComplexConjugateToRealImageFilter.h" #include "../ITKFilter/ITKUltrasound/itkFFT1DRealToComplexConjugateImageFilter.h" #include "mitkImageCast.h" #include "itkMultiplyImageFilter.h" #include "mitkITKImageImport.h" #include "mitkIOUtil.h" +#include "itkComplexToModulusImageFilter.h." mitk::BandpassFilter::BandpassFilter() : m_HighPass(0), m_LowPass(50), m_HighPassAlpha(1), m_LowPassAlpha(1), m_FilterService(mitk::PhotoacousticFilterService::New()) { MITK_INFO << "Instantiating BandpassFilter..."; SetNumberOfIndexedInputs(1); SetNumberOfIndexedOutputs(1); MITK_INFO << "Instantiating BandpassFilter...[Done]"; } mitk::BandpassFilter::~BandpassFilter() { MITK_INFO << "Destructed BandpassFilter."; } itk::Image::Pointer BPFunction(mitk::Image::Pointer reference, float cutoffFrequencyPixelHighPass, float cutoffFrequencyPixelLowPass, float alphaHighPass, float alphaLowPass) { float* imageData = new float[reference->GetDimension(0)*reference->GetDimension(1)]; float width = cutoffFrequencyPixelLowPass - cutoffFrequencyPixelHighPass; - float center = cutoffFrequencyPixelHighPass / 2.0 + width / 2.0; + float center = cutoffFrequencyPixelHighPass + width / 2.f; for (unsigned int n = 0; n < reference->GetDimension(1); ++n) { imageData[reference->GetDimension(0)*n] = 0; } for (int n = 0; n < width; ++n) { - imageData[reference->GetDimension(0)*n] = 1; - if (n <= (alphaHighPass*(width - 1)) / 2.0) + imageData[reference->GetDimension(0)*(int)(n + center - (width / 2.f))] = 1; + if (n <= (alphaHighPass*(width - 1)) / 2.f) { - if (alphaHighPass > 0.00001) + if (alphaHighPass > 0.00001f) { - imageData[reference->GetDimension(0)*(int)(n + center - (width / 2))] = + imageData[reference->GetDimension(0)*(int)(n + center - (width / 2.f))] = (1 + cos(itk::Math::pi*(2 * n / (alphaHighPass*(width - 1)) - 1))) / 2; } else { - imageData[reference->GetDimension(0)*(int)(n + center - (width / 2))] = 1; + imageData[reference->GetDimension(0)*(int)(n + center - (width / 2.f))] = 1; } } - else if (n >= (width - 1)*(1 - alphaLowPass / 2)) + else if (n >= (width - 1)*(1 - alphaLowPass / 2.f)) { - if (alphaLowPass > 0.00001) + if (alphaLowPass > 0.00001f) { - imageData[reference->GetDimension(0)*(int)(n + center - (width / 2))] = + imageData[reference->GetDimension(0)*(int)(n + center - (width / 2.f))] = (1 + cos(itk::Math::pi*(2 * n / (alphaLowPass*(width - 1)) + 1 - 2 / alphaLowPass))) / 2; } else { - imageData[reference->GetDimension(0)*(int)(n + center - (width / 2))] = 1; + imageData[reference->GetDimension(0)*(int)(n + center - (width / 2.f))] = 1; } } } for (unsigned int n = reference->GetDimension(1) / 2; n < reference->GetDimension(1); ++n) { imageData[reference->GetDimension(0)*n] = imageData[(reference->GetDimension(1) - (n + 1)) * reference->GetDimension(0)]; } for (unsigned int line = 1; line < reference->GetDimension(0); ++line) { for (unsigned int sample = 0; sample < reference->GetDimension(1); ++sample) { imageData[reference->GetDimension(0)*sample + line] = imageData[reference->GetDimension(0)*sample]; } } typedef itk::Image< float, 3U > ImageType; ImageType::RegionType region; ImageType::IndexType start; start.Fill(0); region.SetIndex(start); ImageType::SizeType size; size[0] = reference->GetDimension(0); size[1] = reference->GetDimension(1); size[2] = reference->GetDimension(2); region.SetSize(size); ImageType::SpacingType SpacingItk; SpacingItk[0] = reference->GetGeometry()->GetSpacing()[0]; SpacingItk[1] = reference->GetGeometry()->GetSpacing()[1]; SpacingItk[2] = reference->GetGeometry()->GetSpacing()[2]; ImageType::Pointer image = ImageType::New(); image->SetRegions(region); image->Allocate(); image->FillBuffer(itk::NumericTraits::Zero); image->SetSpacing(SpacingItk); ImageType::IndexType pixelIndex; for (unsigned int slice = 0; slice < reference->GetDimension(2); ++slice) { for (unsigned int line = 0; line < reference->GetDimension(0); ++line) { for (unsigned int sample = 0; sample < reference->GetDimension(1); ++sample) { pixelIndex[0] = line; pixelIndex[1] = sample; pixelIndex[2] = slice; image->SetPixel(pixelIndex, imageData[line + sample*reference->GetDimension(0)]); } } } delete[] imageData; return image; } void mitk::BandpassFilter::GenerateData() { auto input = GetInput(); if (!(input->GetPixelType().GetTypeAsString() == "scalar (float)" || input->GetPixelType().GetTypeAsString() == " (float)")) { MITK_ERROR << "Pixel type is not float, abort"; mitkThrow() << "Pixel type is not float, abort"; } auto output = GetOutput(); - mitk::Image::Pointer resampledInput; + mitk::Image::Pointer resampledInput = input; bool powerOfTwo = false; int finalPower = 0; for (int i = 1; pow(2, i) <= input->GetDimension(1); ++i) { finalPower = i; if (pow(2, i) == input->GetDimension(1)) { powerOfTwo = true; } } if (!powerOfTwo) { double spacing_x = input->GetGeometry()->GetSpacing()[0]; double spacing_y = input->GetGeometry()->GetSpacing()[1] * (((double)input->GetDimensions()[1]) / pow(2, finalPower + 1)); double dim[2] = { spacing_x, spacing_y }; resampledInput = m_FilterService->ApplyResampling(input, dim); } // do a fourier transform, multiply with an appropriate window for the filter, and transform back typedef itk::Image< float, 3 > RealImageType; RealImageType::Pointer image; mitk::CastToItkImage(resampledInput, image); typedef itk::FFT1DRealToComplexConjugateImageFilter ForwardFFTFilterType; typedef ForwardFFTFilterType::OutputImageType ComplexImageType; ForwardFFTFilterType::Pointer forwardFFTFilter = ForwardFFTFilterType::New(); forwardFFTFilter->SetInput(image); forwardFFTFilter->SetDirection(1); try { forwardFFTFilter->UpdateOutputInformation(); } catch (itk::ExceptionObject & error) { std::cerr << "Error: " << error << std::endl; MITK_ERROR << "Bandpass could not be applied"; mitkThrow() << "bandpass error"; } + if (m_HighPass > m_LowPass) + mitkThrow() << "High pass frequency higher than low pass frequency, abort"; + + float singleVoxel = 2*M_PI / ((float)resampledInput->GetGeometry()->GetSpacing()[1] * resampledInput->GetDimension(1)); // [MHz] + float cutoffPixelHighPass = std::min((m_HighPass / singleVoxel), (float)resampledInput->GetDimension(1) / 2.0f); + float cutoffPixelLowPass = std::min((m_LowPass / singleVoxel), (float)resampledInput->GetDimension(1) / 2.0f); + + MITK_INFO << "SingleVoxel: " << singleVoxel; MITK_INFO << "HighPass: " << m_HighPass; MITK_INFO << "LowPass: " << m_LowPass; - float singleVoxel = 1.0f / ((float)resampledInput->GetGeometry()->GetSpacing()[1] / 1000); - float cutoffPixelHighPass = std::min((m_HighPass / singleVoxel) / 2, (float)resampledInput->GetDimension(1) / 2.0f); - float cutoffPixelLowPass = std::min((m_LowPass / singleVoxel) / 2, (float)resampledInput->GetDimension(1) / 2.0f - cutoffPixelHighPass); - MITK_INFO << "cutoffPixelHighPass: " << cutoffPixelHighPass; MITK_INFO << "cutoffPixelLowPass: " << cutoffPixelLowPass; RealImageType::Pointer fftMultiplicator = BPFunction(resampledInput, cutoffPixelHighPass, cutoffPixelLowPass, m_HighPassAlpha, m_LowPassAlpha); - mitk::Image::Pointer testImage; - mitk::CastToMitkImage(fftMultiplicator, testImage); typedef itk::MultiplyImageFilter< ComplexImageType, RealImageType, ComplexImageType > MultiplyFilterType; MultiplyFilterType::Pointer multiplyFilter = MultiplyFilterType::New(); multiplyFilter->SetInput1(forwardFFTFilter->GetOutput()); multiplyFilter->SetInput2(fftMultiplicator); + /* + GrabItkImageMemory(fftMultiplicator, output); + mitk::IOUtil::Save(output, "D:/fftImage.nrrd"); + + typedef itk::ComplexToModulusImageFilter< ComplexImageType, RealImageType > modulusType; + modulusType::Pointer modul = modulusType::New(); + + modul->SetInput(multiplyFilter->GetOutput()); + GrabItkImageMemory(modul->GetOutput(), output); + mitk::IOUtil::Save(output, "D:/fftout.nrrd"); + + modul->SetInput(forwardFFTFilter->GetOutput()); + GrabItkImageMemory(modul->GetOutput(), output); + mitk::IOUtil::Save(output, "D:/fftin.nrrd");*/ + typedef itk::FFT1DComplexConjugateToRealImageFilter< ComplexImageType, RealImageType > InverseFilterType; InverseFilterType::Pointer inverseFFTFilter = InverseFilterType::New(); inverseFFTFilter->SetInput(multiplyFilter->GetOutput()); inverseFFTFilter->SetDirection(1); GrabItkImageMemory(inverseFFTFilter->GetOutput(), output); } diff --git a/Modules/PhotoacousticsAlgorithms/test/mitkBandpassFilterTest.cpp b/Modules/PhotoacousticsAlgorithms/test/mitkBandpassFilterTest.cpp index 2d85d68c47..cf43d6efc3 100644 --- a/Modules/PhotoacousticsAlgorithms/test/mitkBandpassFilterTest.cpp +++ b/Modules/PhotoacousticsAlgorithms/test/mitkBandpassFilterTest.cpp @@ -1,145 +1,194 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ +#define _USE_MATH_DEFINES +#include #include #include #include #include #include -#include #include +#include + +#include "../ITKFilter/ITKUltrasound/itkFFT1DRealToComplexConjugateImageFilter.h" +#include "mitkImageCast.h" +#include "mitkITKImageImport.h" +#include "itkComplexToModulusImageFilter.h." class mitkBandpassFilterTestSuite : public mitk::TestFixture { CPPUNIT_TEST_SUITE(mitkBandpassFilterTestSuite); MITK_TEST(testHighPass); MITK_TEST(testLowPass); CPPUNIT_TEST_SUITE_END(); private: mitk::BandpassFilter::Pointer m_BandpassFilter; const unsigned int NUM_ITERATIONS = 15; - const unsigned int DATA_DIM = 500; - const float TIME_SPACING = 0.0001; // [us] - const float HIGHPASS_FREQENCY = 1.f / TIME_SPACING * DATA_DIM * 0.7f; // [Hz] - const float LOWPASS_FREQENCY = 1.f / TIME_SPACING * DATA_DIM * 0.3f; // [Hz] + const unsigned int DATA_XY_DIM = 512; + const unsigned int DATA_Z_DIM = 8; + const float TIME_SPACING = 0.00625; // [us] + const float FREQUENCY_RESOLUTION = 2 * M_PI / (TIME_SPACING * DATA_XY_DIM); // [MHz] + const float MAX_FREQUENCY = FREQUENCY_RESOLUTION * DATA_XY_DIM / 2.f; // [MHz] + const float HIGHPASS_FREQENCY = MAX_FREQUENCY * 0.8f; // [MHz] + const float LOWPASS_FREQENCY = MAX_FREQUENCY * 0.1f; // [MHz] const float ALPHA = 0; // 0 = box, 1 = von Hann; changing this may make the test invalid public: void setUp() override { m_BandpassFilter = mitk::BandpassFilter::New(); } - void test(float signalFreq, float HighPass, float LowPass, float HighPassAlpha, float LowPassAlpha) + void test(float HighPass, float LowPass, float HighPassAlpha, float LowPassAlpha, bool useLow, bool useHigh) { std::random_device r; std::default_random_engine randGen(r()); - std::uniform_real_distribution randDistrHighPass(0, HighPass * 0.9f); - std::uniform_real_distribution randDistrLowPass(LowPass * 1.1f, LowPass * 10); - - float* data = new float[DATA_DIM*DATA_DIM*DATA_DIM]; + std::uniform_real_distribution randDistrHighPass(HighPass * 0.01f, HighPass * 0.2f); + std::uniform_real_distribution randDistrLowPass(LowPass * 1.5f, LowPass * 2.f); + float* data = new float[DATA_XY_DIM*DATA_XY_DIM*DATA_Z_DIM]; mitk::Image::Pointer inputImage = mitk::Image::New(); - unsigned int dimension[3]{ DATA_DIM, DATA_DIM, DATA_DIM }; + unsigned int dimension[3]{ DATA_XY_DIM, DATA_XY_DIM, DATA_Z_DIM }; inputImage->Initialize(mitk::MakeScalarPixelType(), 3, dimension); mitk::Vector3D spacing; spacing[0] = 1; spacing[1] = TIME_SPACING; spacing[2] = 1; inputImage->SetSpacing(spacing); for (unsigned int iteration = 0; iteration < NUM_ITERATIONS; ++iteration) { // fill image with zeros - for (unsigned int i = 0; i < DATA_DIM*DATA_DIM*DATA_DIM; ++i) + for (unsigned int i = 0; i < DATA_XY_DIM*DATA_XY_DIM*DATA_Z_DIM; ++i) { data[i] = 0; } // write specific frequencies to the image - if (HighPass != -1) + if (useHigh) addFrequency(randDistrHighPass(randGen), TIME_SPACING, data, dimension); - if (LowPass != -1) + if (useLow) addFrequency(randDistrLowPass(randGen), TIME_SPACING, data, dimension); inputImage->SetImportVolume(data, 0, 0, mitk::Image::ImportMemoryManagementType::CopyMemory); m_BandpassFilter->SetInput(inputImage); + if (!useHigh) + HighPass = 0; m_BandpassFilter->SetHighPass(HighPass); + + if(!useLow) + LowPass = MAX_FREQUENCY; m_BandpassFilter->SetLowPass(LowPass); + m_BandpassFilter->SetHighPassAlpha(HighPassAlpha); m_BandpassFilter->SetLowPassAlpha(LowPassAlpha); m_BandpassFilter->Update(); mitk::Image::Pointer outputImage = m_BandpassFilter->GetOutput(); - mitk::ImageReadAccessor readAccess(outputImage); - const float* outputData = (const float*)readAccess.GetData(); + // do a fourier transform, and check whether the part of the image that has been filtered is zero + typedef itk::Image< float, 3 > RealImageType; + RealImageType::Pointer image; + + mitk::CastToItkImage(outputImage, image); + + typedef itk::FFT1DRealToComplexConjugateImageFilter ForwardFFTFilterType; + typedef ForwardFFTFilterType::OutputImageType ComplexImageType; + ForwardFFTFilterType::Pointer forwardFFTFilter = ForwardFFTFilterType::New(); + forwardFFTFilter->SetInput(image); + forwardFFTFilter->SetDirection(1); + forwardFFTFilter->UpdateOutputInformation(); + forwardFFTFilter->Update(); + + auto fftResult = forwardFFTFilter->GetOutput(); // the resulting image should consist only of zeros, as we filtered the frequencies out - for (unsigned int z = 0; z < DATA_DIM; ++z) + for (unsigned int z = 0; z < DATA_Z_DIM; ++z) { - for (unsigned int y = 0; y < DATA_DIM; ++y) + for (unsigned int y = 0; y < DATA_XY_DIM / 2; ++y) + { + if (y < (unsigned int)(HighPass / FREQUENCY_RESOLUTION) || y > (unsigned int)(LowPass / FREQUENCY_RESOLUTION)) + { + for (unsigned int x = 0; x < DATA_XY_DIM; ++x) + { + unsigned int outPos = x + y * DATA_XY_DIM + z * DATA_XY_DIM * DATA_XY_DIM; + std::complex value = fftResult->GetPixel({ x,y,z }); + 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)), + (abs(value.real()) < 0.1) && (abs(value.imag() < 0.1))); + } + } + } + + for (unsigned int y = DATA_XY_DIM / 2; y < DATA_XY_DIM; ++y) { - for (unsigned int x = 0; x < DATA_DIM; ++x) + if (y > DATA_XY_DIM - (unsigned int)(HighPass / FREQUENCY_RESOLUTION) || y < DATA_XY_DIM - (unsigned int)(LowPass / FREQUENCY_RESOLUTION)) { - unsigned int outPos = x + y * DATA_DIM + z * DATA_DIM * DATA_DIM; - CPPUNIT_ASSERT_MESSAGE(std::string("expected not :(, got " + std::to_string(outputData[outPos])), - abs(1) < mitk::eps); + for (unsigned int x = 0; x < DATA_XY_DIM; ++x) + { + unsigned int outPos = x + y * DATA_XY_DIM + z * DATA_XY_DIM * DATA_XY_DIM; + std::complex value = fftResult->GetPixel({ x,y,z }); + 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)), + (abs(value.real()) < 0.1) && (abs(value.imag() < 0.1))); + } } } } } delete[] data; } // write a fixed-frequency signal to the image void addFrequency(float freq, float timeSpacing, float* data, unsigned int* dim) { for (unsigned int z = 0; z < dim[2]; ++z) { - for (unsigned int y = 0; y < dim[2]; ++y) + for (unsigned int y = 0; y < dim[1]; ++y) { - for (unsigned int x = 0; x < dim[2]; ++x) + for (unsigned int x = 0; x < dim[0]; ++x) { data[x + y*dim[0] + z*dim[0] * dim[1]] += std::sin(freq * timeSpacing * y); } } } } void testHighPass() { - test(7.5, HIGHPASS_FREQENCY, -1, ALPHA, 0); + MITK_INFO << "Performing HighPass test"; + test(HIGHPASS_FREQENCY, 0, ALPHA, ALPHA, false, true); } void testLowPass() { - test(7.5, -1, LOWPASS_FREQENCY, 0, ALPHA); + return; + MITK_INFO << "Performing LowPass test"; + test(0, LOWPASS_FREQENCY, ALPHA, ALPHA, true, false); } void tearDown() override { m_BandpassFilter = nullptr; } }; MITK_TEST_SUITE_REGISTRATION(mitkBandpassFilter)