diff --git a/Modules/PhotoacousticsAlgorithms/test/mitkBandpassFilterTest.cpp b/Modules/PhotoacousticsAlgorithms/test/mitkBandpassFilterTest.cpp index 97dfebf107..06b1246bae 100644 --- a/Modules/PhotoacousticsAlgorithms/test/mitkBandpassFilterTest.cpp +++ b/Modules/PhotoacousticsAlgorithms/test/mitkBandpassFilterTest.cpp @@ -1,193 +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 "../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_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 + const float EPSILON_FFT = 0.00001f; public: void setUp() override { m_BandpassFilter = mitk::BandpassFilter::New(); } 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(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_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_XY_DIM*DATA_XY_DIM*DATA_Z_DIM; ++i) { data[i] = 0; } // write specific frequencies to the image if (useHigh) addFrequency(randDistrHighPass(randGen), TIME_SPACING, data, dimension); 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(); // 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_Z_DIM; ++z) { for (unsigned int y = 0; y < DATA_XY_DIM / 2; ++y) { if (y < (unsigned int)std::floor(HighPass / FREQUENCY_RESOLUTION) || y > (unsigned int)std::ceil(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.00001) && (abs(value.imag() < 0.00001))); + (abs(value.real()) < EPSILON_FFT) && (abs(value.imag() < EPSILON_FFT))); } } } for (unsigned int y = DATA_XY_DIM / 2; y < DATA_XY_DIM; ++y) { if (y > DATA_XY_DIM - (unsigned int)std::floor(HighPass / FREQUENCY_RESOLUTION) || y < DATA_XY_DIM - (unsigned int)std::ceil(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.00001) && (abs(value.imag() < 0.00001))); + (abs(value.real()) < EPSILON_FFT) && (abs(value.imag() < EPSILON_FFT))); } } } } } 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[1]; ++y) { 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() { MITK_INFO << "Performing HighPass test"; test(HIGHPASS_FREQENCY, 0, ALPHA, ALPHA, false, true); } void testLowPass() { MITK_INFO << "Performing LowPass test"; test(0, LOWPASS_FREQENCY, ALPHA, ALPHA, true, false); } void tearDown() override { m_BandpassFilter = nullptr; } }; MITK_TEST_SUITE_REGISTRATION(mitkBandpassFilter) diff --git a/Modules/PhotoacousticsAlgorithms/test/mitkBeamformingFilterTest.cpp b/Modules/PhotoacousticsAlgorithms/test/mitkBeamformingFilterTest.cpp index 0668999780..6e0dc027e7 100644 --- a/Modules/PhotoacousticsAlgorithms/test/mitkBeamformingFilterTest.cpp +++ b/Modules/PhotoacousticsAlgorithms/test/mitkBeamformingFilterTest.cpp @@ -1,82 +1,266 @@ /*=================================================================== 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. ===================================================================*/ #include #include #include #include #include +#include +#include +#include + +class SyntheticPAImageData +{ +public: + SyntheticPAImageData(float spacing_x, float spacing_y, unsigned int samples, unsigned int num_transducers, float speedOfSound) + { + m_Spacing_x = spacing_x; + m_Spacing_y = spacing_y; + m_TransducerElements = num_transducers; + m_Samples = samples; + m_SpeedOfSound = speedOfSound; + m_Data = new float[m_Samples*m_TransducerElements]; + } + + ~SyntheticPAImageData() + { + delete[] m_Data; + } + + void AddWave(float origin_depth, float origin_x, float base_value= 10000) + { + AddLine(origin_depth, origin_x, base_value); + AddLine(origin_depth + 0.0001, origin_x, -base_value); + } + + const float* GetData() + { + return (const float*)m_Data; + } + +private: + void AddLine(float origin_depth, unsigned int origin_x, float base_value = 10000) + { + for (unsigned int x = 0; x < m_TransducerElements; ++x) + { + float distance = std::abs((int)x - (int)origin_x); + float delay_in_seconds = std::sqrt(std::pow(origin_depth, 2) + std::pow(distance*m_Spacing_x, 2)) / m_SpeedOfSound; + unsigned int pixels = int(delay_in_seconds / m_Spacing_y); + + for (int index = -2; index < 5; ++index) + { + if ((int)pixels + index < (int)m_Samples && (int)pixels + index > 0) + { + //MITK_INFO << "yep: " << base_value / std::sqrt(distance + 1); + m_Data[x + (pixels + index)*m_TransducerElements] += base_value / std::sqrt(distance + 1); + } + } + } + } + + float* m_Data; + float m_Spacing_x; + float m_Spacing_y; + unsigned int m_TransducerElements; + unsigned int m_Samples; + float m_SpeedOfSound; +}; class mitkBeamformingFilterTestSuite : public mitk::TestFixture { CPPUNIT_TEST_SUITE(mitkBeamformingFilterTestSuite); - MITK_TEST(testBeamformingCPU); - MITK_TEST(testBeamformingGPU); + MITK_TEST(testBeamformingCPU_DAS); + MITK_TEST(testBeamformingGPU_DAS); CPPUNIT_TEST_SUITE_END(); private: mitk::BeamformingFilter::Pointer m_BeamformingFilter; const unsigned int NUM_ITERATIONS = 15; - const unsigned int DATA_DIM = 10; + const unsigned int SAMPLES = 5000; + const unsigned int RECONSTRUCTED_SAMPLES = 2048; + const unsigned int ELEMENTS = 128; + const unsigned int RECONSTRUCTED_LINES = 128; + const float SPEED_OF_SOUND = 1540; // m/s + const float SPACING_X = 0.3; // mm + const float SPACING_Y = 0.00625; // us public: void setUp() override { - m_BeamformingFilter = mitk::BeamformingFilter::New(nullptr); + } - void test(bool GPU) + void test() { - float* data = new float[DATA_DIM*DATA_DIM*DATA_DIM]; + std::random_device r; + std::default_random_engine randGen(r()); + float maxDepthInMeters = SAMPLES * (SPACING_Y / 1000000) * SPEED_OF_SOUND / 2; + std::uniform_real_distribution randDistrDepth(maxDepthInMeters*0.1, maxDepthInMeters*0.8); - for (unsigned int i = 0; i < DATA_DIM*DATA_DIM*DATA_DIM; ++i) + for (unsigned int iteration = 0; iteration < NUM_ITERATIONS; ++iteration) { - data[i] = 0; - } + // create some synthetic input data + float depth_in_meters1 = randDistrDepth(randGen); + float depth_in_meters2 = randDistrDepth(randGen); + float depth_in_meters3 = randDistrDepth(randGen); - mitk::Image::Pointer inputImage = mitk::Image::New(); - unsigned int dimension[3]{ DATA_DIM, DATA_DIM, DATA_DIM }; - inputImage->Initialize(mitk::MakeScalarPixelType(), 3, dimension); - inputImage->SetImportVolume(data); + unsigned int element1 = 29; + unsigned int element2 = 63; + unsigned int element3 = 98; - for (unsigned int iteration = 0; iteration < NUM_ITERATIONS; ++iteration) - { + SyntheticPAImageData image(SPACING_X, SPACING_Y, SAMPLES, ELEMENTS, SPEED_OF_SOUND); + image.AddWave(depth_in_meters1, 98); + image.AddWave(depth_in_meters2, 29); + image.AddWave(depth_in_meters3, 63); + + mitk::Image::Pointer inputImage = mitk::Image::New(); + unsigned int dimension[3]{ ELEMENTS, SAMPLES, 1 }; + inputImage->Initialize(mitk::MakeScalarPixelType(), 3, dimension); + mitk::Vector3D spacing; + spacing[0] = SPACING_X; + spacing[1] = SPACING_Y; + spacing[2] = 1; + inputImage->SetSpacing(spacing); + inputImage->SetImportVolume(image.GetData()); + mitk::IOUtil::Save(inputImage, "D:/dff.nrrd"); + + // setup the beamforming filter m_BeamformingFilter->SetInput(inputImage); m_BeamformingFilter->Update(); + mitk::Image::Pointer outputImage = m_BeamformingFilter->GetOutput(); mitk::ImageReadAccessor readAccess(outputImage); const float* outputData = (const float*)readAccess.GetData(); - for (unsigned int i = 0; i < DATA_DIM*DATA_DIM*DATA_DIM; ++i) + + unsigned int pos1[3] = { element1, (unsigned int)std::round(depth_in_meters1 * 1000 / outputImage->GetGeometry()->GetSpacing()[1]), 0 }; + unsigned int pos2[3] = { element2, (unsigned int)std::round(depth_in_meters2 * 1000 / outputImage->GetGeometry()->GetSpacing()[1]), 0 }; + unsigned int pos3[3] = { element3, (unsigned int)std::round(depth_in_meters3 * 1000 / outputImage->GetGeometry()->GetSpacing()[1]), 0 }; + + double average = 0; + + for (unsigned int i = 0; i < RECONSTRUCTED_LINES*RECONSTRUCTED_SAMPLES; ++i) { - CPPUNIT_ASSERT_MESSAGE(std::string("expected " + std::to_string(data[i]) + " but was " + std::to_string(outputData[i])), abs(outputData[i] - data[i]) < mitk::eps); + average += outputData[i] / RECONSTRUCTED_LINES*RECONSTRUCTED_SAMPLES; } + + MITK_INFO << "itssssssssssssss" << average; + CPPUNIT_ASSERT_MESSAGE(std::string("first point source incorrectly reconstructed"), abs(outputData[pos1[0] + pos1[1]*RECONSTRUCTED_LINES] / average) > 100); + CPPUNIT_ASSERT_MESSAGE(std::string("second point source incorrectly reconstructed"), abs(outputData[pos2[0] + pos2[1] * RECONSTRUCTED_LINES] / average) > 100); + CPPUNIT_ASSERT_MESSAGE(std::string("third point source incorrectly reconstructed"), abs(outputData[pos3[0] + pos3[1] * RECONSTRUCTED_LINES] / average) > 100); + } } - void testBeamformingCPU() + mitk::BeamformingSettings::Pointer createConfig(bool UseGPU, unsigned int* inputDim, mitk::BeamformingSettings::BeamformingAlgorithm alg) { - test(false); + return mitk::BeamformingSettings::New(SPACING_X / 1000, + SPEED_OF_SOUND, + SPACING_Y, + 27.f, + true, + RECONSTRUCTED_SAMPLES, + RECONSTRUCTED_LINES, + 0, + false, + nullptr, + inputDim, + UseGPU, + mitk::BeamformingSettings::DelayCalc::Spherical, + mitk::BeamformingSettings::Apodization::Box, + ELEMENTS * 2, + alg, + false, + 0, + 0); } - void testBeamformingGPU() + void testBeamformingCPU_DAS() { - test(true); + unsigned int* inputDim = new unsigned int[3]; + inputDim[0] = ELEMENTS; + inputDim[1] = RECONSTRUCTED_SAMPLES; + inputDim[2] = 1; + m_BeamformingFilter = mitk::BeamformingFilter::New(createConfig(false, inputDim, mitk::BeamformingSettings::BeamformingAlgorithm::DAS)); + + test(); + delete[] inputDim; + } + + void testBeamformingGPU_DAS() + { + unsigned int* inputDim = new unsigned int[3]; + inputDim[0] = ELEMENTS; + inputDim[1] = RECONSTRUCTED_SAMPLES; + inputDim[2] = 1; + m_BeamformingFilter = mitk::BeamformingFilter::New(createConfig(true, inputDim, mitk::BeamformingSettings::BeamformingAlgorithm::DAS)); + + test(); + delete[] inputDim; + } + + void testBeamformingCPU_sDMAS() + { + unsigned int* inputDim = new unsigned int[3]; + inputDim[0] = ELEMENTS; + inputDim[1] = RECONSTRUCTED_SAMPLES; + inputDim[2] = 1; + m_BeamformingFilter = mitk::BeamformingFilter::New(createConfig(false, inputDim, mitk::BeamformingSettings::BeamformingAlgorithm::sDMAS)); + + test(); + delete[] inputDim; + } + + void testBeamformingGPU_sDMAS() + { + unsigned int* inputDim = new unsigned int[3]; + inputDim[0] = ELEMENTS; + inputDim[1] = RECONSTRUCTED_SAMPLES; + inputDim[2] = 1; + m_BeamformingFilter = mitk::BeamformingFilter::New(createConfig(true, inputDim, mitk::BeamformingSettings::BeamformingAlgorithm::sDMAS)); + + test(); + delete[] inputDim; + } + + void testBeamformingCPU_DMAS() + { + unsigned int* inputDim = new unsigned int[3]; + inputDim[0] = ELEMENTS; + inputDim[1] = RECONSTRUCTED_SAMPLES; + inputDim[2] = 1; + m_BeamformingFilter = mitk::BeamformingFilter::New(createConfig(false, inputDim, mitk::BeamformingSettings::BeamformingAlgorithm::DMAS)); + + test(); + delete[] inputDim; + } + + void testBeamformingGPU_DMAS() + { + unsigned int* inputDim = new unsigned int[3]; + inputDim[0] = ELEMENTS; + inputDim[1] = RECONSTRUCTED_SAMPLES; + inputDim[2] = 1; + m_BeamformingFilter = mitk::BeamformingFilter::New(createConfig(true, inputDim, mitk::BeamformingSettings::BeamformingAlgorithm::DMAS)); + + test(); + delete[] inputDim; } }; MITK_TEST_SUITE_REGISTRATION(mitkBeamformingFilter) diff --git a/Modules/PhotoacousticsAlgorithms/test/mitkCropImageFilterTest.cpp b/Modules/PhotoacousticsAlgorithms/test/mitkCropImageFilterTest.cpp index f2d29e66c4..808382eb97 100644 --- a/Modules/PhotoacousticsAlgorithms/test/mitkCropImageFilterTest.cpp +++ b/Modules/PhotoacousticsAlgorithms/test/mitkCropImageFilterTest.cpp @@ -1,127 +1,127 @@ /*=================================================================== 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. ===================================================================*/ #include #include #include #include #include #include class mitkCropImageFilterTestSuite : public mitk::TestFixture { CPPUNIT_TEST_SUITE(mitkCropImageFilterTestSuite); MITK_TEST(testCropImage); CPPUNIT_TEST_SUITE_END(); private: mitk::CropImageFilter::Pointer m_CropImageFilter; const unsigned int NUM_ITERATIONS = 5; const unsigned int DATA_DIM = 25; public: void setUp() override { m_CropImageFilter = mitk::CropImageFilter::New(); } void test() { std::random_device r; std::default_random_engine randGen(r()); std::uniform_int_distribution randDistr(0, (DATA_DIM / 2) - 1); float* data = new float[DATA_DIM*DATA_DIM*DATA_DIM]; for (unsigned int i = 0; i < DATA_DIM*DATA_DIM*DATA_DIM; ++i) { data[i] = (float)i; } mitk::Image::Pointer inputImage = mitk::Image::New(); unsigned int dimension[3]{ DATA_DIM, DATA_DIM, DATA_DIM }; inputImage->Initialize(mitk::MakeScalarPixelType(), 3, dimension); inputImage->SetImportVolume(data, 0, 0, mitk::Image::ImportMemoryManagementType::CopyMemory); for (unsigned int iteration = 0; iteration < NUM_ITERATIONS; ++iteration) { unsigned int XPixelsCropStart = randDistr(randGen); unsigned int YPixelsCropStart = randDistr(randGen); unsigned int ZPixelsCropStart = randDistr(randGen); unsigned int XPixelsCropEnd = randDistr(randGen); unsigned int YPixelsCropEnd = randDistr(randGen); unsigned int ZPixelsCropEnd = randDistr(randGen); unsigned int newXDim = DATA_DIM - XPixelsCropStart - XPixelsCropEnd; unsigned int newYDim = DATA_DIM - YPixelsCropStart - YPixelsCropEnd; unsigned int newZDim = DATA_DIM - ZPixelsCropStart - ZPixelsCropEnd; m_CropImageFilter->SetInput(inputImage); m_CropImageFilter->SetXPixelsCropStart(XPixelsCropStart); m_CropImageFilter->SetYPixelsCropStart(YPixelsCropStart); m_CropImageFilter->SetZPixelsCropStart(ZPixelsCropStart); m_CropImageFilter->SetXPixelsCropEnd(XPixelsCropEnd); m_CropImageFilter->SetYPixelsCropEnd(YPixelsCropEnd); m_CropImageFilter->SetZPixelsCropEnd(ZPixelsCropEnd); m_CropImageFilter->Update(); mitk::Image::Pointer outputImage = m_CropImageFilter->GetOutput(); mitk::ImageReadAccessor readAccess(outputImage); const float* outputData = (const float*)readAccess.GetData(); CPPUNIT_ASSERT_MESSAGE(std::string("expected x size to be " + std::to_string(newXDim) + " but was " + std::to_string(outputImage->GetDimension(0))), newXDim == outputImage->GetDimension(0)); CPPUNIT_ASSERT_MESSAGE(std::string("expected y size to be " + std::to_string(newYDim) + " but was " + std::to_string(outputImage->GetDimension(1))), newYDim == outputImage->GetDimension(1)); CPPUNIT_ASSERT_MESSAGE(std::string("expected z size to be " + std::to_string(newZDim) + " but was " + std::to_string(outputImage->GetDimension(2))), newZDim == outputImage->GetDimension(2)); for (unsigned int z = 0; z < newZDim; ++z) { for (unsigned int y = 0; y < newYDim; ++y) { for (unsigned int x = 0; x < newXDim; ++x) { unsigned int origPos = (x + XPixelsCropStart) + (y + YPixelsCropStart) * DATA_DIM + (z + ZPixelsCropStart) * DATA_DIM * DATA_DIM; unsigned int outPos = x + y * newXDim + z * newXDim * newYDim; CPPUNIT_ASSERT_MESSAGE(std::string("expected " + std::to_string(data[origPos]) + " but was " + std::to_string(outputData[outPos])), - abs(outputData[outPos] - data[origPos]) < mitk::eps); + abs((float)outputData[outPos] - (float)data[origPos]) < mitk::eps); } } } } delete[] data; } void testCropImage() { for (int repetition = 0; repetition < 20; ++repetition) { MITK_INFO << "[" << (repetition + 1) << "/20]"; test(); } } void tearDown() override { m_CropImageFilter = nullptr; } }; MITK_TEST_SUITE_REGISTRATION(mitkCropImageFilter)