diff --git a/Modules/PhotoacousticsAlgorithms/MitkPABeamformingTool/PABeamformingTool.cpp b/Modules/PhotoacousticsAlgorithms/MitkPABeamformingTool/PABeamformingTool.cpp index fd029145e8..e773b3467f 100644 --- a/Modules/PhotoacousticsAlgorithms/MitkPABeamformingTool/PABeamformingTool.cpp +++ b/Modules/PhotoacousticsAlgorithms/MitkPABeamformingTool/PABeamformingTool.cpp @@ -1,233 +1,229 @@ /*=================================================================== 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 #include struct InputParameters { mitk::Image::Pointer inputImage; std::string outputFilename; bool verbose; float speedOfSound; unsigned int cutoff; float angle; unsigned int samples; mitk::BeamformingSettings::BeamformingAlgorithm algorithm; }; InputParameters parseInput(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setCategory("MITK-Photoacoustics"); parser.setTitle("Mitk Photoacoustics Beamforming Tool"); parser.setDescription("Reads a nrrd file as an input and applies a beamforming method as set with the parameters."); parser.setContributor("Computer Assisted Medical Interventions, DKFZ"); parser.setArgumentPrefix("--", "-"); parser.beginGroup("Required parameters"); parser.addArgument( "inputImage", "i", mitkCommandLineParser::InputImage, "Input image (mitk::Image)", "input image (.nrrd file)", us::Any(), false); parser.addArgument( "output", "o", mitkCommandLineParser::OutputFile, "Output filename", "output image (.nrrd file)", us::Any(), false); parser.endGroup(); parser.beginGroup("Optional parameters"); parser.addArgument( "verbose", "v", mitkCommandLineParser::Bool, "Verbose Output", "Whether to produce verbose, or rather debug output. (default: false)"); parser.addArgument( "speed-of-sound", "sos", mitkCommandLineParser::Float, "Speed of Sound [m/s]", "The average speed of sound as assumed for the reconstruction in [m/s]. (default: 1500)"); parser.addArgument( "cutoff", "co", mitkCommandLineParser::Int, "cutoff margin on the top of the image [pixels]", "The number of pixels to be ignored for this filter in [pixels] (default: 0)."); parser.addArgument( "angle", "a", mitkCommandLineParser::Float, "field of view of the transducer elements [degrees]", "The field of view of each individual transducer element [degrees] (default: 27)."); parser.addArgument( "samples", "s", mitkCommandLineParser::Int, "samples per reconstruction line [pixels]", "The pixels along the y axis in the beamformed image [pixels] (default: 2048)."); parser.addArgument( "algorithm", "alg", mitkCommandLineParser::String, "one of [\"DAS\", \"DMAS\", \"sDMAS\"]", "The beamforming algorithm to be used for reconstruction (default: DAS)."); parser.endGroup(); InputParameters input; std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size() == 0) exit(-1); if (parsedArgs.count("verbose")) { input.verbose = true; } else { input.verbose = false; } MITK_INFO(input.verbose) << "### VERBOSE OUTPUT ENABLED ###"; if (parsedArgs.count("inputImage")) { MITK_INFO(input.verbose) << "Reading input image..."; input.inputImage = mitk::IOUtil::Load(us::any_cast(parsedArgs["inputImage"])); MITK_INFO(input.verbose) << "Reading input image...[Done]"; } else { mitkThrow() << "No input image given."; } if (parsedArgs.count("output")) { input.outputFilename = us::any_cast(parsedArgs["output"]); } else { mitkThrow() << "No output image path given.."; } if (parsedArgs.count("speed-of-sound")) { input.speedOfSound = us::any_cast(parsedArgs["speed-of-sound"]); } else { input.speedOfSound = 1500; } if (parsedArgs.count("cutoff")) { input.cutoff = us::any_cast(parsedArgs["cutoff"]); } else { input.cutoff = 0; } if (parsedArgs.count("angle")) { input.angle = us::any_cast(parsedArgs["angle"]); } else { input.angle = 27; } if (parsedArgs.count("samples")) { input.samples = us::any_cast(parsedArgs["samples"]); } else { input.samples = 2048; } if (parsedArgs.count("algorithm")) { std::string algorithm = us::any_cast(parsedArgs["algorithm"]); MITK_INFO(input.verbose) << "Parsing algorithm: " << algorithm; if (algorithm == "DAS") input.algorithm = mitk::BeamformingSettings::BeamformingAlgorithm::DAS; else if (algorithm == "DMAS") input.algorithm = mitk::BeamformingSettings::BeamformingAlgorithm::DMAS; else if (algorithm == "sDMAS") input.algorithm = mitk::BeamformingSettings::BeamformingAlgorithm::sDMAS; else { MITK_INFO(input.verbose) << "Not a valid beamforming algorithm: " << algorithm << " Reverting to DAS"; input.algorithm = mitk::BeamformingSettings::BeamformingAlgorithm::DAS; } MITK_INFO(input.verbose) << "Sucessfully set algorithm: " << algorithm; } else { input.algorithm = mitk::BeamformingSettings::BeamformingAlgorithm::DAS; MITK_INFO(input.verbose) << "No matching algorithm found. Using DAS."; } return input; } mitk::BeamformingSettings::Pointer ParseSettings(InputParameters &input) { mitk::BeamformingSettings::Pointer outputSettings = mitk::BeamformingSettings::New( (float)(input.inputImage->GetGeometry()->GetSpacing()[0] / 1000), (float)(input.speedOfSound), (float)(input.inputImage->GetGeometry()->GetSpacing()[1] / 1000000), input.angle, true, input.inputImage->GetDimension(1), input.inputImage->GetDimension(0), - input.cutoff, - false, - (unsigned int*) nullptr, input.inputImage->GetDimensions(), false, + 16, mitk::BeamformingSettings::DelayCalc::Spherical, mitk::BeamformingSettings::Apodization::Box, input.inputImage->GetDimension(0), - input.algorithm, - false, - 0.0f, - 0.0f); + input.algorithm + ); return outputSettings; } int main(int argc, char * argv[]) { auto input = parseInput(argc, argv); MITK_INFO(input.verbose) << "Beamforming input image..."; mitk::PhotoacousticFilterService::Pointer m_BeamformingService = mitk::PhotoacousticFilterService::New(); mitk::BeamformingSettings::Pointer settings = ParseSettings(input); mitk::CastToFloatImageFilter::Pointer castFilter = mitk::CastToFloatImageFilter::New(); castFilter->SetInput(input.inputImage); castFilter->Update(); auto floatImage = castFilter->GetOutput(); auto output = m_BeamformingService->ApplyBeamforming(floatImage, settings); MITK_INFO(input.verbose) << "Applying BModeFilter to image..."; - auto output2 = m_BeamformingService->ApplyBmodeFilter(output, mitk::PhotoacousticFilterService::EnvelopeDetection, false, 0.3); + auto output2 = m_BeamformingService->ApplyBmodeFilter(output, mitk::PhotoacousticFilterService::EnvelopeDetection, false); MITK_INFO(input.verbose) << "Applying BModeFilter to image...[Done]"; MITK_INFO(input.verbose) << "Saving image..."; mitk::IOUtil::Save(output2, input.outputFilename); MITK_INFO(input.verbose) << "Saving image...[Done]"; MITK_INFO(input.verbose) << "Beamforming input image...[Done]"; } diff --git a/Modules/PhotoacousticsAlgorithms/include/OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.h b/Modules/PhotoacousticsAlgorithms/include/OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.h index 508c4d05ca..7490c71f57 100644 --- a/Modules/PhotoacousticsAlgorithms/include/OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.h +++ b/Modules/PhotoacousticsAlgorithms/include/OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.h @@ -1,146 +1,145 @@ /*=================================================================== 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 _MITKPHOTOACOUSTICSOCLBEAMFORMER_H_ #define _MITKPHOTOACOUSTICSOCLBEAMFORMER_H_ #include #if defined(PHOTOACOUSTICS_USE_GPU) || DOXYGEN #include "mitkOclDataSetToDataSetFilter.h" - +#include "mitkBeamformingSettings.h" #include "mitkPhotoacousticOCLDelayCalculation.h" #include "mitkPhotoacousticOCLUsedLinesCalculation.h" -#include "mitkBeamformingSettings.h" #include namespace mitk { /*! * \brief Class implementing a mitk::OclDataSetToDataSetFilter for beamforming on GPU * * The class must be given a configuration class instance of mitk::BeamformingSettings for beamforming parameters through mitk::PhotoacousticOCLBeamformingFilter::SetConfig(BeamformingSettings settings) * Additional configuration of the apodisation function is needed. */ class PhotoacousticOCLBeamformingFilter : public OclDataSetToDataSetFilter, public itk::Object { public: mitkClassMacroItkParent(PhotoacousticOCLBeamformingFilter, itk::Object); mitkNewMacro1Param(Self, BeamformingSettings::Pointer); /** * @brief SetInput Set the input data through an image. Arbitrary images are supported */ void SetInput(Image::Pointer image, unsigned int inputSlices); /** * brief SetInput Manually set the input data while providing 3 dimensions and memory size of the input data (Bytes per element). */ void SetInput(void* data, unsigned int* dimensions, unsigned int BpE); /** * @brief GetOutput Get a pointer to the processed data. The standard datatype is float. */ void* GetOutput(); /** * @brief GetOutputAsImage Returns an mitk::Image constructed from the processed data */ mitk::Image::Pointer GetOutputAsImage(); /** \brief Update the filter */ void Update(); /** \brief Set the Apodisation function to apply when beamforming */ void SetApodisation(const float* apodisation, unsigned short apodArraySize) { m_ApodArraySize = apodArraySize; m_Apodisation = apodisation; } protected: PhotoacousticOCLBeamformingFilter(BeamformingSettings::Pointer settings); virtual ~PhotoacousticOCLBeamformingFilter(); /** \brief Initialize the filter */ bool Initialize(); /** \brief Updated the used data for beamforming depending on whether the configuration has significantly changed */ void UpdateDataBuffers(); /** \brief Execute the filter */ void Execute(); mitk::PixelType GetOutputType() { return mitk::MakeScalarPixelType(); } int GetBytesPerElem() { return sizeof(float); } virtual us::Module* GetModule(); private: /** The OpenCL kernel for the filter */ cl_kernel m_PixelCalculation; unsigned int m_OutputDim[3]; const float* m_Apodisation; unsigned short m_ApodArraySize; unsigned int m_inputSlices; unsigned short m_PAImage; BeamformingSettings::Pointer m_Conf; mitk::Image::Pointer m_InputImage; size_t m_ChunkSize[3]; mitk::OCLUsedLinesCalculation::Pointer m_UsedLinesCalculation; mitk::OCLDelayCalculation::Pointer m_DelayCalculation; cl_mem m_ApodizationBuffer; cl_mem m_MemoryLocationsBuffer; cl_mem m_DelaysBuffer; cl_mem m_UsedLinesBuffer; }; } #else namespace mitk { class PhotoacousticOCLBeamformingFilter : public itk::Object { public: mitkClassMacroItkParent(mitk::PhotoacousticOCLBeamformingFilter, itk::Object); itkNewMacro(Self); protected: /** Constructor */ PhotoacousticOCLBeamformingFilter() {} /** Destructor */ ~PhotoacousticOCLBeamformingFilter() override {} }; } #endif #endif diff --git a/Modules/PhotoacousticsAlgorithms/include/mitkCropImageFilter.h b/Modules/PhotoacousticsAlgorithms/include/mitkCropImageFilter.h index 888cc42eeb..017cc8327f 100644 --- a/Modules/PhotoacousticsAlgorithms/include/mitkCropImageFilter.h +++ b/Modules/PhotoacousticsAlgorithms/include/mitkCropImageFilter.h @@ -1,69 +1,69 @@ /*=================================================================== 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 MITK_CROP_IMAGE_FILTER #define MITK_CROP_IMAGE_FILTER #include #include "mitkImageToImageFilter.h" #include "itkMacro.h" namespace mitk { /*! * \brief Class implementing an mitk::ImageToImageFilter for casting any mitk image to a float image */ class MITKPHOTOACOUSTICSALGORITHMS_EXPORT CropImageFilter : public ImageToImageFilter { public: mitkClassMacro(CropImageFilter, ImageToImageFilter); itkFactorylessNewMacro(Self); itkCloneMacro(Self); itkGetMacro(XPixelsCropStart, unsigned int); itkSetMacro(XPixelsCropStart, unsigned int); itkGetMacro(YPixelsCropStart, unsigned int); itkSetMacro(YPixelsCropStart, unsigned int); itkGetMacro(ZPixelsCropStart, unsigned int); itkSetMacro(ZPixelsCropStart, unsigned int); itkGetMacro(XPixelsCropEnd, unsigned int); itkSetMacro(XPixelsCropEnd, unsigned int); itkGetMacro(YPixelsCropEnd, unsigned int); itkSetMacro(YPixelsCropEnd, unsigned int); itkGetMacro(ZPixelsCropEnd, unsigned int); itkSetMacro(ZPixelsCropEnd, unsigned int); protected: CropImageFilter(); ~CropImageFilter() override; void GenerateData() override; void SanityCheckPreconditions(); unsigned int m_XPixelsCropStart; - unsigned int m_XPixelsCropEnd; unsigned int m_YPixelsCropStart; - unsigned int m_YPixelsCropEnd; unsigned int m_ZPixelsCropStart; - unsigned int m_ZPixelsCropEnd; - }; + unsigned int m_XPixelsCropEnd; + unsigned int m_YPixelsCropEnd; + unsigned int m_ZPixelsCropEnd; + }; } // namespace mitk #endif //MITK_CROP_IMAGE_FILTER diff --git a/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.cpp b/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.cpp index 4d02104ace..3c3759fe24 100644 --- a/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.cpp +++ b/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.cpp @@ -1,259 +1,259 @@ /*=================================================================== 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_inputSlices(1), + m_Conf(settings), m_InputImage(mitk::Image::New()), m_ApodizationBuffer(nullptr), m_MemoryLocationsBuffer(nullptr), m_DelaysBuffer(nullptr), - m_UsedLinesBuffer(nullptr), - m_inputSlices(1), - m_Conf(settings) + m_UsedLinesBuffer(nullptr) { MITK_INFO << "Instantiating 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 << "Instantiating 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 { size_t outputSize = (size_t)m_Conf->GetReconstructionLines() * (size_t)m_Conf->GetSamplesPerLine() * (size_t)m_inputSlices; m_OutputDim[0] = m_Conf->GetReconstructionLines(); m_OutputDim[1] = m_Conf->GetSamplesPerLine(); m_OutputDim[2] = m_inputSlices; 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), &(m_OutputDim[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 2D/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() { 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); return (OclFilter::IsInitialized() && buildErr); } void mitk::PhotoacousticOCLBeamformingFilter::SetInput(mitk::Image::Pointer image, unsigned int inputSlices) { OclDataSetToDataSetFilter::SetInput(image); m_InputImage = image; m_inputSlices = inputSlices; } 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::CopyMemory); - delete[] pData; + free(pData); } MITK_DEBUG << "Image Initialized."; return outputImage; } void* mitk::PhotoacousticOCLBeamformingFilter::GetOutput() { return OclDataSetToDataSetFilter::GetOutput(); } #endif diff --git a/Modules/PhotoacousticsAlgorithms/source/filters/mitkBandpassFilter.cpp b/Modules/PhotoacousticsAlgorithms/source/filters/mitkBandpassFilter.cpp index fe46b66d99..901622cbae 100644 --- a/Modules/PhotoacousticsAlgorithms/source/filters/mitkBandpassFilter.cpp +++ b/Modules/PhotoacousticsAlgorithms/source/filters/mitkBandpassFilter.cpp @@ -1,274 +1,274 @@ /*=================================================================== 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." +#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."; } void mitk::BandpassFilter::SanityCheckPreconditions() { auto input = GetInput(); std::string type = input->GetPixelType().GetTypeAsString(); if (!(type == "scalar (float)" || type == " (float)")) { MITK_ERROR << "This filter can currently only handle float type images."; mitkThrow() << "This filter can currently only handle float type images."; } if(m_HighPass > m_LowPass) { MITK_ERROR << "High Pass is higher than Low Pass; aborting."; mitkThrow() << "High Pass is higher than Low Pass; aborting."; } } 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 + 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)*(int)(n + center - (width / 2.f))] = 1; if (n <= (alphaHighPass*(width - 1)) / 2.f) { if (alphaHighPass > 0.00001f) { 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.f))] = 1; } } else if (n >= (width - 1)*(1 - alphaLowPass / 2.f)) { if (alphaLowPass > 0.00001f) { 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.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() { SanityCheckPreconditions(); 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 = 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; MITK_INFO << "cutoffPixelHighPass: " << cutoffPixelHighPass; MITK_INFO << "cutoffPixelLowPass: " << cutoffPixelLowPass; RealImageType::Pointer fftMultiplicator = BPFunction(resampledInput, cutoffPixelHighPass, cutoffPixelLowPass, m_HighPassAlpha, m_LowPassAlpha); 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/source/filters/mitkBeamformingFilter.cpp b/Modules/PhotoacousticsAlgorithms/source/filters/mitkBeamformingFilter.cpp index ae0b924a76..9c0b2eafb7 100644 --- a/Modules/PhotoacousticsAlgorithms/source/filters/mitkBeamformingFilter.cpp +++ b/Modules/PhotoacousticsAlgorithms/source/filters/mitkBeamformingFilter.cpp @@ -1,294 +1,297 @@ /*=================================================================== mitkBeamformingFilter 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 "mitkProperties.h" #include "mitkImageReadAccessor.h" #include #include #include #include #include #include "mitkImageCast.h" #include "mitkBeamformingFilter.h" #include "mitkBeamformingUtils.h" mitk::BeamformingFilter::BeamformingFilter(mitk::BeamformingSettings::Pointer settings) : m_OutputData(nullptr), m_InputData(nullptr), m_Conf(settings) { MITK_INFO << "Instantiating BeamformingFilter..."; this->SetNumberOfIndexedInputs(1); this->SetNumberOfRequiredInputs(1); m_ProgressHandle = [](int, std::string) {}; - +#if defined(PHOTOACOUSTICS_USE_GPU) m_BeamformingOclFilter = mitk::PhotoacousticOCLBeamformingFilter::New(m_Conf); +#else + m_BeamformingOclFilter = mitk::PhotoacousticOCLBeamformingFilter::New(); +#endif MITK_INFO << "Instantiating BeamformingFilter...[Done]"; } void mitk::BeamformingFilter::SetProgressHandle(std::function progressHandle) { m_ProgressHandle = progressHandle; } mitk::BeamformingFilter::~BeamformingFilter() { MITK_INFO << "Destructed BeamformingFilter"; } void mitk::BeamformingFilter::GenerateInputRequestedRegion() { Superclass::GenerateInputRequestedRegion(); mitk::Image* output = this->GetOutput(); mitk::Image* input = const_cast (this->GetInput()); if (!output->IsInitialized()) { return; } input->SetRequestedRegionToLargestPossibleRegion(); } void mitk::BeamformingFilter::GenerateOutputInformation() { mitk::Image::ConstPointer input = this->GetInput(); mitk::Image::Pointer output = this->GetOutput(); if ((output->IsInitialized()) && (this->GetMTime() <= m_TimeOfHeaderInitialization.GetMTime())) return; unsigned int dim[] = { m_Conf->GetReconstructionLines(), m_Conf->GetSamplesPerLine(), input->GetDimension(2) }; output->Initialize(mitk::MakeScalarPixelType(), 3, dim); mitk::Vector3D spacing; spacing[0] = m_Conf->GetPitchInMeters() * m_Conf->GetTransducerElements() * 1000 / m_Conf->GetReconstructionLines(); spacing[1] = (m_Conf->GetTimeSpacing() * m_Conf->GetInputDim()[1]) / (2-(int)m_Conf->GetIsPhotoacousticImage()) * m_Conf->GetSpeedOfSound() * 1000 / m_Conf->GetSamplesPerLine(); spacing[2] = 1; output->GetGeometry()->SetSpacing(spacing); output->GetGeometry()->Modified(); output->SetPropertyList(input->GetPropertyList()->Clone()); m_TimeOfHeaderInitialization.Modified(); } void mitk::BeamformingFilter::GenerateData() { mitk::Image::Pointer input = this->GetInput(); if (!(input->GetPixelType().GetTypeAsString() == "scalar (float)" || input->GetPixelType().GetTypeAsString() == " (float)")) { MITK_ERROR << "Pixel type of input needs to be float for this filter to work."; mitkThrow() << "Pixel type of input needs to be float for this filter to work."; } GenerateOutputInformation(); mitk::Image::Pointer output = this->GetOutput(); if (!output->IsInitialized()) return; auto begin = std::chrono::high_resolution_clock::now(); // debbuging the performance... if (!m_Conf->GetUseGPU()) { int progInterval = output->GetDimension(2) / 20 > 1 ? output->GetDimension(2) / 20 : 1; // the interval at which we update the gui progress bar float inputDim[2] = { (float)input->GetDimension(0), (float)input->GetDimension(1) }; float outputDim[2] = { (float)output->GetDimension(0), (float)output->GetDimension(1) }; for (unsigned int i = 0; i < output->GetDimension(2); ++i) // seperate Slices should get Beamforming seperately applied { mitk::ImageReadAccessor inputReadAccessor(input, input->GetSliceData(i)); m_InputData = (float*)inputReadAccessor.GetData(); m_OutputData = new float[m_Conf->GetReconstructionLines()*m_Conf->GetSamplesPerLine()]; // fill the image with zeros for (int l = 0; l < outputDim[0]; ++l) { for (int s = 0; s < outputDim[1]; ++s) { m_OutputData[l*(short)outputDim[1] + s] = 0; } } std::thread *threads = new std::thread[(short)outputDim[0]]; // every line will be beamformed in a seperate thread if (m_Conf->GetAlgorithm() == BeamformingSettings::BeamformingAlgorithm::DAS) { if (m_Conf->GetDelayCalculationMethod() == BeamformingSettings::DelayCalc::QuadApprox) { for (short line = 0; line < outputDim[0]; ++line) { threads[line] = std::thread(&BeamformingUtils::DASQuadraticLine, m_InputData, m_OutputData, inputDim, outputDim, line, m_Conf); } } else if (m_Conf->GetDelayCalculationMethod() == BeamformingSettings::DelayCalc::Spherical) { for (short line = 0; line < outputDim[0]; ++line) { threads[line] = std::thread(&BeamformingUtils::DASSphericalLine, m_InputData, m_OutputData, inputDim, outputDim, line, m_Conf); } } } else if (m_Conf->GetAlgorithm() == BeamformingSettings::BeamformingAlgorithm::DMAS) { if (m_Conf->GetDelayCalculationMethod() == BeamformingSettings::DelayCalc::QuadApprox) { for (short line = 0; line < outputDim[0]; ++line) { threads[line] = std::thread(&BeamformingUtils::DMASQuadraticLine, m_InputData, m_OutputData, inputDim, outputDim, line, m_Conf); } } else if (m_Conf->GetDelayCalculationMethod() == BeamformingSettings::DelayCalc::Spherical) { for (short line = 0; line < outputDim[0]; ++line) { threads[line] = std::thread(&BeamformingUtils::DMASSphericalLine, m_InputData, m_OutputData, inputDim, outputDim, line, m_Conf); } } } else if (m_Conf->GetAlgorithm() == BeamformingSettings::BeamformingAlgorithm::sDMAS) { if (m_Conf->GetDelayCalculationMethod() == BeamformingSettings::DelayCalc::QuadApprox) { for (short line = 0; line < outputDim[0]; ++line) { threads[line] = std::thread(&BeamformingUtils::sDMASQuadraticLine, m_InputData, m_OutputData, inputDim, outputDim, line, m_Conf); } } else if (m_Conf->GetDelayCalculationMethod() == BeamformingSettings::DelayCalc::Spherical) { for (short line = 0; line < outputDim[0]; ++line) { threads[line] = std::thread(&BeamformingUtils::sDMASSphericalLine, m_InputData, m_OutputData, inputDim, outputDim, line, m_Conf); } } } // wait for all lines to finish for (short line = 0; line < outputDim[0]; ++line) { threads[line].join(); } output->SetSlice(m_OutputData, i); if (i % progInterval == 0) m_ProgressHandle((int)((i + 1) / (float)output->GetDimension(2) * 100), "performing reconstruction"); delete[] m_OutputData; m_OutputData = nullptr; m_InputData = nullptr; } } #if defined(PHOTOACOUSTICS_USE_GPU) || DOXYGEN else { try { // first, we check whether the data is float, other formats are unsupported 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"; } unsigned long availableMemory = m_BeamformingOclFilter->GetDeviceMemory(); unsigned int batchSize = m_Conf->GetGPUBatchSize(); unsigned int batches = (unsigned int)((float)input->GetDimension(2) / batchSize) + (input->GetDimension(2) % batchSize > 0); unsigned int batchDim[] = { input->GetDimension(0), input->GetDimension(1), batchSize }; unsigned int batchDimLast[] = { input->GetDimension(0), input->GetDimension(1), input->GetDimension(2) % batchSize }; // the following safeguard is probably only needed for absurdly small GPU memory if((unsigned long)batchSize * ((unsigned long)(batchDim[0] * batchDim[1]) * 4 + // single input image (float) (unsigned long)(m_Conf->GetReconstructionLines() * m_Conf->GetSamplesPerLine()) * 4) // single output image (float) > availableMemory - (unsigned long)(m_Conf->GetReconstructionLines() / 2 * m_Conf->GetSamplesPerLine()) * 2 - // Delays buffer (unsigned short) (unsigned long)(m_Conf->GetReconstructionLines() * m_Conf->GetSamplesPerLine()) * 3 * 2 - // UsedLines buffer (unsigned short) 50 * 1024 * 1024)// 50 MB buffer for local data, system purposes etc { MITK_ERROR << "device memory too small for GPU beamforming; try decreasing the batch size"; return; } mitk::ImageReadAccessor copy(input); for (unsigned int i = 0; i < batches; ++i) { m_ProgressHandle(100.f * (float)i / (float)batches, "performing reconstruction"); mitk::Image::Pointer inputBatch = mitk::Image::New(); unsigned int num_Slices = 1; if (i == batches - 1 && (input->GetDimension(2) % batchSize > 0)) { inputBatch->Initialize(mitk::MakeScalarPixelType(), 3, batchDimLast); num_Slices = batchDimLast[2]; } else { inputBatch->Initialize(mitk::MakeScalarPixelType(), 3, batchDim); num_Slices = batchDim[2]; } inputBatch->SetSpacing(input->GetGeometry()->GetSpacing()); inputBatch->SetImportVolume(&(((float*)copy.GetData())[input->GetDimension(0) * input->GetDimension(1) * batchSize * i])); m_BeamformingOclFilter->SetApodisation(m_Conf->GetApodizationFunction(), m_Conf->GetApodizationArraySize()); m_BeamformingOclFilter->SetInput(inputBatch, num_Slices); m_BeamformingOclFilter->Update(); void* out = m_BeamformingOclFilter->GetOutput(); for (unsigned int slice = 0; slice < num_Slices; ++slice) { output->SetImportSlice( &(((float*)out)[m_Conf->GetReconstructionLines() * m_Conf->GetSamplesPerLine() * slice]), batchSize * i + slice, 0, 0, mitk::Image::ImportMemoryManagementType::CopyMemory); } } } catch (mitk::Exception &e) { std::string errorMessage = "Caught unexpected exception "; errorMessage.append(e.what()); MITK_ERROR << errorMessage; float* dummyData = new float[m_Conf->GetReconstructionLines() * m_Conf->GetSamplesPerLine() * m_Conf->GetInputDim()[2]]; output->SetImportVolume(dummyData, 0, 0, mitk::Image::ImportMemoryManagementType::CopyMemory); } } #endif m_TimeOfHeaderInitialization.Modified(); auto end = std::chrono::high_resolution_clock::now(); MITK_INFO << "Beamforming of " << output->GetDimension(2) << " Images completed in " << ((float)std::chrono::duration_cast(end - begin).count()) / 1000000 << "ms" << std::endl; } diff --git a/Modules/PhotoacousticsAlgorithms/source/filters/mitkBeamformingSettings.cpp b/Modules/PhotoacousticsAlgorithms/source/filters/mitkBeamformingSettings.cpp index d956ecfdcd..b030b2c919 100644 --- a/Modules/PhotoacousticsAlgorithms/source/filters/mitkBeamformingSettings.cpp +++ b/Modules/PhotoacousticsAlgorithms/source/filters/mitkBeamformingSettings.cpp @@ -1,100 +1,103 @@ /*=================================================================== mitkBeamformingSettings 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 "mitkBeamformingSettings.h" #include "mitkBeamformingUtils.h" #include "itkMutexLock.h" mitk::BeamformingSettings::BeamformingSettings(std::string xmlFile) { - MITK_ERROR << "Not implemented yet."; - mitkThrow() << "Not implemented yet."; + if((xmlFile.length())>0) + { + MITK_ERROR << "Not implemented yet."; + mitkThrow() << "Not implemented yet."; + } } mitk::BeamformingSettings::BeamformingSettings(float pitchInMeters, float speedOfSound, float timeSpacing, float angle, bool isPhotoacousticImage, unsigned int samplesPerLine, unsigned int reconstructionLines, unsigned int* inputDim, bool useGPU, unsigned int GPUBatchSize, DelayCalc delayCalculationMethod, Apodization apod, unsigned int apodizationArraySize, BeamformingAlgorithm algorithm ) : m_PitchInMeters(pitchInMeters), m_SpeedOfSound(speedOfSound), m_TimeSpacing(timeSpacing), m_Angle(angle), m_IsPhotoacousticImage(isPhotoacousticImage), m_SamplesPerLine(samplesPerLine), m_ReconstructionLines(reconstructionLines), m_UseGPU(useGPU), m_GPUBatchSize(GPUBatchSize), m_DelayCalculationMethod(delayCalculationMethod), m_Apod(apod), m_ApodizationArraySize(apodizationArraySize), m_Algorithm(algorithm) { if (inputDim == nullptr) { MITK_ERROR << "No input dimension given."; mitkThrow() << "No input dimension given."; } switch (GetApod()) { case BeamformingSettings::Apodization::Hann: m_ApodizationFunction = mitk::BeamformingUtils::VonHannFunction(GetApodizationArraySize()); break; case BeamformingSettings::Apodization::Hamm: m_ApodizationFunction = mitk::BeamformingUtils::HammFunction(GetApodizationArraySize()); break; case BeamformingSettings::Apodization::Box: default: m_ApodizationFunction = mitk::BeamformingUtils::BoxFunction(GetApodizationArraySize()); break; } m_InputDim = new unsigned int[3]{ inputDim[0], inputDim[1], inputDim[2] }; m_TransducerElements = m_InputDim[0]; } mitk::BeamformingSettings::~BeamformingSettings() { MITK_INFO << "Destructing beamforming settings..."; //Free memory if (m_ApodizationFunction != nullptr) { MITK_INFO << "Deleting apodization function..."; delete[] m_ApodizationFunction; MITK_INFO << "Deleting apodization function...[Done]"; } if (m_InputDim != nullptr) { MITK_INFO << "Deleting input dim..."; delete[] m_InputDim; MITK_INFO << "Deleting input dim...[Done]"; } MITK_INFO << "Destructing beamforming settings...[Done]"; } diff --git a/Modules/PhotoacousticsAlgorithms/source/filters/mitkCropImageFilter.cpp b/Modules/PhotoacousticsAlgorithms/source/filters/mitkCropImageFilter.cpp index 396a2c98c4..cd35ae5a10 100644 --- a/Modules/PhotoacousticsAlgorithms/source/filters/mitkCropImageFilter.cpp +++ b/Modules/PhotoacousticsAlgorithms/source/filters/mitkCropImageFilter.cpp @@ -1,120 +1,120 @@ /*=================================================================== mitkCropImageFilter 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 "mitkCropImageFilter.h" #include #include mitk::CropImageFilter::CropImageFilter() : - m_XPixelsCropEnd(0), - m_YPixelsCropEnd(0), - m_ZPixelsCropEnd(0), m_XPixelsCropStart(0), m_YPixelsCropStart(0), - m_ZPixelsCropStart(0) + m_ZPixelsCropStart(0), + m_XPixelsCropEnd(0), + m_YPixelsCropEnd(0), + m_ZPixelsCropEnd(0) { MITK_INFO << "Instantiating CropImageFilter..."; SetNumberOfIndexedInputs(1); SetNumberOfIndexedOutputs(1); MITK_INFO << "Instantiating CropImageFilter...[Done]"; } mitk::CropImageFilter::~CropImageFilter() { MITK_INFO << "Destructed CastToFloatImageFilter."; } void mitk::CropImageFilter::SanityCheckPreconditions() { mitk::Image::Pointer inputImage = GetInput(); std::string type = inputImage->GetPixelType().GetTypeAsString(); if (!(type == "scalar (float)" || type == " (float)")) { MITK_ERROR << "This filter can currently only handle float type images."; mitkThrow() << "This filter can currently only handle float type images."; } if (m_XPixelsCropStart + m_XPixelsCropEnd >= inputImage->GetDimension(0)) { MITK_ERROR << "X Crop area too large for selected input image"; mitkThrow() << "X Crop area too large for selected input image"; } if (m_YPixelsCropStart + m_YPixelsCropEnd >= inputImage->GetDimension(1)) { MITK_ERROR << "Y Crop area too large for selected input image"; mitkThrow() << "Y Crop area too large for selected input image"; } if (inputImage->GetDimension() == 3) { if (m_ZPixelsCropStart + m_ZPixelsCropEnd >= inputImage->GetDimension(2)) { MITK_ERROR << "Y Crop area too large for selected input image"; mitkThrow() << "Z Crop area too large for selected input image"; } } else { if (m_ZPixelsCropStart + m_ZPixelsCropEnd > 0) { mitkThrow() << "Z Crop area too large for selected input image"; } } } void mitk::CropImageFilter::GenerateData() { mitk::Image::Pointer inputImage = GetInput(); mitk::Image::Pointer outputImage = GetOutput(); SanityCheckPreconditions(); unsigned int* outputDim; unsigned int* inputDim = inputImage->GetDimensions(); if (inputImage->GetDimension() == 2) outputDim = new unsigned int[2]{ inputImage->GetDimension(0) - m_XPixelsCropStart - m_XPixelsCropEnd, inputImage->GetDimension(1) - m_YPixelsCropStart - m_YPixelsCropEnd }; else outputDim = new unsigned int[3]{ inputImage->GetDimension(0) - m_XPixelsCropStart - m_XPixelsCropEnd, inputImage->GetDimension(1) - m_YPixelsCropStart - m_YPixelsCropEnd, inputImage->GetDimension(2) - m_ZPixelsCropStart - m_ZPixelsCropEnd }; outputImage->Initialize(mitk::MakeScalarPixelType(), inputImage->GetDimension(), outputDim); outputImage->SetSpacing(inputImage->GetGeometry()->GetSpacing()); ImageReadAccessor accR(inputImage); const float* inputData = (const float*)(accR.GetData()); ImageWriteAccessor accW(outputImage); float* outputData = (float*)(accW.GetData()); unsigned int zLength = inputImage->GetDimension() == 3 ? outputDim[2] : 0; for (unsigned int sl = 0; sl < zLength; ++sl) { for (unsigned int l = 0; l < outputDim[0]; ++l) { for (unsigned int s = 0; s < outputDim[1]; ++s) { outputData[l + s*outputDim[0] + sl*outputDim[0] * outputDim[1]] = inputData[(l + m_XPixelsCropStart) + (s + m_YPixelsCropStart)*inputDim[0] + (sl + m_ZPixelsCropStart)*inputDim[0] * inputDim[1]]; } } } delete[] outputDim; } diff --git a/Modules/PhotoacousticsAlgorithms/test/mitkBandpassFilterTest.cpp b/Modules/PhotoacousticsAlgorithms/test/mitkBandpassFilterTest.cpp index 06b1246bae..b3b319f35b 100644 --- a/Modules/PhotoacousticsAlgorithms/test/mitkBandpassFilterTest.cpp +++ b/Modules/PhotoacousticsAlgorithms/test/mitkBandpassFilterTest.cpp @@ -1,194 +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." +#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; + // 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; + // 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()) < 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; + // 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()) < 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 54e4504acf..310770220b 100644 --- a/Modules/PhotoacousticsAlgorithms/test/mitkBeamformingFilterTest.cpp +++ b/Modules/PhotoacousticsAlgorithms/test/mitkBeamformingFilterTest.cpp @@ -1,279 +1,279 @@ /*=================================================================== 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]; for (size_t i = 0; i < m_Samples * m_TransducerElements; ++i) { m_Data[i] = 0; } } ~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.0001f, 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; int pixels = std::round(delay_in_seconds / m_Spacing_y); for (int index = -4; index < 9; ++index) { if ((int)pixels + index < (int)m_Samples && (int)pixels + index > 0) { m_Data[(size_t)(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_DAS); MITK_TEST(testBeamformingGPU_DAS); MITK_TEST(testBeamformingCPU_DMAS); MITK_TEST(testBeamformingGPU_DMAS); MITK_TEST(testBeamformingCPU_sDMAS); MITK_TEST(testBeamformingGPU_sDMAS); CPPUNIT_TEST_SUITE_END(); private: mitk::BeamformingFilter::Pointer m_BeamformingFilter; const unsigned int NUM_ITERATIONS = 15; const unsigned int SAMPLES = 5000 * 2; 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 / 2; // us const unsigned int GPU_BATCH_SIZE = 16; public: void setUp() override { } void test() { 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 iteration = 0; iteration < NUM_ITERATIONS; ++iteration) { // create some synthetic input data float depth_in_meters1 = randDistrDepth(randGen); float depth_in_meters2 = randDistrDepth(randGen); float depth_in_meters3 = randDistrDepth(randGen); unsigned int element1 = 29; unsigned int element2 = 63; unsigned int element3 = 98; SyntheticPAImageData image(SPACING_X / 1000.f, SPACING_Y / 1000000.f, SAMPLES, ELEMENTS, SPEED_OF_SOUND); image.AddWave(depth_in_meters1, 29); image.AddWave(depth_in_meters2, 63); image.AddWave(depth_in_meters3, 98); 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((const void*)image.GetData(), mitk::Image::RtlCopyMemory); + inputImage->SetImportVolume((const void*)image.GetData(), mitk::Image::CopyMemory); // 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(); unsigned int pos1[3] = { element1, (unsigned int)std::round(depth_in_meters1 * 1000.f / outputImage->GetGeometry()->GetSpacing()[1]), 0 }; unsigned int pos2[3] = { element2, (unsigned int)std::round(depth_in_meters2 * 1000.f / outputImage->GetGeometry()->GetSpacing()[1]), 0 }; unsigned int pos3[3] = { element3, (unsigned int)std::round(depth_in_meters3 * 1000.f / outputImage->GetGeometry()->GetSpacing()[1]), 0 }; double average = 0; for (unsigned int i = 0; i < RECONSTRUCTED_LINES*RECONSTRUCTED_SAMPLES; ++i) { average += outputData[i] / (RECONSTRUCTED_LINES*RECONSTRUCTED_SAMPLES); } CPPUNIT_ASSERT_MESSAGE(std::string("Iteration " + std::to_string(iteration) + ": first point source incorrectly reconstructed; should be > average*100, is " + std::to_string(abs(outputData[pos1[0] + pos1[1] * RECONSTRUCTED_LINES]))) + " < " + std::to_string(average) + "*100" , abs(outputData[pos1[0] + pos1[1]*RECONSTRUCTED_LINES] / average) > 100); CPPUNIT_ASSERT_MESSAGE(std::string("Iteration " + std::to_string(iteration) + ": second point source incorrectly reconstructed; should be > average*100, is " + std::to_string(abs(outputData[pos2[0] + pos2[1] * RECONSTRUCTED_LINES]))) + " < " + std::to_string(average) + "*100" , abs(outputData[pos2[0] + pos2[1] * RECONSTRUCTED_LINES] / average) > 100); CPPUNIT_ASSERT_MESSAGE(std::string("Iteration " + std::to_string(iteration) + ": third point source incorrectly reconstructed; should be > average*100, is " + std::to_string(abs(outputData[pos3[0] + pos3[1] * RECONSTRUCTED_LINES]))) + " < " + std::to_string(average) + "*100" , abs(outputData[pos3[0] + pos3[1] * RECONSTRUCTED_LINES] / average) > 100); } } mitk::BeamformingSettings::Pointer createConfig(bool UseGPU, unsigned int* inputDim, mitk::BeamformingSettings::BeamformingAlgorithm alg) { return mitk::BeamformingSettings::New(SPACING_X / 1000, SPEED_OF_SOUND, SPACING_Y / 1000000, 27.f, true, RECONSTRUCTED_SAMPLES, RECONSTRUCTED_LINES, inputDim, UseGPU, GPU_BATCH_SIZE, mitk::BeamformingSettings::DelayCalc::Spherical, mitk::BeamformingSettings::Apodization::Box, ELEMENTS * 2, alg); } void testBeamformingCPU_DAS() { MITK_INFO << "Started DAS test on CPU"; unsigned int* inputDim = new unsigned int[3]; inputDim[0] = ELEMENTS; inputDim[1] = SAMPLES; inputDim[2] = 1; m_BeamformingFilter = mitk::BeamformingFilter::New(createConfig(false, inputDim, mitk::BeamformingSettings::BeamformingAlgorithm::DAS)); test(); delete[] inputDim; } void testBeamformingGPU_DAS() { MITK_INFO << "Started DAS test on GPU"; unsigned int* inputDim = new unsigned int[3]; inputDim[0] = ELEMENTS; inputDim[1] = SAMPLES; inputDim[2] = 1; m_BeamformingFilter = mitk::BeamformingFilter::New(createConfig(true, inputDim, mitk::BeamformingSettings::BeamformingAlgorithm::DAS)); test(); delete[] inputDim; } void testBeamformingCPU_sDMAS() { MITK_INFO << "Started sDMAS test on CPU"; unsigned int* inputDim = new unsigned int[3]; inputDim[0] = ELEMENTS; inputDim[1] = SAMPLES; inputDim[2] = 1; m_BeamformingFilter = mitk::BeamformingFilter::New(createConfig(false, inputDim, mitk::BeamformingSettings::BeamformingAlgorithm::sDMAS)); test(); delete[] inputDim; } void testBeamformingGPU_sDMAS() { MITK_INFO << "Started sDMAS test on GPU"; unsigned int* inputDim = new unsigned int[3]; inputDim[0] = ELEMENTS; inputDim[1] = SAMPLES; inputDim[2] = 1; m_BeamformingFilter = mitk::BeamformingFilter::New(createConfig(true, inputDim, mitk::BeamformingSettings::BeamformingAlgorithm::sDMAS)); test(); delete[] inputDim; } void testBeamformingCPU_DMAS() { MITK_INFO << "Started DMAS test on CPU"; unsigned int* inputDim = new unsigned int[3]; inputDim[0] = ELEMENTS; inputDim[1] = SAMPLES; inputDim[2] = 1; m_BeamformingFilter = mitk::BeamformingFilter::New(createConfig(false, inputDim, mitk::BeamformingSettings::BeamformingAlgorithm::DMAS)); test(); delete[] inputDim; } void testBeamformingGPU_DMAS() { MITK_INFO << "Started DMAS test on GPU"; unsigned int* inputDim = new unsigned int[3]; inputDim[0] = ELEMENTS; inputDim[1] = SAMPLES; inputDim[2] = 1; m_BeamformingFilter = mitk::BeamformingFilter::New(createConfig(true, inputDim, mitk::BeamformingSettings::BeamformingAlgorithm::DMAS)); test(); delete[] inputDim; } }; MITK_TEST_SUITE_REGISTRATION(mitkBeamformingFilter)