diff --git a/Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticOCLBeamformer.cpp b/Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticOCLBeamformer.cpp index b76e8ec278..3ec346268d 100644 --- a/Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticOCLBeamformer.cpp +++ b/Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticOCLBeamformer.cpp @@ -1,222 +1,260 @@ /*=================================================================== 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. ===================================================================*/ #ifdef PHOTOACOUSTICS_USE_GPU #include "mitkPhotoacousticOCLBeamformer.h" + #include "usServiceReference.h" mitk::PhotoacousticOCLBeamformingFilter::PhotoacousticOCLBeamformingFilter() : m_PixelCalculation( NULL ), m_InputImage(mitk::Image::New()) { this->AddSourceFile("DASQuadratic.cl"); this->AddSourceFile("DMASQuadratic.cl"); - this->AddSourceFile("DASspherical.cl"); - this->AddSourceFile("DMASspherical.cl"); + this->AddSourceFile("DMAS.cl"); unsigned int dim[] = { 128, 2048, 2 }; mitk::Vector3D spacing; spacing[0] = 1; spacing[1] = 1; spacing[2] = 1; m_InputImage->Initialize(mitk::MakeScalarPixelType(), 3, dim); m_InputImage->SetSpacing(spacing); - this->m_FilterID = "PixelCalculation"; + m_ChunkSize[0] = 128; + m_ChunkSize[1] = 128; + m_ChunkSize[2] = 8; + + m_SumFilter = mitk::OCLMemoryLocSum::New(); + m_UsedLinesCalculation = mitk::OCLUsedLinesCalculation::New(); + m_DelayCalculation = mitk::OCLDelayCalculation::New(); + + this->m_FilterID = "OpenCLBeamformingFilter"; } mitk::PhotoacousticOCLBeamformingFilter::~PhotoacousticOCLBeamformingFilter() { if ( this->m_PixelCalculation ) { clReleaseKernel( m_PixelCalculation ); } } 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(); } } +#include + void mitk::PhotoacousticOCLBeamformingFilter::Execute() { cl_int clErr = 0; try { size_t outputSize = m_OutputDim[0] * m_OutputDim[1] * m_OutputDim[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; } if (m_Apodisation == nullptr) { MITK_INFO << "No apodisation function set; Beamforming will be done without any apodisation."; m_Apodisation = new float[1]; m_Apodisation[0] = 1; m_ApodArraySize = 1; } us::ServiceReference ref = GetModuleContext()->GetServiceReference(); OclResourceService* resources = GetModuleContext()->GetService(ref); cl_context gpuContext = resources->GetContext(); - cl_mem cl_input = clCreateBuffer(gpuContext, CL_MEM_READ_ONLY | CL_MEM_USE_HOST_PTR, sizeof(float) * m_ApodArraySize, m_Apodisation, &clErr); + cl_mem ApodizationBuffer = clCreateBuffer(gpuContext, CL_MEM_READ_ONLY | CL_MEM_USE_HOST_PTR, sizeof(float) * m_ApodArraySize, m_Apodisation, &clErr); CHECK_OCL_ERR(clErr); m_PAImage = (unsigned short)m_Conf.Photoacoustic; - if (true || m_Conf.Algorithm == BeamformingFilter::beamformingSettings::BeamformingAlgorithm::DAS) + if (m_Conf.Algorithm == BeamformingFilter::beamformingSettings::BeamformingAlgorithm::DAS) { // set kernel arguments - clErr = clSetKernelArg(this->m_PixelCalculation, 2, sizeof(cl_mem), &cl_input); + clErr = clSetKernelArg(this->m_PixelCalculation, 2, sizeof(cl_mem), &ApodizationBuffer); clErr |= clSetKernelArg(this->m_PixelCalculation, 3, sizeof(cl_ushort), &(this->m_ApodArraySize)); clErr |= clSetKernelArg(this->m_PixelCalculation, 4, sizeof(cl_float), &(this->m_Conf.SpeedOfSound)); clErr |= clSetKernelArg(this->m_PixelCalculation, 5, sizeof(cl_float), &(this->m_Conf.TimeSpacing)); clErr |= clSetKernelArg(this->m_PixelCalculation, 6, sizeof(cl_float), &(this->m_Conf.Pitch)); clErr |= clSetKernelArg(this->m_PixelCalculation, 7, sizeof(cl_float), &(this->m_Conf.Angle)); clErr |= clSetKernelArg(this->m_PixelCalculation, 8, sizeof(cl_ushort), &(this->m_PAImage)); clErr |= clSetKernelArg(this->m_PixelCalculation, 9, sizeof(cl_ushort), &(this->m_Conf.TransducerElements)); clErr |= clSetKernelArg(this->m_PixelCalculation, 10, sizeof(cl_uint), &(this->m_InputDim[0])); clErr |= clSetKernelArg(this->m_PixelCalculation, 11, sizeof(cl_uint), &(this->m_InputDim[1])); clErr |= clSetKernelArg(this->m_PixelCalculation, 12, sizeof(cl_uint), &(this->m_InputDim[2])); clErr |= clSetKernelArg(this->m_PixelCalculation, 13, sizeof(cl_uint), &(this->m_OutputDim[0])); clErr |= clSetKernelArg(this->m_PixelCalculation, 14, sizeof(cl_uint), &(this->m_OutputDim[1])); CHECK_OCL_ERR(clErr); - - size_t chunkSize[3] = { 128, 128, 8 }; - - // execute the filter on a 3D NDRange - if (m_OutputDim[2] == 1) - this->ExecuteKernelChunks(m_PixelCalculation, 2, chunkSize); - else - this->ExecuteKernelChunks(m_PixelCalculation, 3, chunkSize); } - else + else if (m_Conf.Algorithm == BeamformingFilter::beamformingSettings::BeamformingAlgorithm::DMAS) { + m_UsedLinesCalculation->SetConfig(m_Conf); + m_UsedLinesCalculation->Update(); + + unsigned int sumDimensions[2] = { m_Conf.ReconstructionLines, m_Conf.SamplesPerLine }; + + m_SumFilter->SetInput(m_UsedLinesCalculation->GetGPUOutput()); + m_SumFilter->SetInputDimensions(sumDimensions); + m_SumFilter->Update(); + unsigned int size = m_SumFilter->GetSum(); + + cl_mem UsedLinesBuffer = m_UsedLinesCalculation->GetGPUOutput()->GetGPUBuffer(); + cl_mem MemoryLocationsBuffer = m_SumFilter->GetGPUOutput()->GetGPUBuffer(); + + m_DelayCalculation->SetConfig(m_Conf); + m_DelayCalculation->SetInputs(UsedLinesBuffer, MemoryLocationsBuffer, size); + m_DelayCalculation->Update(); + + cl_mem DelaysBuffer = m_DelayCalculation->GetGPUOutput()->GetGPUBuffer(); + + clErr = clSetKernelArg(this->m_PixelCalculation, 2, sizeof(cl_mem), &UsedLinesBuffer); + clErr |= clSetKernelArg(this->m_PixelCalculation, 3, sizeof(cl_mem), &MemoryLocationsBuffer); + clErr |= clSetKernelArg(this->m_PixelCalculation, 4, sizeof(cl_mem), &DelaysBuffer); + clErr |= clSetKernelArg(this->m_PixelCalculation, 5, sizeof(cl_mem), &ApodizationBuffer); + clErr |= clSetKernelArg(this->m_PixelCalculation, 6, sizeof(cl_ushort), &(this->m_ApodArraySize)); + clErr |= clSetKernelArg(this->m_PixelCalculation, 7, sizeof(cl_uint), &(this->m_InputDim[0])); + clErr |= clSetKernelArg(this->m_PixelCalculation, 8, sizeof(cl_uint), &(this->m_InputDim[1])); + clErr |= clSetKernelArg(this->m_PixelCalculation, 9, sizeof(cl_uint), &(this->m_InputDim[2])); + clErr |= clSetKernelArg(this->m_PixelCalculation, 10, sizeof(cl_uint), &(this->m_OutputDim[0])); + clErr |= clSetKernelArg(this->m_PixelCalculation, 11, sizeof(cl_uint), &(this->m_OutputDim[1])); + + CHECK_OCL_ERR(clErr); } + + // execute the filter on a 3D NDRange + if (m_OutputDim[2] == 1) + this->ExecuteKernelChunks(m_PixelCalculation, 2, m_ChunkSize); + else + this->ExecuteKernelChunks(m_PixelCalculation, 3, m_ChunkSize); // 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.Algorithm) { case BeamformingFilter::beamformingSettings::BeamformingAlgorithm::DAS: { if(m_Conf.DelayCalculationMethod == BeamformingFilter::beamformingSettings::DelayCalc::QuadApprox) this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "ckDASQuad", &clErr); else if (m_Conf.DelayCalculationMethod == BeamformingFilter::beamformingSettings::DelayCalc::Spherical) this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "ckDASSphe", &clErr); break; } case BeamformingFilter::beamformingSettings::BeamformingAlgorithm::DMAS: { if (m_Conf.DelayCalculationMethod == BeamformingFilter::beamformingSettings::DelayCalc::QuadApprox) - this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "ckDMASQuad", &clErr); + this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "ckDMAS", &clErr); else if (m_Conf.DelayCalculationMethod == BeamformingFilter::beamformingSettings::DelayCalc::Spherical) - this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "ckDASSphe", &clErr); + this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "ckDMAS", &clErr); break; } default: { this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "ckDASQuad", &clErr); break; } } buildErr |= CHECK_OCL_ERR( clErr ); } return (OclFilter::IsInitialized() && buildErr ); } void mitk::PhotoacousticOCLBeamformingFilter::SetInput(mitk::Image::Pointer image) { OclDataSetToDataSetFilter::SetInput(image); m_InputImage = image; m_InputDim[0] = m_InputImage->GetDimension(0); m_InputDim[1] = m_InputImage->GetDimension(1); m_InputDim[2] = m_InputImage->GetDimension(2); } void mitk::PhotoacousticOCLBeamformingFilter::SetInput(void* data, unsigned int* dimensions, unsigned int BpE) { OclDataSetToDataSetFilter::SetInput(data, dimensions[0] * dimensions[1] * dimensions[2], BpE); m_InputDim[0] = dimensions[0]; m_InputDim[1] = dimensions[1]; m_InputDim[2] = dimensions[2]; } 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->SetGeometry(m_InputImage->GetGeometry()); outputImage->SetImportVolume(pData, 0, 0, mitk::Image::ReferenceMemory); } MITK_DEBUG << "Image Initialized."; return outputImage; } #endif diff --git a/Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticOCLBeamformer.h b/Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticOCLBeamformer.h index dd6a59ef94..59a1d9aec1 100644 --- a/Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticOCLBeamformer.h +++ b/Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticOCLBeamformer.h @@ -1,125 +1,133 @@ /*=================================================================== 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_ #ifdef PHOTOACOUSTICS_USE_GPU #include "mitkOclDataSetToDataSetFilter.h" #include #include "mitkPhotoacousticBeamformingFilter.h" -#include +#include "mitkPhotoacousticOCLDelayCalculation.h" +#include "mitkPhotoacousticOCLMemoryLocSum.h" +#include "mitkPhotoacousticOCLUsedLinesCalculation.h" namespace mitk { /** Documentation * * \brief The OclBinaryThresholdImageFilter computes a binary segmentation based on given threshold values. * * The filter requires two threshold values ( the upper and the lower threshold ) and two image values ( inside and outside ). The resulting voxel of the segmentation image is assigned the inside value 1 if the image value is between the given thresholds and the outside value otherwise. */ class PhotoacousticOCLBeamformingFilter : public OclDataSetToDataSetFilter, public itk::Object { public: mitkClassMacroItkParent(PhotoacousticOCLBeamformingFilter, itk::Object); itkNewMacro(Self); /** * @brief SetInput Set the input data through an image. Arbitrary images are supported */ void SetInput(Image::Pointer image); /** * @brief SetInput Manually set the input data while providing dimensions and memory size of the input data. */ void SetInput(void* data, unsigned int* dimensions, unsigned int BpE); /** * @brief GetOutputAsImage Returns an mitk::Image constructed from the processed data */ mitk::Image::Pointer GetOutputAsImage(); /** Update the filter */ void Update(); /** Set the Output dimensions, which are also used for the openCL global worksize */ void SetOutputDim( unsigned int outputDim[3]) { m_OutputDim[0] = outputDim[0]; m_OutputDim[1] = outputDim[1]; m_OutputDim[2] = outputDim[2]; } /** Set the Apodisation function to apply when beamforming */ void SetApodisation(float* apodisation, unsigned short apodArraySize) { m_ApodArraySize = apodArraySize; m_Apodisation = apodisation; } void SetConfig(BeamformingFilter::beamformingSettings settings) { m_Conf = settings; } protected: /** Constructor */ PhotoacousticOCLBeamformingFilter(); /** Destructor */ virtual ~PhotoacousticOCLBeamformingFilter(); /** Initialize the filter */ bool Initialize(); 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]; unsigned int m_InputDim[3]; float* m_Apodisation; unsigned short m_ApodArraySize; unsigned short m_PAImage; BeamformingFilter::beamformingSettings m_Conf; mitk::Image::Pointer m_InputImage; + + size_t m_ChunkSize[3]; + + mitk::OCLMemoryLocSum::Pointer m_SumFilter; + mitk::OCLUsedLinesCalculation::Pointer m_UsedLinesCalculation; + mitk::OCLDelayCalculation::Pointer m_DelayCalculation; }; } #endif #endif diff --git a/Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticOCLDelayCalculation.cpp b/Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticOCLDelayCalculation.cpp new file mode 100644 index 0000000000..04c0e240b3 --- /dev/null +++ b/Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticOCLDelayCalculation.cpp @@ -0,0 +1,117 @@ +/*=================================================================== + +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 "mitkPhotoacousticOCLDelayCalculation.h" +#include "usServiceReference.h" +#include "mitkImageReadAccessor.h" + +mitk::OCLDelayCalculation::OCLDelayCalculation() + : m_PixelCalculation(NULL) +{ + this->AddSourceFile("DelayCalculation.cl"); + this->m_FilterID = "DelayCalculation"; +} + +mitk::OCLDelayCalculation::~OCLDelayCalculation() +{ + if (this->m_PixelCalculation) + { + clReleaseKernel(m_PixelCalculation); + } +} + +void mitk::OCLDelayCalculation::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::OCLDelayCalculation::Execute() +{ + cl_int clErr = 0; + + unsigned int gridDim[3] = { m_Conf.ReconstructionLines, m_Conf.SamplesPerLine, m_Conf.ReconstructionLines }; + + try + { + this->InitExecNoInput(this->m_PixelCalculation, gridDim, m_BufferSize, sizeof(unsigned short)); + } + catch (const mitk::Exception& e) + { + MITK_ERROR << "Caught exception while initializing filter: " << e.what(); + return; + } + + if (m_Conf.DelayCalculationMethod == BeamformingFilter::beamformingSettings::DelayCalc::QuadApprox) + m_DelayMultiplicatorRaw = pow(1 / (m_Conf.TimeSpacing*m_Conf.SpeedOfSound) * m_Conf.Pitch * m_Conf.TransducerElements / m_Conf.inputDim[0], 2) / 2; + else if (m_Conf.DelayCalculationMethod == BeamformingFilter::beamformingSettings::DelayCalc::Spherical) + m_DelayMultiplicatorRaw = 1 / (m_Conf.TimeSpacing*m_Conf.SpeedOfSound) * (m_Conf.Pitch*m_Conf.TransducerElements); + + m_IsPAImage = m_Conf.Photoacoustic; + + clErr |= clSetKernelArg(this->m_PixelCalculation, 1, sizeof(cl_mem), &m_UsedLines); + clErr |= clSetKernelArg(this->m_PixelCalculation, 2, sizeof(cl_mem), &m_MemLoc); + clErr |= clSetKernelArg(this->m_PixelCalculation, 3, sizeof(cl_uint), &(this->m_Conf.inputDim[0])); + clErr |= clSetKernelArg(this->m_PixelCalculation, 4, sizeof(cl_uint), &(this->m_Conf.inputDim[1])); + clErr |= clSetKernelArg(this->m_PixelCalculation, 5, sizeof(cl_uint), &(this->m_Conf.ReconstructionLines)); + clErr |= clSetKernelArg(this->m_PixelCalculation, 6, sizeof(cl_uint), &(this->m_Conf.SamplesPerLine)); + clErr |= clSetKernelArg(this->m_PixelCalculation, 7, sizeof(cl_char), &(this->m_IsPAImage)); + clErr |= clSetKernelArg(this->m_PixelCalculation, 8, sizeof(cl_float), &(this->m_DelayMultiplicatorRaw)); + + CHECK_OCL_ERR(clErr); + + // execute the filter on a 3D NDRange + this->ExecuteKernel(m_PixelCalculation, 3); + + // signalize the GPU-side data changed + m_Output->Modified(GPU_DATA); +} + +us::Module *mitk::OCLDelayCalculation::GetModule() +{ + return us::GetModuleContext()->GetModule(); +} + +bool mitk::OCLDelayCalculation::Initialize() +{ + bool buildErr = true; + cl_int clErr = 0; + + if (OclFilter::Initialize()) + { + if(m_Conf.DelayCalculationMethod == BeamformingFilter::beamformingSettings::DelayCalc::QuadApprox) + this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "ckDelayCalculationQuad", &clErr); + if (m_Conf.DelayCalculationMethod == BeamformingFilter::beamformingSettings::DelayCalc::Spherical) + this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "ckDelayCalculationSphe", &clErr); + buildErr |= CHECK_OCL_ERR(clErr); + } + return (OclFilter::IsInitialized() && buildErr); +} \ No newline at end of file diff --git a/Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticOCLUsedLinesCalculation.h b/Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticOCLDelayCalculation.h similarity index 76% copy from Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticOCLUsedLinesCalculation.h copy to Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticOCLDelayCalculation.h index 8e3ab5d865..fc8ff48fb3 100644 --- a/Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticOCLUsedLinesCalculation.h +++ b/Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticOCLDelayCalculation.h @@ -1,95 +1,105 @@ /*=================================================================== 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 _MITKPHOTOACOUSTICSTESTOCL_H_ -#define _MITKPHOTOACOUSTICSTESTOCL_H_ +#ifndef _MITKPHOTOACOUSTICSDELAYCALC_H_ +#define _MITKPHOTOACOUSTICSDELAYCALC_H_ #ifdef PHOTOACOUSTICS_USE_GPU #include "mitkOclDataSetToDataSetFilter.h" #include #include namespace mitk { /** Documentation * * \brief The OclBinaryThresholdImageFilter computes a binary segmentation based on given threshold values. * * The filter requires two threshold values ( the upper and the lower threshold ) and two image values ( inside and outside ). The resulting voxel of the segmentation image is assigned the inside value 1 if the image value is between the given thresholds and the outside value otherwise. */ - class OCLUsedLinesCalculation : public OclDataSetToDataSetFilter, public itk::Object + class OCLDelayCalculation : public OclDataSetToDataSetFilter, public itk::Object { public: - mitkClassMacroItkParent(OCLUsedLinesCalculation, itk::Object); + mitkClassMacroItkParent(OCLDelayCalculation, itk::Object); itkNewMacro(Self); /** * @brief SetInput Set the input image. Only 3D images are supported for now. * @param image a 3D image. * @throw mitk::Exception if the dimesion is not 3. */ /** Update the filter */ void Update(); void SetConfig(mitk::BeamformingFilter::beamformingSettings conf) { m_Conf = conf; } + void SetInputs(cl_mem usedLines, cl_mem memoryLocations, unsigned int bufferSize) + { + m_MemLoc = memoryLocations; + m_UsedLines = usedLines; + m_BufferSize = bufferSize; + } + protected: /** Constructor */ - OCLUsedLinesCalculation(); + OCLDelayCalculation(); /** Destructor */ - virtual ~OCLUsedLinesCalculation(); + virtual ~OCLDelayCalculation(); /** Initialize the filter */ bool Initialize(); void Execute(); mitk::PixelType GetOutputType() { return mitk::MakeScalarPixelType(); } int GetBytesPerElem() { return sizeof(unsigned short); } virtual us::Module* GetModule(); int m_sizeThis; - private: /** The OpenCL kernel for the filter */ cl_kernel m_PixelCalculation; mitk::BeamformingFilter::beamformingSettings m_Conf; - float m_part; + cl_mem m_MemLoc; + cl_mem m_UsedLines; + unsigned int m_BufferSize; + float m_DelayMultiplicatorRaw; + char m_IsPAImage; }; } #endif #endif \ No newline at end of file diff --git a/Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticOCLMemoryLocSum.cpp b/Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticOCLMemoryLocSum.cpp new file mode 100644 index 0000000000..566a0c0767 --- /dev/null +++ b/Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticOCLMemoryLocSum.cpp @@ -0,0 +1,117 @@ +/*=================================================================== + +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 "mitkPhotoacousticOCLMemoryLocSum.h" +#include "usServiceReference.h" +#include "mitkImageReadAccessor.h" + +mitk::OCLMemoryLocSum::OCLMemoryLocSum() + : m_PixelCalculation(NULL), m_Sum(0) +{ + this->AddSourceFile("MemoryLocSum.cl"); + this->m_FilterID = "MemoryLocSum"; +} + +mitk::OCLMemoryLocSum::~OCLMemoryLocSum() +{ + if (this->m_PixelCalculation) + { + clReleaseKernel(m_PixelCalculation); + } +} + +void mitk::OCLMemoryLocSum::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::OCLMemoryLocSum::Execute() +{ + cl_int clErr = 0; + + unsigned int gridDim[3] = { m_Dim[0], m_Dim[1], 1 }; + size_t outputSize = gridDim[0] * gridDim[1] * 1; + + try + { + this->InitExec(this->m_PixelCalculation, gridDim, outputSize, sizeof(unsigned int)); + } + catch (const mitk::Exception& e) + { + MITK_ERROR << "Caught exception while initializing filter: " << e.what(); + return; + } + + us::ServiceReference ref = GetModuleContext()->GetServiceReference(); + OclResourceService* resources = GetModuleContext()->GetService(ref); + cl_context gpuContext = resources->GetContext(); + + cl_mem cl_Sum = clCreateBuffer(gpuContext, CL_MEM_WRITE_ONLY, sizeof(unsigned int), NULL, &clErr); + CHECK_OCL_ERR(clErr); + + clErr = clSetKernelArg(this->m_PixelCalculation, 2, sizeof(cl_mem), &cl_Sum); + + CHECK_OCL_ERR(clErr); + + // execute the filter on a 3D NDRange + this->ExecuteKernel(m_PixelCalculation, 2); + + // signalize the GPU-side data changed + m_Output->Modified(GPU_DATA); + + char* data = new char[1 * sizeof(unsigned int)]; + + clErr = clEnqueueReadBuffer(m_CommandQue, cl_Sum, CL_FALSE, 0, 1 * sizeof(unsigned int), data, 0, nullptr, nullptr); + CHECK_OCL_ERR(clErr); + + clFlush(m_CommandQue); + + m_Sum = ((unsigned int*)data)[0]; +} + +us::Module *mitk::OCLMemoryLocSum::GetModule() +{ + return us::GetModuleContext()->GetModule(); +} + +bool mitk::OCLMemoryLocSum::Initialize() +{ + bool buildErr = true; + cl_int clErr = 0; + + if (OclFilter::Initialize()) + { + this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "ckMemoryLocSum", &clErr); + buildErr |= CHECK_OCL_ERR(clErr); + } + return (OclFilter::IsInitialized() && buildErr); +} \ No newline at end of file diff --git a/Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticOCLUsedLinesCalculation.h b/Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticOCLMemoryLocSum.h similarity index 77% copy from Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticOCLUsedLinesCalculation.h copy to Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticOCLMemoryLocSum.h index 8e3ab5d865..bce56edbc2 100644 --- a/Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticOCLUsedLinesCalculation.h +++ b/Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticOCLMemoryLocSum.h @@ -1,95 +1,99 @@ /*=================================================================== 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 _MITKPHOTOACOUSTICSTESTOCL_H_ -#define _MITKPHOTOACOUSTICSTESTOCL_H_ +#ifndef _MITKPHOTOACOUSTICSPARALLELSUM_H_ +#define _MITKPHOTOACOUSTICSPARALLELSUM_H_ #ifdef PHOTOACOUSTICS_USE_GPU #include "mitkOclDataSetToDataSetFilter.h" #include -#include namespace mitk { /** Documentation * * \brief The OclBinaryThresholdImageFilter computes a binary segmentation based on given threshold values. * * The filter requires two threshold values ( the upper and the lower threshold ) and two image values ( inside and outside ). The resulting voxel of the segmentation image is assigned the inside value 1 if the image value is between the given thresholds and the outside value otherwise. */ - class OCLUsedLinesCalculation : public OclDataSetToDataSetFilter, public itk::Object + class OCLMemoryLocSum : public OclDataSetToDataSetFilter, public itk::Object { public: - mitkClassMacroItkParent(OCLUsedLinesCalculation, itk::Object); + mitkClassMacroItkParent(OCLMemoryLocSum, itk::Object); itkNewMacro(Self); /** * @brief SetInput Set the input image. Only 3D images are supported for now. * @param image a 3D image. * @throw mitk::Exception if the dimesion is not 3. */ /** Update the filter */ void Update(); - void SetConfig(mitk::BeamformingFilter::beamformingSettings conf) + void SetInputDimensions(unsigned int* dimensions) { - m_Conf = conf; + m_Dim[0] = dimensions[0]; + m_Dim[1] = dimensions[1]; + } + + unsigned int GetSum() + { + return m_Sum; } protected: /** Constructor */ - OCLUsedLinesCalculation(); + OCLMemoryLocSum(); /** Destructor */ - virtual ~OCLUsedLinesCalculation(); + virtual ~OCLMemoryLocSum(); /** Initialize the filter */ bool Initialize(); void Execute(); mitk::PixelType GetOutputType() { return mitk::MakeScalarPixelType(); } int GetBytesPerElem() { return sizeof(unsigned short); } virtual us::Module* GetModule(); int m_sizeThis; - private: /** The OpenCL kernel for the filter */ cl_kernel m_PixelCalculation; - mitk::BeamformingFilter::beamformingSettings m_Conf; - float m_part; + unsigned int m_Sum; + unsigned int m_Dim[3]; }; } #endif #endif \ No newline at end of file diff --git a/Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticOCLUsedLinesCalculation.h b/Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticOCLUsedLinesCalculation.h index 8e3ab5d865..c8c681db67 100644 --- a/Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticOCLUsedLinesCalculation.h +++ b/Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticOCLUsedLinesCalculation.h @@ -1,95 +1,95 @@ /*=================================================================== 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 _MITKPHOTOACOUSTICSTESTOCL_H_ -#define _MITKPHOTOACOUSTICSTESTOCL_H_ +#ifndef _MITKPHOTOACOUSTICSOCLUSEDLINESCALCULATION_H_ +#define _MITKPHOTOACOUSTICSOCLUSEDLINESCALCULATION_H_ #ifdef PHOTOACOUSTICS_USE_GPU #include "mitkOclDataSetToDataSetFilter.h" #include #include namespace mitk { /** Documentation * * \brief The OclBinaryThresholdImageFilter computes a binary segmentation based on given threshold values. * * The filter requires two threshold values ( the upper and the lower threshold ) and two image values ( inside and outside ). The resulting voxel of the segmentation image is assigned the inside value 1 if the image value is between the given thresholds and the outside value otherwise. */ class OCLUsedLinesCalculation : public OclDataSetToDataSetFilter, public itk::Object { public: mitkClassMacroItkParent(OCLUsedLinesCalculation, itk::Object); itkNewMacro(Self); /** * @brief SetInput Set the input image. Only 3D images are supported for now. * @param image a 3D image. * @throw mitk::Exception if the dimesion is not 3. */ /** Update the filter */ void Update(); void SetConfig(mitk::BeamformingFilter::beamformingSettings conf) { m_Conf = conf; } protected: /** Constructor */ OCLUsedLinesCalculation(); /** Destructor */ virtual ~OCLUsedLinesCalculation(); /** Initialize the filter */ bool Initialize(); void Execute(); mitk::PixelType GetOutputType() { return mitk::MakeScalarPixelType(); } int GetBytesPerElem() { return sizeof(unsigned short); } virtual us::Module* GetModule(); int m_sizeThis; private: /** The OpenCL kernel for the filter */ cl_kernel m_PixelCalculation; mitk::BeamformingFilter::beamformingSettings m_Conf; float m_part; }; } #endif #endif \ No newline at end of file diff --git a/Modules/PhotoacousticsAlgorithms/Algorithms/mitkPhotoacousticBeamformingFilter.cpp b/Modules/PhotoacousticsAlgorithms/Algorithms/mitkPhotoacousticBeamformingFilter.cpp index 8a4321f466..0fc0d0bc7a 100644 --- a/Modules/PhotoacousticsAlgorithms/Algorithms/mitkPhotoacousticBeamformingFilter.cpp +++ b/Modules/PhotoacousticsAlgorithms/Algorithms/mitkPhotoacousticBeamformingFilter.cpp @@ -1,542 +1,546 @@ /*=================================================================== mitkPhotoacousticBeamformingFilter 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 "mitkPhotoacousticBeamformingFilter.h" #include "mitkProperties.h" #include "mitkImageReadAccessor.h" #include #include #include #include #include #include #include "mitkImageCast.h" -#include +#include "./OpenCLFilter/mitkPhotoacousticOCLBeamformer.h" +#include "./OpenCLFilter/mitkPhotoacousticOCLDelayCalculation.h" +#include "./OpenCLFilter/mitkPhotoacousticOCLMemoryLocSum.h" +#include "./OpenCLFilter/mitkPhotoacousticOCLUsedLinesCalculation.h" mitk::BeamformingFilter::BeamformingFilter() : m_OutputData(nullptr), m_InputData(nullptr) { this->SetNumberOfIndexedInputs(1); this->SetNumberOfRequiredInputs(1); m_ProgressHandle = [](int, std::string) {}; } void mitk::BeamformingFilter::SetProgressHandle(std::function progressHandle) { m_ProgressHandle = progressHandle; } mitk::BeamformingFilter::~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(); //GenerateTimeInInputRegion(output, input); } 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; itkDebugMacro(<< "GenerateOutputInformation()"); unsigned int dim[] = { m_Conf.ReconstructionLines, m_Conf.SamplesPerLine, input->GetDimension(2) }; output->Initialize(mitk::MakeScalarPixelType(), 3, dim); mitk::Vector3D spacing; spacing[0] = m_Conf.Pitch * m_Conf.TransducerElements * 1000 / m_Conf.ReconstructionLines; spacing[1] = m_Conf.RecordTime / 2 * m_Conf.SpeedOfSound * 1000 / m_Conf.SamplesPerLine; spacing[2] = 1; output->GetGeometry()->SetSpacing(spacing); output->GetGeometry()->Modified(); output->SetPropertyList(input->GetPropertyList()->Clone()); m_TimeOfHeaderInitialization.Modified(); } + void mitk::BeamformingFilter::GenerateData() { GenerateOutputInformation(); mitk::Image::Pointer input = this->GetInput(); mitk::Image::Pointer output = this->GetOutput(); if (!output->IsInitialized()) return; - const int apodArraySize = m_Conf.TransducerElements * 4; // set the resolution of the apodization array + const int apodArraySize = m_Conf.TransducerElements; // set the resolution of the apodization array float* ApodWindow; // calculate the appropiate apodization window switch (m_Conf.Apod) { case beamformingSettings::Apodization::Hann: ApodWindow = VonHannFunction(apodArraySize); break; case beamformingSettings::Apodization::Hamm: ApodWindow = HammFunction(apodArraySize); break; case beamformingSettings::Apodization::Box: ApodWindow = BoxFunction(apodArraySize); break; default: ApodWindow = BoxFunction(apodArraySize); break; } auto begin = std::chrono::high_resolution_clock::now(); // debbuging the performance... if (!m_Conf.UseGPU) { 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)); // first, we check whether the dara is float, other formats are unsupported if (input->GetPixelType().GetTypeAsString() == "scalar (float)" || input->GetPixelType().GetTypeAsString() == " (float)") { m_InputData = (float*)inputReadAccessor.GetData(); } else { MITK_INFO << "Pixel type is not float, abort"; return; } m_OutputData = new float[m_Conf.ReconstructionLines*m_Conf.SamplesPerLine]; // 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.Algorithm == beamformingSettings::BeamformingAlgorithm::DAS) { if (m_Conf.DelayCalculationMethod == beamformingSettings::DelayCalc::QuadApprox) { for (short line = 0; line < outputDim[0]; ++line) { threads[line] = std::thread(&BeamformingFilter::DASQuadraticLine, this, m_InputData, m_OutputData, inputDim, outputDim, line, ApodWindow, apodArraySize); } } else if (m_Conf.DelayCalculationMethod == beamformingSettings::DelayCalc::Spherical) { for (short line = 0; line < outputDim[0]; ++line) { threads[line] = std::thread(&BeamformingFilter::DASSphericalLine, this, m_InputData, m_OutputData, inputDim, outputDim, line, ApodWindow, apodArraySize); } } } else if (m_Conf.Algorithm == beamformingSettings::BeamformingAlgorithm::DMAS) { if (m_Conf.DelayCalculationMethod == beamformingSettings::DelayCalc::QuadApprox) { for (short line = 0; line < outputDim[0]; ++line) { threads[line] = std::thread(&BeamformingFilter::DMASQuadraticLine, this, m_InputData, m_OutputData, inputDim, outputDim, line, ApodWindow, apodArraySize); } } else if (m_Conf.DelayCalculationMethod == beamformingSettings::DelayCalc::Spherical) { for (short line = 0; line < outputDim[0]; ++line) { threads[line] = std::thread(&BeamformingFilter::DMASSphericalLine, this, m_InputData, m_OutputData, inputDim, outputDim, line, ApodWindow, apodArraySize); } } } // 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; } } #ifdef PHOTOACOUSTICS_USE_GPU 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"; return; } m_ProgressHandle(50, "performing reconstruction"); mitk::PhotoacousticOCLBeamformingFilter::Pointer m_oclFilter = mitk::PhotoacousticOCLBeamformingFilter::New(); m_oclFilter->SetApodisation(ApodWindow, apodArraySize); m_oclFilter->SetConfig(m_Conf); m_oclFilter->SetOutputDim(output->GetDimensions()); m_oclFilter->SetInput(input); m_oclFilter->Update(); void* out = m_oclFilter->GetOutput(); output->SetImportVolume(out, 0, 0, mitk::Image::ReferenceMemory); } catch (mitk::Exception &e) { std::string errorMessage = "Caught unexpected exception "; errorMessage.append(e.what()); MITK_ERROR << errorMessage; } } #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; } void mitk::BeamformingFilter::Configure(beamformingSettings settings) { m_Conf = settings; } float* mitk::BeamformingFilter::VonHannFunction(int samples) { float* ApodWindow = new float[samples]; for (int n = 0; n < samples; ++n) { ApodWindow[n] = (1 - cos(2 * M_PI * n / (samples - 1))) / 2; } return ApodWindow; } float* mitk::BeamformingFilter::HammFunction(int samples) { float* ApodWindow = new float[samples]; for (int n = 0; n < samples; ++n) { ApodWindow[n] = 0.54 - 0.46*cos(2 * M_PI*n / (samples - 1)); } return ApodWindow; } float* mitk::BeamformingFilter::BoxFunction(int samples) { float* ApodWindow = new float[samples]; for (int n = 0; n < samples; ++n) { ApodWindow[n] = 1; } return ApodWindow; } void mitk::BeamformingFilter::DASQuadraticLine(float* input, float* output, float inputDim[2], float outputDim[2], const short& line, float* apodisation, const short& apodArraySize) { float& inputS = inputDim[1]; float& inputL = inputDim[0]; float& outputS = outputDim[1]; float& outputL = outputDim[0]; short AddSample = 0; short maxLine = 0; short minLine = 0; float delayMultiplicator = 0; float l_i = 0; float s_i = 0; float part = 0.07 * inputL; float tan_phi = std::tan(m_Conf.Angle / 360 * 2 * M_PI); float part_multiplicator = tan_phi * m_Conf.TimeSpacing * m_Conf.SpeedOfSound / m_Conf.Pitch * inputL / m_Conf.TransducerElements; float apod_mult = 1; short usedLines = (maxLine - minLine); //quadratic delay l_i = line / outputL * inputL; for (short sample = 0; sample < outputS; ++sample) { s_i = (float)sample / outputS * inputS / 2; part = part_multiplicator*s_i; if (part < 1) part = 1; maxLine = (short)std::min((l_i + part) + 1, inputL); minLine = (short)std::max((l_i - part), 0.0f); usedLines = (maxLine - minLine); - apod_mult = apodArraySize / (maxLine - minLine); + apod_mult = (float)apodArraySize / (float)usedLines; delayMultiplicator = pow((1 / (m_Conf.TimeSpacing*m_Conf.SpeedOfSound) * (m_Conf.Pitch*m_Conf.TransducerElements) / inputL), 2) / s_i / 2; for (short l_s = minLine; l_s < maxLine; ++l_s) { AddSample = delayMultiplicator * pow((l_s - l_i), 2) + s_i + (1 - m_Conf.Photoacoustic)*s_i; if (AddSample < inputS && AddSample >= 0) output[sample*(short)outputL + line] += input[l_s + AddSample*(short)inputL] * apodisation[(short)((l_s - minLine)*apod_mult)]; else --usedLines; } output[sample*(short)outputL + line] = output[sample*(short)outputL + line] / usedLines; } } void mitk::BeamformingFilter::DASSphericalLine(float* input, float* output, float inputDim[2], float outputDim[2], const short& line, float* apodisation, const short& apodArraySize) { float& inputS = inputDim[1]; float& inputL = inputDim[0]; float& outputS = outputDim[1]; float& outputL = outputDim[0]; short AddSample = 0; short maxLine = 0; short minLine = 0; float l_i = 0; float s_i = 0; float part = 0.07 * inputL; float tan_phi = std::tan(m_Conf.Angle / 360 * 2 * M_PI); float part_multiplicator = tan_phi * m_Conf.TimeSpacing * m_Conf.SpeedOfSound / m_Conf.Pitch * inputL / m_Conf.TransducerElements; float apod_mult = 1; short usedLines = (maxLine - minLine); //exact delay l_i = (float)line / outputL * inputL; for (short sample = 0; sample < outputS; ++sample) { s_i = (float)sample / outputS * inputS / 2; part = part_multiplicator*s_i; if (part < 1) part = 1; maxLine = (short)std::min((l_i + part) + 1, inputL); minLine = (short)std::max((l_i - part), 0.0f); usedLines = (maxLine - minLine); - apod_mult = apodArraySize / (maxLine - minLine); + apod_mult = (float)apodArraySize / (float)usedLines; for (short l_s = minLine; l_s < maxLine; ++l_s) { AddSample = (int)sqrt( pow(s_i, 2) + pow((1 / (m_Conf.TimeSpacing*m_Conf.SpeedOfSound) * ((l_s - l_i)*m_Conf.Pitch*m_Conf.TransducerElements) / inputL), 2) ) + (1 - m_Conf.Photoacoustic)*s_i; if (AddSample < inputS && AddSample >= 0) output[sample*(short)outputL + line] += input[l_s + AddSample*(short)inputL] * apodisation[(short)((l_s - minLine)*apod_mult)]; else --usedLines; } output[sample*(short)outputL + line] = output[sample*(short)outputL + line] / usedLines; } } void mitk::BeamformingFilter::DMASQuadraticLine(float* input, float* output, float inputDim[2], float outputDim[2], const short& line, float* apodisation, const short& apodArraySize) { float& inputS = inputDim[1]; float& inputL = inputDim[0]; float& outputS = outputDim[1]; float& outputL = outputDim[0]; short maxLine = 0; short minLine = 0; float delayMultiplicator = 0; float l_i = 0; float s_i = 0; float part = 0.07 * inputL; float tan_phi = std::tan(m_Conf.Angle / 360 * 2 * M_PI); float part_multiplicator = tan_phi * m_Conf.TimeSpacing * m_Conf.SpeedOfSound / m_Conf.Pitch * inputL / m_Conf.TransducerElements; float apod_mult = 1; float mult = 0; short usedLines = (maxLine - minLine); //quadratic delay l_i = line / outputL * inputL; for (short sample = 0; sample < outputS; ++sample) { s_i = sample / outputS * inputS / 2; part = part_multiplicator*s_i; if (part < 1) part = 1; maxLine = (short)std::min((l_i + part) + 1, inputL); minLine = (short)std::max((l_i - part), 0.0f); usedLines = (maxLine - minLine); - apod_mult = apodArraySize / (maxLine - minLine); + apod_mult = (float)apodArraySize / (float)usedLines; delayMultiplicator = pow((1 / (m_Conf.TimeSpacing*m_Conf.SpeedOfSound) * (m_Conf.Pitch*m_Conf.TransducerElements) / inputL), 2) / s_i / 2; //calculate the AddSamples beforehand to save some time short* AddSample = new short[maxLine - minLine]; for (short l_s = 0; l_s < maxLine - minLine; ++l_s) { AddSample[l_s] = (short)(delayMultiplicator * pow((minLine + l_s - l_i), 2) + s_i) + (1 - m_Conf.Photoacoustic)*s_i; } for (short l_s1 = minLine; l_s1 < maxLine - 1; ++l_s1) { if (AddSample[l_s1 - minLine] < (short)inputS && AddSample[l_s1 - minLine] >= 0) { for (short l_s2 = l_s1 + 1; l_s2 < maxLine; ++l_s2) { if (AddSample[l_s2 - minLine] < inputS && AddSample[l_s2 - minLine] >= 0) { mult = input[l_s2 + AddSample[l_s2 - minLine] * (short)inputL] * apodisation[(short)((l_s2 - minLine)*apod_mult)] * input[l_s1 + AddSample[l_s1 - minLine] * (short)inputL] * apodisation[(short)((l_s1 - minLine)*apod_mult)]; output[sample*(short)outputL + line] += sqrt(abs(mult)) * ((mult > 0) - (mult < 0)); } } } else --usedLines; } output[sample*(short)outputL + line] = 10 * output[sample*(short)outputL + line] / (pow(usedLines, 2) - (usedLines - 1)); delete[] AddSample; } } void mitk::BeamformingFilter::DMASSphericalLine(float* input, float* output, float inputDim[2], float outputDim[2], const short& line, float* apodisation, const short& apodArraySize) { float& inputS = inputDim[1]; float& inputL = inputDim[0]; float& outputS = outputDim[1]; float& outputL = outputDim[0]; short maxLine = 0; short minLine = 0; float l_i = 0; float s_i = 0; float part = 0.07 * inputL; float tan_phi = std::tan(m_Conf.Angle / 360 * 2 * M_PI); float part_multiplicator = tan_phi * m_Conf.TimeSpacing * m_Conf.SpeedOfSound / m_Conf.Pitch * inputL / m_Conf.TransducerElements; float apod_mult = 1; float mult = 0; short usedLines = (maxLine - minLine); //exact delay l_i = line / outputL * inputL; for (short sample = 0; sample < outputS; ++sample) { s_i = sample / outputS * inputS / 2; part = part_multiplicator*s_i; if (part < 1) part = 1; maxLine = (short)std::min((l_i + part) + 1, inputL); minLine = (short)std::max((l_i - part), 0.0f); usedLines = (maxLine - minLine); - apod_mult = apodArraySize / (maxLine - minLine); + apod_mult = (float)apodArraySize / (float)usedLines; //calculate the AddSamples beforehand to save some time short* AddSample = new short[maxLine - minLine]; for (short l_s = 0; l_s < maxLine - minLine; ++l_s) { AddSample[l_s] = (short)sqrt( pow(s_i, 2) + pow((1 / (m_Conf.TimeSpacing*m_Conf.SpeedOfSound) * ((minLine + l_s - l_i)*m_Conf.Pitch*m_Conf.TransducerElements) / inputL), 2) ) + (1 - m_Conf.Photoacoustic)*s_i; } for (short l_s1 = minLine; l_s1 < maxLine - 1; ++l_s1) { if (AddSample[l_s1 - minLine] < inputS && AddSample[l_s1 - minLine] >= 0) { for (short l_s2 = l_s1 + 1; l_s2 < maxLine; ++l_s2) { if (AddSample[l_s2 - minLine] < inputS && AddSample[l_s2 - minLine] >= 0) { mult = input[l_s2 + AddSample[l_s2 - minLine] * (short)inputL] * apodisation[(int)((l_s2 - minLine)*apod_mult)] * input[l_s1 + AddSample[l_s1 - minLine] * (short)inputL] * apodisation[(int)((l_s1 - minLine)*apod_mult)]; output[sample*(short)outputL + line] += sqrt(abs(mult)) * ((mult > 0) - (mult < 0)); } } } else --usedLines; } output[sample*(short)outputL + line] = 10 * output[sample*(short)outputL + line] / (pow(usedLines, 2) - (usedLines - 1)); delete[] AddSample; } } diff --git a/Modules/PhotoacousticsAlgorithms/Resources/DASQuadratic.cl b/Modules/PhotoacousticsAlgorithms/Resources/DASQuadratic.cl index fe873c7b42..12a1456807 100644 --- a/Modules/PhotoacousticsAlgorithms/Resources/DASQuadratic.cl +++ b/Modules/PhotoacousticsAlgorithms/Resources/DASQuadratic.cl @@ -1,73 +1,72 @@ /*=================================================================== 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. ===================================================================*/ __kernel void ckDASQuad( __global float* dSource, // input image __global float* dDest, // output buffer __global float* apodArray, unsigned short apodArraySize, float SpeedOfSound, float TimeSpacing, float Pitch, float Angle, unsigned short PAImage, unsigned short TransducerElements, unsigned int inputL, unsigned int inputS, unsigned int Slices, unsigned int outputL, unsigned int outputS // parameters ) { // get thread identifier unsigned int globalPosX = get_global_id(0); unsigned int globalPosY = get_global_id(1); unsigned int globalPosZ = get_global_id(2); // terminate non-valid threads if ( globalPosX < outputL && globalPosY < outputS && globalPosZ < Slices ) { float l_i = (float)globalPosX / outputL * inputL; float s_i = (float)globalPosY / outputS * inputS / 2; float part = (tan(Angle / 360 * 2 * M_PI) * ((SpeedOfSound * TimeSpacing) * s_i) / (Pitch * TransducerElements)) * inputL; if (part < 1) part = 1; short maxLine = min((l_i + part) + 1, (float)inputL); short minLine = max((l_i - part), 0.0f); short usedLines = (maxLine - minLine); - float apod_mult = apodArraySize / (maxLine - minLine); + float apod_mult = (float)apodArraySize / (float)usedLines; short AddSample = 0; - float output = 0; float delayMultiplicator = pow(1 / (TimeSpacing*SpeedOfSound) * Pitch * TransducerElements / inputL, 2) / s_i / 2; for (short l_s = minLine; l_s < maxLine; ++l_s) { AddSample = delayMultiplicator * pow((l_s - l_i), 2) + s_i + (1-PAImage)*s_i; if (AddSample < inputS && AddSample >= 0) { - output += apodArray[(short)((l_s - minLine)*apod_mult)] * dSource[(unsigned int)(globalPosZ * inputL * inputS + AddSample * inputL + l_s)]; + output += apodArray[(int)((l_s - minLine)*apod_mult)] * dSource[(int)(globalPosZ * inputL * inputS + AddSample * inputL + l_s)]; } else --usedLines; } dDest[ globalPosZ * outputL * outputS + globalPosY * outputL + globalPosX ] = output / usedLines; } } \ No newline at end of file diff --git a/Modules/PhotoacousticsAlgorithms/Resources/DASspherical.cl b/Modules/PhotoacousticsAlgorithms/Resources/DASspherical.cl index f6974aa5a4..7cedb732a0 100644 --- a/Modules/PhotoacousticsAlgorithms/Resources/DASspherical.cl +++ b/Modules/PhotoacousticsAlgorithms/Resources/DASspherical.cl @@ -1,73 +1,73 @@ /*=================================================================== 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. ===================================================================*/ __kernel void ckDASSphe( __global float* dSource, // input image __global float* dDest, // output buffer __global float* apodArray, unsigned short apodArraySize, float SpeedOfSound, float TimeSpacing, float Pitch, float Angle, unsigned short PAImage, unsigned short TransducerElements, unsigned int inputL, unsigned int inputS, unsigned int Slices, unsigned int outputL, unsigned int outputS // parameters ) { // get thread identifier unsigned int globalPosX = get_global_id(0); unsigned int globalPosY = get_global_id(1); unsigned int globalPosZ = get_global_id(2); // terminate non-valid threads if ( globalPosX < outputL && globalPosY < outputS && globalPosZ < Slices ) { float l_i = (float)globalPosX / outputL * inputL; float s_i = (float)globalPosY / outputS * inputS / 2; float part = (tan(Angle / 360 * 2 * M_PI) * ((SpeedOfSound * TimeSpacing) * s_i) / (Pitch * TransducerElements)) * inputL; if (part < 1) part = 1; short maxLine = min((l_i + part) + 1, (float)inputL); short minLine = max((l_i - part), 0.0f); short usedLines = (maxLine - minLine); - float apod_mult = apodArraySize / (maxLine - minLine); + float apod_mult = (float)apodArraySize / (float)usedLines; short AddSample = 0; float output = 0; for (short l_s = minLine; l_s < maxLine; ++l_s) { AddSample = sqrt( pow(s_i, 2) + pow((1 / (TimeSpacing*SpeedOfSound) * ((l_s - l_i)*Pitch*TransducerElements) / inputL), 2) ) + (1-PAImage)*s_i; if (AddSample < inputS && AddSample >= 0) output += apodArray[(short)((l_s - minLine)*apod_mult)] * dSource[(unsigned int)(globalPosZ * inputL * inputS + AddSample * inputL + l_s)]; else --usedLines; } dDest[ globalPosZ * outputL * outputS + globalPosY * outputL + globalPosX ] = output / usedLines; } } \ No newline at end of file diff --git a/Modules/PhotoacousticsAlgorithms/Resources/DMAS.cl b/Modules/PhotoacousticsAlgorithms/Resources/DMAS.cl new file mode 100644 index 0000000000..1d0b707658 --- /dev/null +++ b/Modules/PhotoacousticsAlgorithms/Resources/DMAS.cl @@ -0,0 +1,75 @@ +/*=================================================================== + +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. + +===================================================================*/ + +__kernel void ckDMAS( + __global float* dSource, // input image + __global float* dDest, // output buffer + __global unsigned short* usedLines, + __global unsigned int* memoryLocations, + __global unsigned short* AddSamples, + __global float* apodArray, + unsigned short apodArraySize, + unsigned int inputL, + unsigned int inputS, + unsigned int Slices, + unsigned int outputL, + unsigned int outputS // parameters +) +{ + // get thread identifier + unsigned int globalPosX = get_global_id(0); + unsigned int globalPosY = get_global_id(1); + unsigned int globalPosZ = get_global_id(2); + + // terminate non-valid threads + if ( globalPosX < outputL && globalPosY < outputS && globalPosZ < Slices ) + { + unsigned short minLine = usedLines[globalPosY * outputL + globalPosX + 1 * outputL * outputS]; + unsigned short maxLine = usedLines[globalPosY * outputL + globalPosX + 2 * outputL * outputS]; + unsigned short curUsedLines = usedLines[globalPosY * outputL + globalPosX]; + + float apod_mult = (float)apodArraySize / (float)curUsedLines; + + unsigned short AddSample1 = 0; + unsigned short AddSample2 = 0; + + float output = 0; + float mult = 0; + + for (short l_s1 = minLine; l_s1 < maxLine; ++l_s1) + { + AddSample1 = AddSamples[memoryLocations[globalPosY * outputL + globalPosX] + l_s1 - minLine]; + if (AddSample1 < inputS && AddSample1 >= 0) { + for (short l_s2 = l_s1 + 1; l_s2 < maxLine; ++l_s2) + { + AddSample2 = AddSamples[memoryLocations[globalPosY * outputL + globalPosX] + l_s2 - minLine]; + if (AddSample1 < inputS && AddSample1 >= 0) { + mult = apodArray[(int)((l_s2 - minLine)*apod_mult)] * + dSource[(int)(globalPosZ * inputL * inputS + AddSample2 * inputL + l_s2)] + * apodArray[(int)((l_s1 - minLine)*apod_mult)] * + dSource[(int)(globalPosZ * inputL * inputS + AddSample1 * inputL + l_s1)]; + + output += sqrt(mult * ((float)(mult>0)-(float)(mult<0))) * ((mult > 0) - (mult < 0)); + } + } + } + else + --curUsedLines; + } + + dDest[ globalPosZ * outputL * outputS + globalPosY * outputL + globalPosX ] = output / (pow((float)curUsedLines, 2.0f) - (curUsedLines - 1)); + } +} \ No newline at end of file diff --git a/Modules/PhotoacousticsAlgorithms/Resources/DelayCalculation.cl b/Modules/PhotoacousticsAlgorithms/Resources/DelayCalculation.cl new file mode 100644 index 0000000000..ae175db70d --- /dev/null +++ b/Modules/PhotoacousticsAlgorithms/Resources/DelayCalculation.cl @@ -0,0 +1,74 @@ +/*=================================================================== + +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. + +===================================================================*/ + +__kernel void ckDelayCalculationQuad( __global unsigned short *gDest, + __global unsigned short *usedLines, + __global unsigned int *memoryLocations, + unsigned int inputL, + unsigned int inputS, + unsigned int outputL, + unsigned int outputS, + char isPAImage, + float delayMultiplicatorRaw // parameters + ) + { + uint globalPosX = get_global_id(0); + uint globalPosY = get_global_id(1); + uint globalPosZ = get_global_id(2); + + if(globalPosZ < usedLines[globalPosY * outputL + globalPosX]) + { + float l_i = (float)globalPosX / outputL * inputL; + float s_i = (float)globalPosY / outputS * inputS / 2; + + float l_s = usedLines[outputL * outputS + globalPosY * outputL + globalPosX] + globalPosZ; + + float delayMultiplicator = delayMultiplicatorRaw / s_i; + float AddSample = delayMultiplicator * pow((l_s - l_i), 2) + s_i + (1-isPAImage)*s_i; + gDest[memoryLocations[globalPosY * outputL + globalPosX] + globalPosZ] = AddSample; + } + } + + __kernel void ckDelayCalculationSphe( __global unsigned short *gDest, + __global unsigned short *usedLines, + __global unsigned int *memoryLocations, + unsigned int inputL, + unsigned int inputS, + unsigned int outputL, + unsigned int outputS, + char isPAImage, + float delayMultiplicatorRaw // parameters + ) + { + uint globalPosX = get_global_id(0); + uint globalPosY = get_global_id(1); + uint globalPosZ = get_global_id(2); + + if(globalPosZ < usedLines[globalPosY * outputL + globalPosX]) + { + float l_i = (float)globalPosX / outputL * inputL; + float s_i = (float)globalPosY / outputS * inputS / 2; + + float l_s = usedLines[outputL * outputS + globalPosY * outputL + globalPosX] + globalPosZ; + + gDest[memoryLocations[globalPosY * outputL + globalPosX] + globalPosZ] = + sqrt( + pow(s_i, 2) + + + pow((delayMultiplicatorRaw * ((l_s - l_i)) / inputL), 2) + ) + (1-isPAImage)*s_i; + } + } \ No newline at end of file diff --git a/Modules/PhotoacousticsAlgorithms/Resources/MemoryLocSum.cl b/Modules/PhotoacousticsAlgorithms/Resources/MemoryLocSum.cl new file mode 100644 index 0000000000..d25d9a31fd --- /dev/null +++ b/Modules/PhotoacousticsAlgorithms/Resources/MemoryLocSum.cl @@ -0,0 +1,35 @@ +/*=================================================================== + +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. + +===================================================================*/ + +__kernel void ckMemoryLocSum( __global unsigned short *input, + __global unsigned int *sums, + __global unsigned int *finalSum + ) + { + uint id = get_global_id(0) + get_global_id(1) * get_global_size(0); + + uint sum = 0; + + for (uint pos = 0; pos < id; ++pos) + { + sum += input[pos]; + } + + sums[id] = sum; + + if (id == (get_global_size(0) - 1) + (get_global_size(1) - 1) * get_global_size(0)) + finalSum[0] = sum + input[id]; + } \ No newline at end of file diff --git a/Modules/PhotoacousticsAlgorithms/Resources/ParallelSum.cl b/Modules/PhotoacousticsAlgorithms/Resources/ParallelSum.cl new file mode 100644 index 0000000000..8d640492b9 --- /dev/null +++ b/Modules/PhotoacousticsAlgorithms/Resources/ParallelSum.cl @@ -0,0 +1,36 @@ +/*=================================================================== + +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. + +===================================================================*/ + +__kernel void ckParallelSum( __const const unsigned short *input, + __global unsigned int *partialSums, + __local unsigned int *localSums) + { + uint local_id = get_local_id(0) + get_local_id(1) * get_local_size(2); + uint group_size = get_local_size(0) * get_local_size(1); + + localSums[local_id] = input[get_global_id(0) + get_global_id(1) * get_global_size(0)]; + + for (uint stride = group_size/2; stride>0; stride /=2) + { + barrier(CLK_LOCAL_MEM_FENCE); + + if (local_id < stride) + localSums[local_id] += localSums[local_id + stride]; + } + + if (local_id == 0) + partialSums[get_group_id(0)] = localSums[0]; + } \ No newline at end of file diff --git a/Modules/PhotoacousticsAlgorithms/Resources/UsedLinesCalculation.cl b/Modules/PhotoacousticsAlgorithms/Resources/UsedLinesCalculation.cl index 3b7ea59743..d284511dd4 100644 --- a/Modules/PhotoacousticsAlgorithms/Resources/UsedLinesCalculation.cl +++ b/Modules/PhotoacousticsAlgorithms/Resources/UsedLinesCalculation.cl @@ -1,48 +1,48 @@ /*=================================================================== 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. ===================================================================*/ __kernel void ckUsedLines( __global unsigned short* dDest, // output buffer float partMult, unsigned int inputL, unsigned int inputS ) { // get thread identifier unsigned int globalPosX = get_global_id(0); unsigned int globalPosY = get_global_id(1); unsigned short outputS = get_global_size(1); unsigned short outputL = get_global_size(0); // terminate non-valid threads if ( globalPosX < outputL && globalPosY < outputS) { float l_i = (float)globalPosX / outputL * inputL; float s_i = (float)globalPosY / outputS * inputS / 2; float part = partMult * s_i; if (part < 1) part = 1; unsigned short maxLine = min((l_i + part) + 1, (float)inputL); unsigned short minLine = max((l_i - part), 0.0f); - dDest[globalPosY * outputL + globalPosX] = minLine; //minLine - dDest[globalPosY * outputL + globalPosX + 1 * outputL * outputS] = maxLine; //maxLine - dDest[globalPosY * outputL + globalPosX + 2 * outputL * outputS] = (maxLine - minLine); //usedLines + dDest[globalPosY * outputL + globalPosX] = (maxLine - minLine); //usedLines + dDest[globalPosY * outputL + globalPosX + 1 * outputL * outputS] = minLine; //minLine + dDest[globalPosY * outputL + globalPosX + 2 * outputL * outputS] = maxLine; //maxLine } } \ No newline at end of file diff --git a/Modules/PhotoacousticsAlgorithms/files.cmake b/Modules/PhotoacousticsAlgorithms/files.cmake index f01c928c35..a20bd86761 100644 --- a/Modules/PhotoacousticsAlgorithms/files.cmake +++ b/Modules/PhotoacousticsAlgorithms/files.cmake @@ -1,21 +1,24 @@ set(CPP_FILES mitkPhotoacousticImage.cpp Algorithms/mitkPhotoacousticBeamformingFilter.cpp - Algorithms/OCL/mitkPhotoacousticOCLBeamformer.cpp - - Algorithms/OCL/mitkPhotoacousticBModeFilter.cpp - - Algorithms/OCL/mitkPhotoacousticOCLUsedLinesCalculation.cpp + Algorithms/OpenCLFilter/mitkPhotoacousticOCLBeamformer.cpp + Algorithms/OpenCLFilter/mitkPhotoacousticBModeFilter.cpp + Algorithms/OpenCLFilter/mitkPhotoacousticOCLUsedLinesCalculation.cpp + Algorithms/OpenCLFilter/mitkPhotoacousticOCLDelayCalculation.cpp + Algorithms/OpenCLFilter/mitkPhotoacousticOCLMemoryLocSum.cpp ) set(RESOURCE_FILES DASQuadratic.cl DMASQuadratic.cl DASSpherical.cl DMASSpherical.cl BModeAbs.cl BModeAbsLog.cl UsedLinesCalculation.cl + MemoryLocSum.cl + DelayCalculation.cl + DMAS.cl ) \ No newline at end of file diff --git a/Modules/PhotoacousticsAlgorithms/mitkPhotoacousticImage.cpp b/Modules/PhotoacousticsAlgorithms/mitkPhotoacousticImage.cpp index 0138062f02..7c50b25791 100644 --- a/Modules/PhotoacousticsAlgorithms/mitkPhotoacousticImage.cpp +++ b/Modules/PhotoacousticsAlgorithms/mitkPhotoacousticImage.cpp @@ -1,520 +1,520 @@ /*=================================================================== 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 "mitkPhotoacousticImage.h" #include "itkBModeImageFilter.h" #include "itkPhotoacousticBModeImageFilter.h" #include "mitkImageCast.h" #include "mitkITKImageImport.h" #include "mitkPhotoacousticBeamformingFilter.h" #include #include -#include +#include "./OpenCLFilter/mitkPhotoacousticBModeFilter.h" // 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 "itkFFT1DComplexConjugateToRealImageFilter.h" #include "itkFFT1DRealToComplexConjugateImageFilter.h" mitk::PhotoacousticImage::PhotoacousticImage() { MITK_INFO << "[PhotoacousticImage Debug] created that image"; } mitk::PhotoacousticImage::~PhotoacousticImage() { MITK_INFO << "[PhotoacousticImage Debug] destroyed that image"; } mitk::Image::Pointer mitk::PhotoacousticImage::ApplyBmodeFilter(mitk::Image::Pointer inputImage, BModeMethod method, bool UseGPU, bool UseLogFilter, float resampleSpacing) { // 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; typedef itk::IdentityTransform TransformType; if (method == BModeMethod::Abs) { mitk::Image::Pointer input; mitk::Image::Pointer out; if (inputImage->GetPixelType().GetTypeAsString() == "scalar (float)" || inputImage->GetPixelType().GetTypeAsString() == " (float)") input = inputImage; else input = ApplyCropping(inputImage, 0, 0, 0, 0, 0, inputImage->GetDimension(2) - 1); if (!UseGPU) { PhotoacousticBModeFilter::Pointer filter = PhotoacousticBModeFilter::New(); filter->SetParameters(UseLogFilter); filter->SetInput(input); filter->Update(); out = filter->GetOutput(); if (resampleSpacing == 0) return out; } #ifdef PHOTOACOUSTICS_USE_GPU else { PhotoacousticOCLBModeFilter::Pointer filter = PhotoacousticOCLBModeFilter::New(); filter->SetParameters(UseLogFilter); filter->SetInput(input); filter->Update(); out = filter->GetOutput(); if (resampleSpacing == 0) return out; } #endif typedef itk::ResampleImageFilter < itkFloatImageType, itkFloatImageType > ResampleImageFilter; ResampleImageFilter::Pointer resampleImageFilter = ResampleImageFilter::New(); itkFloatImageType::Pointer itkImage; 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]; outputSize[1] = inputSize[1] * itkImage->GetSpacing()[1] / outputSpacing[1]; typedef itk::IdentityTransform TransformType; resampleImageFilter->SetInput(itkImage); resampleImageFilter->SetSize(outputSize); resampleImageFilter->SetOutputSpacing(outputSpacing); resampleImageFilter->SetTransform(TransformType::New()); resampleImageFilter->UpdateLargestPossibleRegion(); return mitk::GrabItkImageMemory(resampleImageFilter->GetOutput()); } else if (method == BModeMethod::ShapeDetection) { 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; 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(); } // resampleSpacing == 0 means: do no resampling if (resampleSpacing == 0) { return mitk::GrabItkImageMemory(bmode); } 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]; outputSize[1] = inputSize[1] * itkImage->GetSpacing()[1] / outputSpacing[1]; resampleImageFilter->SetInput(bmode); resampleImageFilter->SetSize(outputSize); resampleImageFilter->SetOutputSpacing(outputSpacing); resampleImageFilter->SetTransform(TransformType::New()); resampleImageFilter->UpdateLargestPossibleRegion(); return mitk::GrabItkImageMemory(resampleImageFilter->GetOutput()); } return nullptr; } /*mitk::Image::Pointer mitk::PhotoacousticImage::ApplyScatteringCompensation(mitk::Image::Pointer inputImage, int scattering) { typedef itk::Image< float, 3 > itkFloatImageType; typedef itk::MultiplyImageFilter MultiplyImageFilterType; itkFloatImageType::Pointer itkImage; mitk::CastToItkImage(inputImage, itkImage); MultiplyImageFilterType::Pointer multiplyFilter = MultiplyImageFilterType::New(); multiplyFilter->SetInput1(itkImage); multiplyFilter->SetInput2(m_FluenceCompResizedItk.at(m_ScatteringCoefficient)); return mitk::GrabItkImageMemory(multiplyFilter->GetOutput()); }*/ mitk::Image::Pointer mitk::PhotoacousticImage::ApplyResampling(mitk::Image::Pointer inputImage, unsigned int outputSize[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[2] = itkImage->GetSpacing()[2]; typedef itk::IdentityTransform TransformType; T_Interpolator::Pointer _pInterpolator = T_Interpolator::New(); 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::PhotoacousticImage::ApplyCropping(mitk::Image::Pointer inputImage, int above, int below, int right, int left, int minSlice, int maxSlice) { unsigned int inputDim[3] = { inputImage->GetDimension(0), inputImage->GetDimension(1), inputImage->GetDimension(2) }; unsigned int outputDim[3] = { inputImage->GetDimension(0) - left - right, inputImage->GetDimension(1) - (unsigned int)above - (unsigned int)below, (unsigned int)maxSlice - (unsigned int)minSlice + 1 }; void* inputData; float* outputData = new float[outputDim[0] * outputDim[1] * outputDim[2]]; ImageReadAccessor acc(inputImage); inputData = const_cast(acc.GetData()); // convert the data to float by default // as of now only those float, short, float are used at all... though it's easy to add other ones if (inputImage->GetPixelType().GetTypeAsString() == "scalar (float)" || inputImage->GetPixelType().GetTypeAsString() == " (float)") { // copy the data into the cropped image for (unsigned short sl = 0; sl < outputDim[2]; ++sl) { for (unsigned short l = 0; l < outputDim[0]; ++l) { for (unsigned short s = 0; s < outputDim[1]; ++s) { outputData[l + s*(unsigned short)outputDim[0] + sl*outputDim[0] * outputDim[1]] = (float)((float*)inputData)[(l + left) + (s + above)*(unsigned short)inputDim[0] + (sl + minSlice)*inputDim[0] * inputDim[1]]; } } } } else if (inputImage->GetPixelType().GetTypeAsString() == "scalar (short)" || inputImage->GetPixelType().GetTypeAsString() == " (short)") { // copy the data unsigned shorto the cropped image for (unsigned short sl = 0; sl < outputDim[2]; ++sl) { for (unsigned short l = 0; l < outputDim[0]; ++l) { for (unsigned short s = 0; s < outputDim[1]; ++s) { outputData[l + s*(unsigned short)outputDim[0] + sl*outputDim[0] * outputDim[1]] = (float)((short*)inputData)[(l + left) + (s + above)*(unsigned short)inputDim[0] + (sl + minSlice)*inputDim[0] * inputDim[1]]; } } } } else if (inputImage->GetPixelType().GetTypeAsString() == "scalar (double)" || inputImage->GetPixelType().GetTypeAsString() == " (double)") { // copy the data unsigned shorto the cropped image for (unsigned short sl = 0; sl < outputDim[2]; ++sl) { for (unsigned short l = 0; l < outputDim[0]; ++l) { for (unsigned short s = 0; s < outputDim[1]; ++s) { outputData[l + s*(unsigned short)outputDim[0] + sl*outputDim[0] * outputDim[1]] = (float)((double*)inputData)[(l + left) + (s + above)*(unsigned short)inputDim[0] + (sl + minSlice)*inputDim[0] * inputDim[1]]; } } } } else { MITK_INFO << "Could not determine pixel type"; } mitk::Image::Pointer output = mitk::Image::New(); output->Initialize(mitk::MakeScalarPixelType(), 3, outputDim); output->SetSpacing(inputImage->GetGeometry()->GetSpacing()); output->SetImportVolume(outputData, 0, 0, mitk::Image::ReferenceMemory); return output; } mitk::Image::Pointer mitk::PhotoacousticImage::ApplyBeamforming(mitk::Image::Pointer inputImage, BeamformingFilter::beamformingSettings config, int cutoff, std::function progressHandle) { config.RecordTime = config.RecordTime - (cutoff) / inputImage->GetDimension(1) * config.RecordTime; // adjust the recorded time lost by cropping progressHandle(0, "cropping image"); if (!config.partial) { config.CropBounds[0] = 0; config.CropBounds[1] = inputImage->GetDimension(2) - 1; } Image::Pointer processedImage = ApplyCropping(inputImage, cutoff, 0, 0, 0, config.CropBounds[0], config.CropBounds[1]); config.inputDim = processedImage->GetDimensions(); // perform the beamforming BeamformingFilter::Pointer Beamformer = BeamformingFilter::New(); Beamformer->SetInput(processedImage); Beamformer->Configure(config); Beamformer->SetProgressHandle(progressHandle); Beamformer->UpdateLargestPossibleRegion(); processedImage = Beamformer->GetOutput(); return processedImage; } mitk::Image::Pointer mitk::PhotoacousticImage::BandpassFilter(mitk::Image::Pointer data, float recordTime, float BPHighPass, float BPLowPass, float alpha) { 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)}; data = ApplyResampling(data, dim); } MITK_INFO << data->GetDimension(0); // do a fourier transform, multiply with an appropriate window for the filter, and transform back typedef float PixelType; typedef itk::Image< PixelType, 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; } 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, alpha); 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::PhotoacousticImage::BPFunction(mitk::Image::Pointer reference, int cutoffFrequencyPixelHighPass, int cutoffFrequencyPixelLowPass, float alpha) { float* imageData = new float[reference->GetDimension(0)*reference->GetDimension(1)]; // tukey window float width = reference->GetDimension(1) / 2 - (float)cutoffFrequencyPixelHighPass - (float)cutoffFrequencyPixelLowPass; float center = (float)cutoffFrequencyPixelHighPass / 2 + width / 2; MITK_INFO << width << "width " << center << "center " << alpha; for (unsigned int n = 0; n < reference->GetDimension(1); ++n) { imageData[reference->GetDimension(0)*n] = 0; } for (int n = 0; n < width; ++n) { if (n <= (alpha*(width - 1)) / 2) { imageData[reference->GetDimension(0)*(int)(n + center - (width / 2))] = (1 + cos(M_PI*(2 * n / (alpha*(width - 1)) - 1))) / 2; } else if (n >= (width - 1)*(1 - alpha / 2) && n <= (width - 1)) { imageData[reference->GetDimension(0)*(int)(n + center - (width / 2))] = (1 + cos(M_PI*(2 * n / (alpha*(width - 1)) + 1 - 2 / alpha))) / 2; } else { imageData[reference->GetDimension(0)*(int)(n + center - (width / 2))] = 1; } } // Butterworth-Filter /* // first, write the HighPass if (cutoffFrequencyPixelHighPass != reference->GetDimension(1) / 2) { for (int n = 0; n < reference->GetDimension(1) / 2; ++n) { imageData[reference->GetDimension(0)*n] = 1 / (1 + pow( (float)n / (float)(reference->GetDimension(1) / 2 - cutoffFrequencyPixelHighPass) , 2 * butterworthOrder)); } } else { for (int n = 0; n < reference->GetDimension(1) / 2; ++n) { imageData[reference->GetDimension(0)*n] = 1; } } // now, the LowPass for (int n = 0; n < reference->GetDimension(1) / 2; ++n) { imageData[reference->GetDimension(0)*n] *= 1 / (1 + pow( (float)(reference->GetDimension(1) / 2 - 1 - n) / (float)(reference->GetDimension(1) / 2 - cutoffFrequencyPixelLowPass) , 2 * butterworthOrder)); } */ // 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 (ImageType::IndexValueType slice = 0; slice < reference->GetDimension(2); ++slice) { for (ImageType::IndexValueType line = 0; line < reference->GetDimension(0); ++line) { for (ImageType::IndexValueType 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; } \ No newline at end of file