diff --git a/Modules/PhotoacousticsAlgorithms/include/OpenCLFilter/mitkPhotoacousticBModeFilter.h b/Modules/PhotoacousticsAlgorithms/include/OpenCLFilter/mitkPhotoacousticBModeFilter.h index 07413ed8e2..2f692e24fd 100644 --- a/Modules/PhotoacousticsAlgorithms/include/OpenCLFilter/mitkPhotoacousticBModeFilter.h +++ b/Modules/PhotoacousticsAlgorithms/include/OpenCLFilter/mitkPhotoacousticBModeFilter.h @@ -1,67 +1,67 @@ /*=================================================================== 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. ===================================================================*/ #ifndef _MITKPHOTOACOUSTICSBMODEFILTER_H_ #define _MITKPHOTOACOUSTICSBMODEFILTER_H_ #include #include "mitkImageToImageFilter.h" namespace mitk { /*! * \brief Class implementing a mitk::ImageToImageFilter for BMode filtering on CPU * * The only parameter that needs to be provided is whether it should use a logfilter. * Currently this class only performs an absolute BMode filter. */ class PhotoacousticBModeFilter : public ImageToImageFilter { public: mitkClassMacro(PhotoacousticBModeFilter, ImageToImageFilter); - itkFactorylessNewMacro(Self) - itkCloneMacro(Self) - - /** \brief Set parameters for the filter - * - * @param useLogFilter If true, the filter will apply a logfilter on the processed image - */ - void SetParameters(bool useLogFilter) + itkFactorylessNewMacro(Self); + itkCloneMacro(Self); + + /** \brief Set parameters for the filter + * + * @param useLogFilter If true, the filter will apply a logfilter on the processed image + */ + void UseLogFilter(bool useLogFilter) { m_UseLogFilter = useLogFilter; } protected: PhotoacousticBModeFilter(); ~PhotoacousticBModeFilter() override; void GenerateInputRequestedRegion() override; void GenerateOutputInformation() override; void GenerateData() override; //##Description //## @brief Time when Header was last initialized itk::TimeStamp m_TimeOfHeaderInitialization; bool m_UseLogFilter; }; } #endif /* _MITKPHOTOACOUSTICSBMODEFILTER_H_ */ diff --git a/Modules/PhotoacousticsAlgorithms/include/mitkPhotoacousticFilterService.h b/Modules/PhotoacousticsAlgorithms/include/mitkPhotoacousticFilterService.h index d64ff8f64d..fa5d0b81be 100644 --- a/Modules/PhotoacousticsAlgorithms/include/mitkPhotoacousticFilterService.h +++ b/Modules/PhotoacousticsAlgorithms/include/mitkPhotoacousticFilterService.h @@ -1,132 +1,131 @@ /*=================================================================== 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. ===================================================================*/ #ifndef mitkPhotoacousticFilterService_H_HEADER_INCLUDED #define mitkPhotoacousticFilterService_H_HEADER_INCLUDED #include "itkObject.h" #include "mitkCommon.h" #include "mitkImage.h" #include #include "mitkBeamformingSettings.h" #include "mitkBeamformingFilter.h" #include "MitkPhotoacousticsAlgorithmsExports.h" namespace mitk { /*! * \brief Class holding methods to apply all Filters within the Photoacoustics Algorithms Module * * Implemented are: * - A B-Mode Filter * - A Resampling Filter * - Beamforming on GPU and CPU * - A Bandpass Filter */ class MITKPHOTOACOUSTICSALGORITHMS_EXPORT PhotoacousticFilterService : public itk::Object { public: mitkClassMacroItkParent(mitk::PhotoacousticFilterService, itk::Object); itkFactorylessNewMacro(Self); /** \brief Defines the methods for the B-Mode filter * Currently implemented are an Envelope Detection filter and a simple Absolute filter. */ enum BModeMethod { EnvelopeDetection, Abs }; /** \brief Applies a B-Mode Filter * * Applies a B-Mode filter using the given parameters. * @param inputImage The image to be processed. * @param method The kind of B-Mode Filter to be used. * @param UseLogFilter Setting this to true will apply a simple logarithm to the image after the B-Mode Filter has been applied. * @param resampleSpacing If this is set to 0, nothing will be done; otherwise, the image is resampled to a spacing of resampleSpacing mm per pixel. * @return The processed image is returned after the filter has finished. */ mitk::Image::Pointer ApplyBmodeFilter(mitk::Image::Pointer inputImage, BModeMethod method = BModeMethod::Abs, - bool UseLogFilter = false, - float resampleSpacing = 0.15); + bool UseLogFilter = false); /** \brief Resamples the given image * * Resamples an image using the given parameters. * @param inputImage The image to be processed. * @param outputSize An array of dimensions the image should be resampled to. * @return The processed image is returned after the filter has finished. */ - mitk::Image::Pointer ApplyResampling(mitk::Image::Pointer inputImage, unsigned int outputSize[2]); + mitk::Image::Pointer ApplyResampling(mitk::Image::Pointer inputImage, double outputSpacing[2]); /** \brief Beamforms the given image * * Resamples an image using the given parameters. * @param inputImage The image to be processed. * @param config The configuration set to be used for beamforming. * @param progressHandle An std::function, through which progress of the currently updating filter is reported. * The integer argument is a number between 0 an 100 to indicate how far completion has been achieved, the std::string argument indicates what the filter is currently doing. * @return The processed image is returned after the filter has finished. */ mitk::Image::Pointer ApplyBeamforming(mitk::Image::Pointer inputImage, BeamformingSettings::Pointer config, std::function progressHandle = [](int, std::string) {}); /** \brief Crops the given image * * Crops an image in 3 dimension using the given parameters. * @param inputImage The image to be processed. * @param above How many voxels will be cut from the top of the image. * @param below How many voxels will be cut from the bottom of the image. * @param right How many voxels will be cut from the right side of the image. * @param left How many voxels will be cut from the left side of the image. * @param minSlice The first slice to be present in the resulting image. * @param maxSlice The last slice to be present in the resulting image. * @return The processed image is returned after the filter has finished. For the purposes of this module, the returned image is always of type float. */ mitk::Image::Pointer ApplyCropping(mitk::Image::Pointer inputImage, int above, int below, int right, int left, int minSlice, int maxSlice); /** \brief Applies a Bandpass filter to the given image * * Applies a bandpass filter to the given image using the given parameters. * @param data The image to be processed. * @param recordTime The depth of the image in seconds. * @param BPHighPass The position at which Lower frequencies are completely cut off in Hz. * @param BPLowPass The position at which Higher frequencies are completely cut off in Hz. * @param alphaHighPass The high pass tukey window parameter to control the shape of the bandpass filter: 0 will make it a Box function, 1 a Hann function. alpha can be set between those two bounds. * @param alphaLowPass The low passtukey window parameter to control the shape of the bandpass filter: 0 will make it a Box function, 1 a Hann function. alpha can be set between those two bounds. * @return The processed image is returned after the filter has finished. */ mitk::Image::Pointer BandpassFilter(mitk::Image::Pointer data, float recordTime, float BPHighPass, float BPLowPass, float alphaHighPass, float alphaLowPass); protected: PhotoacousticFilterService(); ~PhotoacousticFilterService() override; /** \brief For performance reasons, an instance of the Beamforming filter is initialized as soon as possible and kept for all further uses. */ mitk::BeamformingFilter::Pointer m_BeamformingFilter; /** \brief Function that creates a Tukey function for the bandpass */ itk::Image::Pointer BPFunction(mitk::Image::Pointer reference, - int cutoffFrequencyPixelHighPass, - int cutoffFrequencyPixelLowPass, + float cutoffFrequencyPixelHighPass, + float cutoffFrequencyPixelLowPass, float alphaHighPass, float alphaLowPass); }; } // namespace mitk #endif /* mitkPhotoacousticFilterService_H_HEADER_INCLUDED */ diff --git a/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticBModeFilter.cpp b/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticBModeFilter.cpp index b2267d3a97..8408229d5a 100644 --- a/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticBModeFilter.cpp +++ b/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticBModeFilter.cpp @@ -1,94 +1,91 @@ /*=================================================================== 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 "./OpenCLFilter/mitkPhotoacousticBModeFilter.h" #include "usServiceReference.h" #include mitk::PhotoacousticBModeFilter::PhotoacousticBModeFilter() : m_UseLogFilter(false) { this->SetNumberOfIndexedInputs(1); this->SetNumberOfRequiredInputs(1); } mitk::PhotoacousticBModeFilter::~PhotoacousticBModeFilter() { } void mitk::PhotoacousticBModeFilter::GenerateInputRequestedRegion() { Superclass::GenerateInputRequestedRegion(); - mitk::Image* output = this->GetOutput(); - mitk::Image* input = const_cast (this->GetInput()); + mitk::Image::Pointer output = this->GetOutput(); + mitk::Image::Pointer input = this->GetInput(); + if (!output->IsInitialized()) - { return; - } input->SetRequestedRegionToLargestPossibleRegion(); } void mitk::PhotoacousticBModeFilter::GenerateOutputInformation() { mitk::Image::ConstPointer input = this->GetInput(); mitk::Image::Pointer output = this->GetOutput(); if ((output->IsInitialized()) && (this->GetMTime() <= m_TimeOfHeaderInitialization.GetMTime())) return; - itkDebugMacro(<< "GenerateOutputInformation()"); - output->Initialize(input->GetPixelType(), input->GetDimension(), input->GetDimensions()); output->GetGeometry()->SetSpacing(input->GetGeometry()->GetSpacing()); output->GetGeometry()->Modified(); output->SetPropertyList(input->GetPropertyList()->Clone()); m_TimeOfHeaderInitialization.Modified(); } void mitk::PhotoacousticBModeFilter::GenerateData() { GenerateOutputInformation(); mitk::Image::Pointer input = this->GetInput(); mitk::Image::Pointer output = this->GetOutput(); if (!output->IsInitialized()) return; mitk::ImageReadAccessor reader(input); unsigned int size = output->GetDimension(0) * output->GetDimension(1) * output->GetDimension(2); const float* InputData = (const float*)(reader.GetData()); float* OutputData = new float[size]; if (!m_UseLogFilter) for (unsigned int i = 0; i < size; ++i) { OutputData[i] = abs(InputData[i]); } else { for (unsigned int i = 0; i < size; ++i) { OutputData[i] = log(abs(InputData[i])); } } output->SetImportVolume(OutputData, 0, 0, mitk::Image::ImportMemoryManagementType::CopyMemory); delete[] OutputData; m_TimeOfHeaderInitialization.Modified(); } diff --git a/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.cpp b/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.cpp index 6a5eb4a76c..588047ab63 100644 --- a/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.cpp +++ b/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.cpp @@ -1,260 +1,261 @@ /*=================================================================== 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. ===================================================================*/ #if defined(PHOTOACOUSTICS_USE_GPU) || DOXYGEN #include "./OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.h" #include "usServiceReference.h" mitk::PhotoacousticOCLBeamformingFilter::PhotoacousticOCLBeamformingFilter(BeamformingSettings::Pointer settings) : m_PixelCalculation(NULL), m_InputImage(mitk::Image::New()), m_ApodizationBuffer(nullptr), m_MemoryLocationsBuffer(nullptr), m_DelaysBuffer(nullptr), m_UsedLinesBuffer(nullptr), m_Conf(settings) { MITK_INFO << "Initializing OCL beamforming Filter..."; this->AddSourceFile("DAS.cl"); this->AddSourceFile("DMAS.cl"); this->AddSourceFile("sDMAS.cl"); this->m_FilterID = "OpenCLBeamformingFilter"; this->Initialize(); unsigned int dim[] = { 128, 2048, 2 }; m_InputImage->Initialize(mitk::MakeScalarPixelType(), 3, dim); m_ChunkSize[0] = 128; m_ChunkSize[1] = 128; m_ChunkSize[2] = 8; m_UsedLinesCalculation = mitk::OCLUsedLinesCalculation::New(m_Conf); m_DelayCalculation = mitk::OCLDelayCalculation::New(m_Conf); MITK_INFO << "Initializing OCL beamforming Filter...[Done]"; } mitk::PhotoacousticOCLBeamformingFilter::~PhotoacousticOCLBeamformingFilter() { if (this->m_PixelCalculation) { clReleaseKernel(m_PixelCalculation); } if (m_ApodizationBuffer) clReleaseMemObject(m_ApodizationBuffer); } void mitk::PhotoacousticOCLBeamformingFilter::Update() { //Check if context & program available if (!this->Initialize()) { us::ServiceReference ref = GetModuleContext()->GetServiceReference(); OclResourceService* resources = GetModuleContext()->GetService(ref); // clean-up also the resources resources->InvalidateStorage(); mitkThrow() << "Filter is not initialized. Cannot update."; } else { // Execute this->Execute(); } } void mitk::PhotoacousticOCLBeamformingFilter::UpdateDataBuffers() { /*us::ServiceReference ref = GetModuleContext()->GetServiceReference(); OclResourceService* resources = GetModuleContext()->GetService(ref); cl_ulong globalMemSize = oclGetGlobalMemSize(resources->GetCurrentDevice());*/ //Initialize the Output try { MITK_INFO << "Updating Workgroup size for new dimensions"; size_t outputSize = (size_t)m_Conf->GetReconstructionLines() * (size_t)m_Conf->GetSamplesPerLine() * (size_t)m_Conf->GetInputDim()[2]; m_OutputDim[0] = m_Conf->GetReconstructionLines(); m_OutputDim[1] = m_Conf->GetSamplesPerLine(); m_OutputDim[2] = m_Conf->GetInputDim()[2]; this->InitExec(this->m_PixelCalculation, m_OutputDim, outputSize, sizeof(float)); } catch (const mitk::Exception& e) { MITK_ERROR << "Caught exception while initializing filter: " << e.what(); return; } //TODO FIXME cl_int clErr = 0; MITK_DEBUG << "Updating GPU Buffers for new configuration"; // create the apodisation buffer if (m_Apodisation == nullptr) { MITK_INFO << "No apodisation function set; Beamforming will be done without any apodisation."; m_Apodisation = new float[1]{ 1 }; m_ApodArraySize = 1; } us::ServiceReference ref = GetModuleContext()->GetServiceReference(); OclResourceService* resources = GetModuleContext()->GetService(ref); cl_context gpuContext = resources->GetContext(); if (m_ApodizationBuffer) clReleaseMemObject(m_ApodizationBuffer); this->m_ApodizationBuffer = clCreateBuffer(gpuContext, CL_MEM_READ_ONLY | CL_MEM_USE_HOST_PTR, sizeof(float) * m_ApodArraySize, const_cast(m_Apodisation), &clErr); CHECK_OCL_ERR(clErr); // calculate used lines m_UsedLinesCalculation->Update(); m_UsedLinesBuffer = m_UsedLinesCalculation->GetGPUOutput()->GetGPUBuffer(); // calculate the Delays m_DelayCalculation->SetInputs(m_UsedLinesBuffer); m_DelayCalculation->Update(); m_DelaysBuffer = m_DelayCalculation->GetGPUOutput()->GetGPUBuffer(); //m_ConfOld = m_Conf; } void mitk::PhotoacousticOCLBeamformingFilter::Execute() { cl_int clErr = 0; UpdateDataBuffers(); unsigned int reconstructionLines = this->m_Conf->GetReconstructionLines(); unsigned int samplesPerLine = this->m_Conf->GetSamplesPerLine(); clErr = clSetKernelArg(this->m_PixelCalculation, 2, sizeof(cl_mem), &(this->m_UsedLinesBuffer)); clErr |= clSetKernelArg(this->m_PixelCalculation, 3, sizeof(cl_mem), &(this->m_DelaysBuffer)); clErr |= clSetKernelArg(this->m_PixelCalculation, 4, sizeof(cl_mem), &(this->m_ApodizationBuffer)); clErr |= clSetKernelArg(this->m_PixelCalculation, 5, sizeof(cl_ushort), &(this->m_ApodArraySize)); clErr |= clSetKernelArg(this->m_PixelCalculation, 6, sizeof(cl_uint), &(this->m_Conf->GetInputDim()[0])); clErr |= clSetKernelArg(this->m_PixelCalculation, 7, sizeof(cl_uint), &(this->m_Conf->GetInputDim()[1])); clErr |= clSetKernelArg(this->m_PixelCalculation, 8, sizeof(cl_uint), &(this->m_Conf->GetInputDim()[2])); clErr |= clSetKernelArg(this->m_PixelCalculation, 9, sizeof(cl_uint), &(reconstructionLines)); clErr |= clSetKernelArg(this->m_PixelCalculation, 10, sizeof(cl_uint), &(samplesPerLine)); // execute the filter on a 3D NDRange if (m_OutputDim[2] == 1 || m_ChunkSize[2] == 1) { if (!this->ExecuteKernelChunksInBatches(m_PixelCalculation, 2, m_ChunkSize, 16, 50)) mitkThrow() << "openCL Error when executing Kernel"; } else { if (!this->ExecuteKernelChunksInBatches(m_PixelCalculation, 3, m_ChunkSize, 16, 50)) mitkThrow() << "openCL Error when executing Kernel"; } // signalize the GPU-side data changed m_Output->Modified(GPU_DATA); } us::Module *mitk::PhotoacousticOCLBeamformingFilter::GetModule() { return us::GetModuleContext()->GetModule(); } bool mitk::PhotoacousticOCLBeamformingFilter::Initialize() { MITK_INFO << "Initializing OCL Beamforming Filter..."; bool buildErr = true; cl_int clErr = 0; if (OclFilter::Initialize()) { switch (m_Conf->GetAlgorithm()) { case BeamformingSettings::BeamformingAlgorithm::DAS: { this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "ckDAS", &clErr); break; } case BeamformingSettings::BeamformingAlgorithm::DMAS: { this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "ckDMAS", &clErr); break; } case BeamformingSettings::BeamformingAlgorithm::sDMAS: { this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "cksDMAS", &clErr); break; } default: { MITK_INFO << "No beamforming algorithm specified, setting to DAS"; this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "ckDAS", &clErr); break; } } buildErr |= CHECK_OCL_ERR(clErr); } CHECK_OCL_ERR(clErr); MITK_INFO << "Initializing OCL Beamforming Filter...[Done]"; return (OclFilter::IsInitialized() && buildErr); } void mitk::PhotoacousticOCLBeamformingFilter::SetInput(mitk::Image::Pointer image) { OclDataSetToDataSetFilter::SetInput(image); m_InputImage = image; } void mitk::PhotoacousticOCLBeamformingFilter::SetInput(void* data, unsigned int* dimensions, unsigned int BpE) { OclDataSetToDataSetFilter::SetInput(data, dimensions[0] * dimensions[1] * dimensions[2], BpE); } mitk::Image::Pointer mitk::PhotoacousticOCLBeamformingFilter::GetOutputAsImage() { mitk::Image::Pointer outputImage = mitk::Image::New(); if (m_Output->IsModified(GPU_DATA)) { void* pData = m_Output->TransferDataToCPU(m_CommandQue); const unsigned int dimension = 3; unsigned int dimensions[3] = { (unsigned int)m_OutputDim[0], (unsigned int)m_OutputDim[1], (unsigned int)m_OutputDim[2] }; const mitk::SlicedGeometry3D::Pointer p_slg = m_InputImage->GetSlicedGeometry(); MITK_DEBUG << "Creating new MITK Image."; outputImage->Initialize(this->GetOutputType(), dimension, dimensions); outputImage->SetSpacing(p_slg->GetSpacing()); - outputImage->SetImportVolume(pData, 0, 0, mitk::Image::ImportMemoryManagementType::ManageMemory); + outputImage->SetImportVolume(pData, 0, 0, mitk::Image::ImportMemoryManagementType::CopyMemory); + delete[] pData; } MITK_DEBUG << "Image Initialized."; return outputImage; } void* mitk::PhotoacousticOCLBeamformingFilter::GetOutput() { return OclDataSetToDataSetFilter::GetOutput(); } #endif diff --git a/Modules/PhotoacousticsAlgorithms/source/filters/mitkCastToFloatImageFilter.cpp b/Modules/PhotoacousticsAlgorithms/source/filters/mitkCastToFloatImageFilter.cpp index 4809971115..cc26284be8 100644 --- a/Modules/PhotoacousticsAlgorithms/source/filters/mitkCastToFloatImageFilter.cpp +++ b/Modules/PhotoacousticsAlgorithms/source/filters/mitkCastToFloatImageFilter.cpp @@ -1,92 +1,92 @@ /*=================================================================== 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. ===================================================================*/ #include "mitkCastToFloatImageFilter.h" #include "mitkImageReadAccessor.h" #include "mitkImageWriteAccessor.h" mitk::CastToFloatImageFilter::CastToFloatImageFilter() { MITK_INFO << "Instantiating CastToFloatImageFilter..."; SetNumberOfIndexedInputs(1); SetNumberOfIndexedOutputs(1); MITK_INFO << "Instantiating CastToFloatImageFilter...[Done]"; } mitk::CastToFloatImageFilter::~CastToFloatImageFilter() { MITK_INFO << "Destructed CastToFloatImageFilter."; } template float* CastData(const void* inputArray, unsigned long length) { float* outputArray = new float[length]; for (unsigned long i = 0; i < length; ++i) { outputArray[i] = (float)((TPixelType*)inputArray)[i]; } return outputArray; } void mitk::CastToFloatImageFilter::GenerateData() { mitk::Image::Pointer inputImage = GetInput(); mitk::Image::Pointer outputImage = GetOutput(); std::string type = inputImage->GetPixelType().GetTypeAsString(); if (type == "scalar (float)" || type == " (float)") { outputImage->Initialize(mitk::MakeScalarPixelType(), inputImage->GetDimension(), inputImage->GetDimensions()); outputImage->SetSpacing(inputImage->GetGeometry()->GetSpacing()); - outputImage->SetImportVolume(mitk::ImageWriteAccessor(inputImage).GetData(), 0, 0, mitk::Image::ImportMemoryManagementType::RtlCopyMemory); + outputImage->SetImportVolume(mitk::ImageWriteAccessor(inputImage).GetData(), 0, 0, mitk::Image::ImportMemoryManagementType::CopyMemory); MITK_INFO << "Input is already float type. Nothing to do here."; return; } if (outputImage.IsNull()) { outputImage = mitk::Image::New(); } unsigned long length = 1; for (unsigned int i = 0, limit = inputImage->GetDimension(); i < limit; ++i) length *= inputImage->GetDimensions()[i]; float* outputData; mitk::ImageReadAccessor readAccess(inputImage); if (type == "scalar (short)" || type == " (short)") outputData = CastData(readAccess.GetData(), length); else if (type == "scalar (unsigned short)" || type == " (unsigned short)" || type == " (unsigned_short)" || type == "scalar (unsigned_short)") outputData = CastData(readAccess.GetData(), length); else if (type == "scalar (int)" || type == " (int)") outputData = CastData(readAccess.GetData(), length); else if (type == "scalar (unsigned int)" || type == " (unsigned int)" || type == " (unsigned_int)" || type == "scalar (unsigned_int)") outputData = CastData(readAccess.GetData(), length); else if (type == "scalar (double)" || type == " (double)") outputData = CastData(readAccess.GetData(), length); else if (type == "scalar (long)" || type == " (long)") outputData = CastData(readAccess.GetData(), length); else mitkThrow() << "The given image has a datatype that is not yet supported. Given datatype: " << type; outputImage->Initialize(mitk::MakeScalarPixelType(), inputImage->GetDimension(), inputImage->GetDimensions()); outputImage->SetSpacing(inputImage->GetGeometry()->GetSpacing()); - outputImage->SetImportVolume(outputData, 0, 0, mitk::Image::ImportMemoryManagementType::RtlCopyMemory); + outputImage->SetImportVolume(outputData, 0, 0, mitk::Image::ImportMemoryManagementType::CopyMemory); delete[] outputData; } diff --git a/Modules/PhotoacousticsAlgorithms/source/utils/mitkPhotoacousticFilterService.cpp b/Modules/PhotoacousticsAlgorithms/source/utils/mitkPhotoacousticFilterService.cpp index 3fa2126eb8..8c2207663f 100644 --- a/Modules/PhotoacousticsAlgorithms/source/utils/mitkPhotoacousticFilterService.cpp +++ b/Modules/PhotoacousticsAlgorithms/source/utils/mitkPhotoacousticFilterService.cpp @@ -1,400 +1,364 @@ /*=================================================================== 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 "mitkPhotoacousticFilterService.h" #include "../ITKFilter/ITKUltrasound/itkBModeImageFilter.h" #include "../ITKFilter/itkPhotoacousticBModeImageFilter.h" #include "mitkImageCast.h" #include "mitkITKImageImport.h" #include "mitkBeamformingFilter.h" #include #include #include #include "./OpenCLFilter/mitkPhotoacousticBModeFilter.h" #include "mitkConvert2Dto3DImageFilter.h" #include // itk dependencies #include "itkImage.h" #include "itkResampleImageFilter.h" #include "itkCastImageFilter.h" #include "itkCropImageFilter.h" #include "itkRescaleIntensityImageFilter.h" #include "itkIntensityWindowingImageFilter.h" #include #include "itkMultiplyImageFilter.h" #include "itkBSplineInterpolateImageFunction.h" #include // needed itk image filters #include "mitkITKImageImport.h" #include "itkFFTShiftImageFilter.h" #include "itkMultiplyImageFilter.h" #include "itkComplexToModulusImageFilter.h" #include #include "../ITKFilter/ITKUltrasound/itkFFT1DComplexConjugateToRealImageFilter.h" #include "../ITKFilter/ITKUltrasound/itkFFT1DRealToComplexConjugateImageFilter.h" mitk::PhotoacousticFilterService::PhotoacousticFilterService() { MITK_INFO << "[PhotoacousticFilterService] created filter service"; } mitk::PhotoacousticFilterService::~PhotoacousticFilterService() { MITK_INFO << "[PhotoacousticFilterService] destructed filter service"; } -mitk::Image::Pointer mitk::PhotoacousticFilterService::ApplyBmodeFilter(mitk::Image::Pointer inputImage, BModeMethod method, bool UseLogFilter, float resampleSpacing) +mitk::Image::Pointer mitk::PhotoacousticFilterService::ApplyBmodeFilter(mitk::Image::Pointer inputImage, BModeMethod method, bool UseLogFilter) { // the image needs to be of floating point type for the envelope filter to work; the casting is done automatically by the CastToItkImage typedef itk::Image< float, 3 > itkFloatImageType; - if (inputImage.IsNull() || !(inputImage->GetPixelType().GetTypeAsString() == "scalar (float)" || inputImage->GetPixelType().GetTypeAsString() == " (float)")) + if (inputImage.IsNull() || + !(inputImage->GetPixelType().GetTypeAsString() == "scalar (float)" + || inputImage->GetPixelType().GetTypeAsString() == " (float)")) { MITK_ERROR << "BMode Filter can only handle float image types."; mitkThrow() << "BMode Filter can only handle float image types."; } - mitk::Image::Pointer out; - if (method == BModeMethod::Abs) { PhotoacousticBModeFilter::Pointer filter = PhotoacousticBModeFilter::New(); - filter->SetParameters(UseLogFilter); + filter->UseLogFilter(UseLogFilter); filter->SetInput(inputImage); filter->Update(); - out = filter->GetOutput(); + return filter->GetOutput(); } - else if (method == BModeMethod::EnvelopeDetection) - { - typedef itk::BModeImageFilter < itkFloatImageType, itkFloatImageType > BModeFilterType; - BModeFilterType::Pointer bModeFilter = BModeFilterType::New(); // LogFilter - - typedef itk::PhotoacousticBModeImageFilter < itkFloatImageType, itkFloatImageType > PhotoacousticBModeImageFilter; - PhotoacousticBModeImageFilter::Pointer photoacousticBModeFilter = PhotoacousticBModeImageFilter::New(); // No LogFilter - - typedef itk::ResampleImageFilter < itkFloatImageType, itkFloatImageType > ResampleImageFilter; - ResampleImageFilter::Pointer resampleImageFilter = ResampleImageFilter::New(); - itkFloatImageType::Pointer itkImage; + typedef itk::BModeImageFilter < itkFloatImageType, itkFloatImageType > BModeFilterType; + BModeFilterType::Pointer bModeFilter = BModeFilterType::New(); // LogFilter - mitk::CastToItkImage(inputImage, itkImage); - - itkFloatImageType::Pointer bmode; - - if (UseLogFilter) - { - bModeFilter->SetInput(itkImage); - bModeFilter->SetDirection(1); - bmode = bModeFilter->GetOutput(); - } - else - { - photoacousticBModeFilter->SetInput(itkImage); - photoacousticBModeFilter->SetDirection(1); - bmode = photoacousticBModeFilter->GetOutput(); - } - - out = mitk::GrabItkImageMemory(bmode); - } - - if (resampleSpacing == 0) - return out; + typedef itk::PhotoacousticBModeImageFilter < itkFloatImageType, itkFloatImageType > PhotoacousticBModeImageFilter; + PhotoacousticBModeImageFilter::Pointer photoacousticBModeFilter = PhotoacousticBModeImageFilter::New(); // No LogFilter typedef itk::ResampleImageFilter < itkFloatImageType, itkFloatImageType > ResampleImageFilter; ResampleImageFilter::Pointer resampleImageFilter = ResampleImageFilter::New(); - itkFloatImageType::Pointer itkImage = itkFloatImageType::New(); - mitk::CastToItkImage(out, itkImage); - itkFloatImageType::SpacingType outputSpacing; - itkFloatImageType::SizeType inputSize = itkImage->GetLargestPossibleRegion().GetSize(); - itkFloatImageType::SizeType outputSize = inputSize; - - outputSpacing[0] = itkImage->GetSpacing()[0]; - outputSpacing[1] = resampleSpacing; - outputSpacing[2] = itkImage->GetSpacing()[2]; + itkFloatImageType::Pointer itkImage; - outputSize[1] = inputSize[1] * itkImage->GetSpacing()[1] / outputSpacing[1]; + mitk::CastToItkImage(inputImage, itkImage); - typedef itk::IdentityTransform TransformType; - resampleImageFilter->SetInput(itkImage); - resampleImageFilter->SetSize(outputSize); - resampleImageFilter->SetOutputSpacing(outputSpacing); - resampleImageFilter->SetTransform(TransformType::New()); + if (UseLogFilter) + { + bModeFilter->SetInput(itkImage); + bModeFilter->SetDirection(1); + itkImage = bModeFilter->GetOutput(); + } + else + { + photoacousticBModeFilter->SetInput(itkImage); + photoacousticBModeFilter->SetDirection(1); + itkImage = photoacousticBModeFilter->GetOutput(); + } - resampleImageFilter->UpdateLargestPossibleRegion(); - return mitk::GrabItkImageMemory(resampleImageFilter->GetOutput()); + return mitk::GrabItkImageMemory(itkImage); } -mitk::Image::Pointer mitk::PhotoacousticFilterService::ApplyResampling(mitk::Image::Pointer inputImage, unsigned int outputSize[2]) +mitk::Image::Pointer mitk::PhotoacousticFilterService::ApplyResampling( + mitk::Image::Pointer inputImage, + double outputSpacing[2]) { typedef itk::Image< float, 3 > itkFloatImageType; typedef itk::ResampleImageFilter < itkFloatImageType, itkFloatImageType > ResampleImageFilter; ResampleImageFilter::Pointer resampleImageFilter = ResampleImageFilter::New(); - typedef itk::LinearInterpolateImageFunction T_Interpolator; + itkFloatImageType::Pointer itkImage; mitk::CastToItkImage(inputImage, itkImage); itkFloatImageType::SpacingType outputSpacingItk; itkFloatImageType::SizeType inputSizeItk = itkImage->GetLargestPossibleRegion().GetSize(); itkFloatImageType::SizeType outputSizeItk = inputSizeItk; - outputSizeItk[0] = outputSize[0]; - outputSizeItk[1] = outputSize[1]; - outputSizeItk[2] = inputSizeItk[2]; - - outputSpacingItk[0] = itkImage->GetSpacing()[0] * (static_cast(inputSizeItk[0]) / static_cast(outputSizeItk[0])); - outputSpacingItk[1] = itkImage->GetSpacing()[1] * (static_cast(inputSizeItk[1]) / static_cast(outputSizeItk[1])); + outputSpacingItk[0] = outputSpacing[0]; + outputSpacingItk[1] = outputSpacing[1]; outputSpacingItk[2] = itkImage->GetSpacing()[2]; - typedef itk::IdentityTransform TransformType; - T_Interpolator::Pointer _pInterpolator = T_Interpolator::New(); + outputSizeItk[0] = outputSizeItk[0] * (inputImage->GetGeometry()->GetSpacing()[0] / outputSpacing[0]); + outputSizeItk[1] = outputSizeItk[1] * (inputImage->GetGeometry()->GetSpacing()[1] / outputSpacing[1]); + resampleImageFilter->SetInput(itkImage); resampleImageFilter->SetSize(outputSizeItk); resampleImageFilter->SetOutputSpacing(outputSpacingItk); - resampleImageFilter->SetTransform(TransformType::New()); - resampleImageFilter->SetInterpolator(_pInterpolator); resampleImageFilter->UpdateLargestPossibleRegion(); return mitk::GrabItkImageMemory(resampleImageFilter->GetOutput()); } -mitk::Image::Pointer mitk::PhotoacousticFilterService::ApplyCropping(mitk::Image::Pointer inputImage, +mitk::Image::Pointer mitk::PhotoacousticFilterService::ApplyCropping( + mitk::Image::Pointer inputImage, int above, int below, int right, int left, int zStart, int zEnd) { mitk::CropImageFilter::Pointer cropImageFilter = mitk::CropImageFilter::New(); cropImageFilter->SetInput(inputImage); cropImageFilter->SetXPixelsCropStart(left); cropImageFilter->SetXPixelsCropEnd(right); cropImageFilter->SetYPixelsCropStart(above); cropImageFilter->SetYPixelsCropEnd(below); cropImageFilter->SetZPixelsCropStart(zStart); cropImageFilter->SetZPixelsCropEnd(zEnd); cropImageFilter->Update(); return cropImageFilter->GetOutput(); } -mitk::Image::Pointer mitk::PhotoacousticFilterService::ApplyBeamforming(mitk::Image::Pointer inputImage, - BeamformingSettings::Pointer config, std::function progressHandle) +mitk::Image::Pointer mitk::PhotoacousticFilterService::ApplyBeamforming( + mitk::Image::Pointer inputImage, + BeamformingSettings::Pointer config, + std::function progressHandle) { Image::Pointer processedImage = mitk::Image::New(); if (inputImage->GetDimension() != 3) { mitk::Convert2Dto3DImageFilter::Pointer dimensionImageFilter = mitk::Convert2Dto3DImageFilter::New(); dimensionImageFilter->SetInput(inputImage); dimensionImageFilter->Update(); processedImage = dimensionImageFilter->GetOutput(); } else { processedImage = inputImage; } // perform the beamforming m_BeamformingFilter = mitk::BeamformingFilter::New(config); m_BeamformingFilter->SetInput(processedImage); m_BeamformingFilter->SetProgressHandle(progressHandle); m_BeamformingFilter->UpdateLargestPossibleRegion(); processedImage = m_BeamformingFilter->GetOutput(); return processedImage; } -mitk::Image::Pointer mitk::PhotoacousticFilterService::BandpassFilter(mitk::Image::Pointer data, float recordTime, +mitk::Image::Pointer mitk::PhotoacousticFilterService::BandpassFilter( + mitk::Image::Pointer data, float recordTime, float BPHighPass, float BPLowPass, float alphaHighPass, float alphaLowPass) { bool powerOfTwo = false; int finalPower = 0; for (int i = 1; pow(2, i) <= data->GetDimension(1); ++i) { finalPower = i; if (pow(2, i) == data->GetDimension(1)) { powerOfTwo = true; } } if (!powerOfTwo) { - unsigned int dim[2] = { data->GetDimension(0), (unsigned int)pow(2,finalPower + 1) }; + double spacing_x = data->GetGeometry()->GetSpacing()[0]; + double spacing_y = data->GetGeometry()->GetSpacing()[1] * (((double)data->GetDimensions()[1]) / pow(2, finalPower + 1)); + double dim[2] = { spacing_x, spacing_y }; data = ApplyResampling(data, dim); } MITK_INFO << data->GetDimension(0); + MITK_INFO << data->GetDimension(1); // do a fourier transform, multiply with an appropriate window for the filter, and transform back - typedef float PixelType; - typedef itk::Image< PixelType, 3 > RealImageType; + typedef itk::Image< float, 3 > RealImageType; RealImageType::Pointer image; mitk::CastToItkImage(data, 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_WARN << "Bandpass could not be applied"; - return data; + MITK_ERROR << "Bandpass could not be applied"; + mitkThrow() << "bandpass error"; } float singleVoxel = 1 / (recordTime / data->GetDimension(1)) / 2 / 1000; float cutoffPixelHighPass = std::min(BPHighPass / singleVoxel, (float)data->GetDimension(1) / 2); float cutoffPixelLowPass = std::min(BPLowPass / singleVoxel, (float)data->GetDimension(1) / 2 - cutoffPixelHighPass); RealImageType::Pointer fftMultiplicator = BPFunction(data, cutoffPixelHighPass, cutoffPixelLowPass, alphaHighPass, alphaLowPass); typedef itk::MultiplyImageFilter< ComplexImageType, RealImageType, ComplexImageType > MultiplyFilterType; MultiplyFilterType::Pointer multiplyFilter = MultiplyFilterType::New(); multiplyFilter->SetInput1(forwardFFTFilter->GetOutput()); multiplyFilter->SetInput2(fftMultiplicator); - /*itk::ComplexToModulusImageFilter::Pointer toReal = itk::ComplexToModulusImageFilter::New(); - toReal->SetInput(forwardFFTFilter->GetOutput()); - return GrabItkImageMemory(toReal->GetOutput()); - return GrabItkImageMemory(fftMultiplicator); *///DEBUG - typedef itk::FFT1DComplexConjugateToRealImageFilter< ComplexImageType, RealImageType > InverseFilterType; InverseFilterType::Pointer inverseFFTFilter = InverseFilterType::New(); inverseFFTFilter->SetInput(multiplyFilter->GetOutput()); inverseFFTFilter->SetDirection(1); return GrabItkImageMemory(inverseFFTFilter->GetOutput()); } itk::Image::Pointer mitk::PhotoacousticFilterService::BPFunction(mitk::Image::Pointer reference, - int cutoffFrequencyPixelHighPass, - int cutoffFrequencyPixelLowPass, + float cutoffFrequencyPixelHighPass, + float cutoffFrequencyPixelLowPass, float alphaHighPass, float alphaLowPass) { float* imageData = new float[reference->GetDimension(0)*reference->GetDimension(1)]; - float width = reference->GetDimension(1) / 2.0 - (float)cutoffFrequencyPixelHighPass - (float)cutoffFrequencyPixelLowPass; - float center = (float)cutoffFrequencyPixelHighPass / 2.0 + width / 2.0; + float width = reference->GetDimension(1) / 2.0 - cutoffFrequencyPixelHighPass - cutoffFrequencyPixelLowPass; + float center = cutoffFrequencyPixelHighPass / 2.0 + width / 2.0; 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) { if (alphaHighPass > 0.00001) { imageData[reference->GetDimension(0)*(int)(n + center - (width / 2))] = (1 + cos(itk::Math::pi*(2 * n / (alphaHighPass*(width - 1)) - 1))) / 2; } else { imageData[reference->GetDimension(0)*(int)(n + center - (width / 2))] = 1; } } - else if (n >= (width - 1)*(1 - alphaLowPass / 2)) //??? + else if (n >= (width - 1)*(1 - alphaLowPass / 2)) { if (alphaLowPass > 0.00001) { imageData[reference->GetDimension(0)*(int)(n + center - (width / 2))] = (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; } } - //MITK_INFO << "n:" << n << " is " << imageData[reference->GetDimension(0)*(int)(n + center - (width / 2))]; } - MITK_INFO << "width: " << width << ", center: " << center << ", alphaHighPass: " << alphaHighPass << ", alphaLowPass: " << alphaLowPass; - // mirror the first half of the image 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)]; } - // copy and paste to all lines 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; } diff --git a/Modules/PhotoacousticsAlgorithms/test/mitkPAFilterServiceTest.cpp b/Modules/PhotoacousticsAlgorithms/test/mitkPAFilterServiceTest.cpp index 6bf86cfb2f..281d9ee49d 100644 --- a/Modules/PhotoacousticsAlgorithms/test/mitkPAFilterServiceTest.cpp +++ b/Modules/PhotoacousticsAlgorithms/test/mitkPAFilterServiceTest.cpp @@ -1,111 +1,111 @@ /*=================================================================== 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 // us #include #include #include #include #include #include #include class mitkPAFilterServiceTestSuite : public mitk::TestFixture { CPPUNIT_TEST_SUITE(mitkPAFilterServiceTestSuite); MITK_TEST(testRunning); CPPUNIT_TEST_SUITE_END(); private: mitk::PhotoacousticFilterService::Pointer m_PhotoacousticFilterService; public: void setUp() override { m_PhotoacousticFilterService = mitk::PhotoacousticFilterService::New(); } mitk::BeamformingSettings::Pointer CreateEmptyBeamformingSettings() { mitk::BeamformingSettings::Pointer config = mitk::BeamformingSettings::New(); return config; } mitk::BeamformingSettings::Pointer CreateFunctionalBeamformingSettings() { mitk::BeamformingSettings::Pointer config = mitk::BeamformingSettings::New(); return config; } void testRunning() { int xDim = 3; int yDim = 3; int length = yDim * xDim; float* testArray = new float[length]; for (int i = 0; i < length; ++i) { testArray[i] = 0; } unsigned int NUMBER_OF_SPATIAL_DIMENSIONS = 2; unsigned int* dimensions = new unsigned int[NUMBER_OF_SPATIAL_DIMENSIONS]; dimensions[0] = yDim; dimensions[1] = xDim; mitk::PixelType pixelType = mitk::MakeScalarPixelType(); mitk::Image::Pointer testImage = mitk::Image::New(); testImage->Initialize(pixelType, NUMBER_OF_SPATIAL_DIMENSIONS, dimensions); - testImage->SetImportSlice(testArray, 0, 0, 0, mitk::Image::ImportMemoryManagementType::RtlCopyMemory); + testImage->SetImportSlice(testArray, 0, 0, 0, mitk::Image::ImportMemoryManagementType::CopyMemory); delete[] testArray; mitk::ImageReadAccessor readAccessInput(testImage); const float* inputArray = (const float*)readAccessInput.GetData(); for (int i = 0; i < length; ++i) { CPPUNIT_ASSERT_MESSAGE(std::string("Input array already not correct: " + std::to_string(inputArray[i])), abs(inputArray[i]) < 1e-5f); } mitk::BeamformingSettings::Pointer config = CreateEmptyBeamformingSettings(); auto output = m_PhotoacousticFilterService->ApplyBeamforming(testImage, config); mitk::ImageReadAccessor readAccess(output); const float* outputArray = (const float*)readAccess.GetData(); for (int i = 0; i < length; ++i) { CPPUNIT_ASSERT_MESSAGE(std::string("Ouput array not correct: " + std::to_string(abs(outputArray[i]))), abs(outputArray[i]) < 1e-5f); } } void tearDown() override { m_PhotoacousticFilterService = nullptr; } }; MITK_TEST_SUITE_REGISTRATION(mitkPAFilterService)