diff --git a/Modules/OpenCL/mitkOclDataSet.cpp b/Modules/OpenCL/mitkOclDataSet.cpp index 0d00e4938d..34d0fc87d3 100644 --- a/Modules/OpenCL/mitkOclDataSet.cpp +++ b/Modules/OpenCL/mitkOclDataSet.cpp @@ -1,155 +1,158 @@ /*=================================================================== 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 "mitkOclDataSet.h" #include "mitkCommon.h" #include "mitkLogMacros.h" #include "mitkOclUtils.h" #include mitk::OclDataSet::OclDataSet() : m_gpuBuffer(nullptr), m_context(nullptr), m_bufferSize(0), m_gpuModified(false), m_cpuModified(false), - m_Data(nullptr), m_size(0), m_BpE(1) + m_Data(nullptr), m_BpE(1) { } mitk::OclDataSet::~OclDataSet() { MITK_INFO << "OclDataSet Destructor"; //release GMEM Image buffer if (m_gpuBuffer) clReleaseMemObject(m_gpuBuffer); } -cl_mem mitk::OclDataSet::CreateGPUBuffer(unsigned int _size, unsigned int _bpe) +cl_mem mitk::OclDataSet::CreateGPUBuffer() { - MITK_INFO << "InitializeGPUBuffer call with: BPE=" << _bpe; - - m_bufferSize = _size; - - m_BpE = _bpe; + MITK_INFO << "InitializeGPUBuffer call with: BPE=" << m_BpE; us::ServiceReference ref = GetModuleContext()->GetServiceReference(); OclResourceService* resources = GetModuleContext()->GetService(ref); m_context = resources->GetContext(); int clErr; m_gpuBuffer = clCreateBuffer(m_context, CL_MEM_READ_WRITE, m_bufferSize * m_BpE, nullptr, &clErr); + MITK_INFO << "Created GPU Buffer Object of size: " << m_BpE* m_bufferSize << " Bytes"; + CHECK_OCL_ERR(clErr); return m_gpuBuffer; } bool mitk::OclDataSet::IsModified(int _type) { if (_type) return m_cpuModified; else return m_gpuModified; } void mitk::OclDataSet::Modified(int _type) { // defines... GPU: 0, CPU: 1 m_cpuModified = _type; m_gpuModified = !_type; } int mitk::OclDataSet::TransferDataToGPU(cl_command_queue gpuComQueue) { cl_int clErr = 0; // check whether an image present if (m_Data == nullptr){ MITK_ERROR("ocl.DataSet") << "(mitk) No data present!\n"; return -1; } // there is a need for copy only if RAM-Data newer then GMEM data if (m_cpuModified) { //check the buffer if(m_gpuBuffer == nullptr) { - MITK_ERROR("ocl.DataSet") << "(mitk) No GPU buffer present!\n"; - return -1; + CreateGPUBuffer(); } if (m_gpuBuffer != nullptr) { clErr = clEnqueueWriteBuffer(gpuComQueue, m_gpuBuffer, CL_TRUE, 0, m_bufferSize * m_BpE, m_Data, 0, NULL, NULL); + MITK_INFO << "Wrote Data to GPU Buffer Object."; } else { MITK_ERROR << "No GPU buffer present!"; } CHECK_OCL_ERR(clErr); m_gpuModified = true; } return clErr; } cl_mem mitk::OclDataSet::GetGPUBuffer() { // clGetMemObjectInfo() cl_mem_object_type memInfo; cl_int clErr = 0; // query image object info only if already initialized if( this->m_gpuBuffer ) { clErr = clGetMemObjectInfo(this->m_gpuBuffer, CL_MEM_TYPE, sizeof(cl_mem_object_type), &memInfo, nullptr ); CHECK_OCL_ERR(clErr); } MITK_INFO << "Querying info for object, recieving: " << memInfo; return m_gpuBuffer; } void* mitk::OclDataSet::TransferDataToCPU(cl_command_queue gpuComQueue) { cl_int clErr = 0; // if image created on GPU, needs to create mitk::Image if( m_gpuBuffer == nullptr ){ MITK_ERROR("ocl.DataSet") << "(mitk) No buffer present!\n"; return nullptr; } // check buffersize char* data = new char[m_bufferSize * m_BpE]; // debug info oclPrintMemObjectInfo( m_gpuBuffer ); clErr = clEnqueueReadBuffer( gpuComQueue, m_gpuBuffer, CL_FALSE, 0, m_bufferSize * m_BpE, data ,0, nullptr, nullptr); CHECK_OCL_ERR(clErr); clFlush( gpuComQueue ); // the cpu data is same as gpu this->m_gpuModified = false; return (void*) data; } -void mitk::OclDataSet::SetSize(unsigned short size) +void mitk::OclDataSet::SetBufferSize(unsigned int size) +{ + m_bufferSize = size; +} + +void mitk::OclDataSet::SetBpE(unsigned short BpE) { - m_size = size; + m_BpE = BpE; } \ No newline at end of file diff --git a/Modules/OpenCL/mitkOclDataSet.h b/Modules/OpenCL/mitkOclDataSet.h index 93926ace86..123e39a7a9 100644 --- a/Modules/OpenCL/mitkOclDataSet.h +++ b/Modules/OpenCL/mitkOclDataSet.h @@ -1,133 +1,135 @@ /*=================================================================== 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 __mitkOclDataSet_h #define __mitkOclDataSet_h #define GPU_DATA 0 #define CPU_DATA 1 #include #include "MitkOpenCLExports.h" #include "mitkOclBaseData.h" #include "mitkOpenCLActivator.h" #include namespace mitk { /*! * \brief Class implementing processing of arbitrary data sets for GPU Image Processing * * The class holds a pointer to the data stored in RAM and performs an * on-demand-copy to the graphics memory. It is the basic data structure for all * mitk::oclDataSetToDataSetFilter classes */ class MITKOPENCL_EXPORT OclDataSet : public OclBaseData { public: mitkClassMacro(OclDataSet, OclBaseData); itkFactorylessNewMacro(Self) itkCloneMacro(Self) /*! \brief Copies the RAM-stored data to GMEM */ virtual int TransferDataToGPU(cl_command_queue); /*! \brief Copies the in GMEM stored data to RAM */ virtual void* TransferDataToCPU(cl_command_queue); /*! \brief Returns the pointer to the referenced data */ void* GetData() { return m_Data; } /** Returns the pointer to the GPU buffer */ cl_mem GetGPUBuffer(); /** Create the GPU buffer for image * */ - cl_mem CreateGPUBuffer(unsigned int _size, unsigned int _bpe); + cl_mem CreateGPUBuffer(); /** \brief Returns the status of the image buffer * * @param _type The flag to specify the buffer type ( GPU / CPU ) */ bool IsModified(int _type); using OclBaseData::Modified; /** \brief Set the modified flag for one of the buffers * * @param _type The flag to specify the buffer type ( GPU / CPU ) */ void Modified(int _type); /** \brief Initialze the OclDataSet with data. */ void SetData(void* data) { this->m_cpuModified = true; this->m_gpuModified = false; m_Data = data; } - /*! \brief returns the size of the DataSet */ - int GetSize() const + /*! \brief returns the amount of elements in the DataSet */ + int GetBufferSize() const { - return this->m_size; + return this->m_bufferSize; } short GetBytesPerElement() const { return this->m_BpE; } - /** @brief Set the DataSet size*/ - void SetSize(unsigned short size); + /** @brief Set the amount of elements in buffer*/ + void SetBufferSize(unsigned int size); + + /** @brief Set the DataSet memory Size per Element in Bytes*/ + void SetBpE(unsigned short BpE); protected: /*! \brief Constructor */ OclDataSet(); /** @brief Destructor */ virtual ~OclDataSet(); /*! GMEM Image buffer */ cl_mem m_gpuBuffer; /*! GPU Context the Buffer was created in, needed for access */ cl_context m_context; - /*! GMEM Buffer Size */ - unsigned int m_bufferSize; - private: bool m_gpuModified; bool m_cpuModified; /*! Reference to the data */ void* m_Data; - unsigned short m_size; + /*! GMEM Buffer Size in elements*/ + unsigned int m_bufferSize; + /*! Bytes per Element in Buffer*/ unsigned short m_BpE; }; } #endif //__mitkOclDataSet_h diff --git a/Modules/OpenCL/mitkOclDataSetFilter.cpp b/Modules/OpenCL/mitkOclDataSetFilter.cpp index 83112f42ae..f30719b5ae 100644 --- a/Modules/OpenCL/mitkOclDataSetFilter.cpp +++ b/Modules/OpenCL/mitkOclDataSetFilter.cpp @@ -1,41 +1,55 @@ /*=================================================================== 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 "mitkOclDataSetFilter.h" #include "mitkOclFilter.h" #include "mitkOclDataSet.h" +#include "mitkImageReadAccessor.h" mitk::OclDataSetFilter::OclDataSetFilter() { } mitk::OclDataSetFilter::~OclDataSetFilter() { } void mitk::OclDataSetFilter::SetInput(mitk::OclDataSet::Pointer DataSet) { m_Input = DataSet; } void mitk::OclDataSetFilter::SetInput(void* DataSet, unsigned int size, unsigned int BpE) { m_Input = mitk::OclDataSet::New(); m_Input->SetData(DataSet); - m_Input->CreateGPUBuffer(size, BpE); - m_CurrentSize = size; + m_Input->SetBufferSize(size); + m_Input->SetBpE(BpE); + m_CurrentSize = BpE; } +void mitk::OclDataSetFilter::SetInput(mitk::Image::Pointer image) +{ + m_Input = mitk::OclDataSet::New(); + mitk::ImageReadAccessor reader(image); + + m_Input->SetData(const_cast(reader.GetData())); + m_CurrentSize = (unsigned int)(image->GetPixelType().GetBitsPerComponent() / 8); + unsigned int elements = image->GetDimension(0) * image->GetDimension(1) * image->GetDimension(2); + + m_Input->SetBufferSize(elements); + m_Input->SetBpE(m_CurrentSize); +} \ No newline at end of file diff --git a/Modules/OpenCL/mitkOclDataSetFilter.h b/Modules/OpenCL/mitkOclDataSetFilter.h index 1f09669a78..881291bbd0 100644 --- a/Modules/OpenCL/mitkOclDataSetFilter.h +++ b/Modules/OpenCL/mitkOclDataSetFilter.h @@ -1,63 +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 __mitkOclDataSetFilter_h #define __mitkOclDataSetFilter_h #include "mitkOclFilter.h" #include "mitkOclDataSet.h" #define FILTER_UCHAR 0 #define FILTER_SHORT 1 namespace mitk { class OclFilter; class OclDataSetFilter; /** * \brief The OclDataSetFilter is the topmost class for all filter which take DataSets as input. * * The input DataSet can be intialized via an oclDataSet or a pointer to the data * This makes it possible to create a filter pipeline of GPU-based filters * and to bind this part into the CPU (ITK) filter pipeline. */ class MITKOPENCL_EXPORT OclDataSetFilter: public OclFilter { public: /** * @brief SetInput SetInput Set the input DataSet (as mitk::OclDataSet). * @param DataSet The DataSet in mitk::OclDataSet. */ void SetInput(mitk::OclDataSet::Pointer DataSet); /** - * @brief SetInput Set the input DataSet (as mitk::DataSet). - * @param DataSet The DataSet in mitk::DataSet. + * @brief SetInput Set the input DataSet (as a pointer to the data). + * @param DataSet The DataSet in mitk::OclDataSet. */ void SetInput(void* DataSet, unsigned int size, unsigned int BpE); + /** + * @brief SetInput Set the input DataSet (as mitk::Image). + * @param DataSet The DataSet in mitk::OclDataSet. + */ + void SetInput(mitk::Image::Pointer image); + protected: OclDataSetFilter(); virtual ~OclDataSetFilter(); /** The input DataSet */ mitk::OclDataSet::Pointer m_Input; unsigned int m_CurrentSize; }; } #endif // __mitkOclDataSetFilter_h diff --git a/Modules/OpenCL/mitkOclDataSetToDataSetFilter.cpp b/Modules/OpenCL/mitkOclDataSetToDataSetFilter.cpp index bbf4faf052..c1a68cd7be 100644 --- a/Modules/OpenCL/mitkOclDataSetToDataSetFilter.cpp +++ b/Modules/OpenCL/mitkOclDataSetToDataSetFilter.cpp @@ -1,93 +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. ===================================================================*/ #include "mitkOclDataSetToDataSetFilter.h" #include "mitkOclDataSet.h" #include "mitkException.h" mitk::OclDataSetToDataSetFilter::OclDataSetToDataSetFilter() { m_Output = mitk::OclDataSet::New(); } mitk::OclDataSetToDataSetFilter::~OclDataSetToDataSetFilter() { } mitk::OclDataSet::Pointer mitk::OclDataSetToDataSetFilter::GetGPUOutput() { return this->m_Output; } void* mitk::OclDataSetToDataSetFilter::GetOutput() { void* pData = m_Output->TransferDataToCPU(m_CommandQue); return pData; } int mitk::OclDataSetToDataSetFilter::GetBytesPerElem() { return this->m_CurrentSizeOutput; } bool mitk::OclDataSetToDataSetFilter::InitExec(cl_kernel ckKernel, unsigned int* dimensions, unsigned int outputBpE) { cl_int clErr = 0; if (m_Input.IsNull()) mitkThrow() << "Input DataSet is null."; // get DataSet size once const unsigned int uiDataSetWidth = dimensions[0]; const unsigned int uiDataSetHeight = dimensions[1]; const unsigned int uiDataSetDepth = dimensions[2]; // compute work sizes this->SetWorkingSize(8, uiDataSetWidth, 8, uiDataSetHeight, 8, uiDataSetDepth); cl_mem clBuffIn = m_Input->GetGPUBuffer(); cl_mem clBuffOut = m_Output->GetGPUBuffer(); if (!clBuffIn) { if (m_Input->TransferDataToGPU(m_CommandQue) != CL_SUCCESS) { mitkThrow() << "Could not create / initialize gpu DataSet."; } clBuffIn = m_Input->GetGPUBuffer(); } // output DataSet not initialized if (!clBuffOut) { MITK_DEBUG << "Create GPU DataSet call " << uiDataSetWidth << "x" << uiDataSetHeight << "x" << uiDataSetDepth; - clBuffOut = m_Output->CreateGPUBuffer(uiDataSetWidth*uiDataSetHeight*uiDataSetDepth, outputBpE); + m_Output->SetBpE(outputBpE); + m_Output->SetBufferSize(uiDataSetWidth*uiDataSetHeight*uiDataSetDepth); + clBuffOut = m_Output->CreateGPUBuffer(); m_CurrentSizeOutput = outputBpE; } clErr = 0; clErr = clSetKernelArg(ckKernel, 0, sizeof(cl_mem), &clBuffIn); clErr |= clSetKernelArg(ckKernel, 1, sizeof(cl_mem), &clBuffOut); CHECK_OCL_ERR(clErr); if (clErr != CL_SUCCESS) mitkThrow() << "OpenCL Part initialization failed with " << GetOclErrorAsString(clErr); return(clErr == CL_SUCCESS); } \ No newline at end of file diff --git a/Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticOCLBeamformer.cpp b/Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticOCLBeamformer.cpp index 114810ddc9..e5d1c0684f 100644 --- a/Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticOCLBeamformer.cpp +++ b/Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticOCLBeamformer.cpp @@ -1,177 +1,204 @@ /*=================================================================== 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 "mitkPhotoacousticOCLBeamformer.h" #include "usServiceReference.h" mitk::PhotoacousticOCLBeamformer::PhotoacousticOCLBeamformer() -: m_PixelCalculation( NULL ) +: m_PixelCalculation( NULL ), m_InputImage(mitk::Image::New()) { this->AddSourceFile("DASQuadratic.cl"); this->AddSourceFile("DMASQuadratic.cl"); this->AddSourceFile("DASspherical.cl"); this->AddSourceFile("DMASspherical.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"; } mitk::PhotoacousticOCLBeamformer::~PhotoacousticOCLBeamformer() { if ( this->m_PixelCalculation ) { clReleaseKernel( m_PixelCalculation ); } } void mitk::PhotoacousticOCLBeamformer::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::PhotoacousticOCLBeamformer::Execute() { cl_int clErr = 0; try { - this->InitExec(this->m_PixelCalculation, m_OutputDim); + this->InitExec(this->m_PixelCalculation, m_OutputDim, 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); CHECK_OCL_ERR(clErr); // set kernel arguments clErr = clSetKernelArg( this->m_PixelCalculation, 2, sizeof(cl_mem), &cl_input); clErr |= clSetKernelArg( this->m_PixelCalculation, 3, sizeof(cl_ushort), &(this->m_ApodArraySize) ); clErr |= clSetKernelArg( this->m_PixelCalculation, 4, sizeof(cl_float), &(this->m_SpeedOfSound) ); clErr |= clSetKernelArg( this->m_PixelCalculation, 5, sizeof(cl_float), &(this->m_TimeSpacing) ); clErr |= clSetKernelArg( this->m_PixelCalculation, 6, sizeof(cl_float), &(this->m_Pitch) ); clErr |= clSetKernelArg( this->m_PixelCalculation, 7, sizeof(cl_float), &(this->m_Angle) ); clErr |= clSetKernelArg( this->m_PixelCalculation, 8, sizeof(cl_ushort), &(this->m_PAImage)); clErr |= clSetKernelArg(this->m_PixelCalculation, 9, sizeof(cl_ushort), &(this->m_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])); CHECK_OCL_ERR( clErr ); // execute the filter on a 3D NDRange - this->ExecuteKernel( m_PixelCalculation, 3); + if (m_OutputDim[2] == 1) + this->ExecuteKernel(m_PixelCalculation, 2); + else + this->ExecuteKernel(m_PixelCalculation, 3); // signalize the GPU-side data changed m_Output->Modified( GPU_DATA ); } us::Module *mitk::PhotoacousticOCLBeamformer::GetModule() { return us::GetModuleContext()->GetModule(); } bool mitk::PhotoacousticOCLBeamformer::Initialize() { bool buildErr = true; cl_int clErr = 0; if ( OclFilter::Initialize() ) { switch (m_Algorithm) { case BeamformingAlgorithm::DASQuad: { this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "ckDASQuad", &clErr); break; } case BeamformingAlgorithm::DMASQuad: { this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "ckDMASQuad", &clErr); break; } case BeamformingAlgorithm::DASSphe: { this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "ckDASSphe", &clErr); break; } case BeamformingAlgorithm::DMASSphe: { this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "ckDMASSphe", &clErr); break; } } buildErr |= CHECK_OCL_ERR( clErr ); } return (OclFilter::IsInitialized() && buildErr ); } void mitk::PhotoacousticOCLBeamformer::SetInput(mitk::Image::Pointer image) { - if(image->GetDimension() != 3) - { - mitkThrowException(mitk::Exception) << "Input for " << this->GetNameOfClass() << - " is not 3D. The filter only supports 3D. Please change your input."; - } - OclImageToImageFilter::SetInput(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); } -mitk::Image::Pointer mitk::PhotoacousticOCLBeamformer::GetOutput() +void mitk::PhotoacousticOCLBeamformer::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::PhotoacousticOCLBeamformer::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 = m_OutputDim; + 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_Input->GetMITKImage()->GetSlicedGeometry(); + const mitk::SlicedGeometry3D::Pointer p_slg = m_InputImage->GetSlicedGeometry(); MITK_DEBUG << "Creating new MITK Image."; - m_Output->GetMITKImage()->Initialize(this->GetOutputType(), dimension, dimensions); - m_Output->GetMITKImage()->SetSpacing(p_slg->GetSpacing()); - m_Output->GetMITKImage()->SetGeometry(m_Input->GetMITKImage()->GetGeometry()); - m_Output->GetMITKImage()->SetImportVolume(pData, 0, 0, mitk::Image::ReferenceMemory); + 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 m_Output->GetMITKImage(); + return outputImage; } diff --git a/Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticOCLBeamformer.h b/Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticOCLBeamformer.h index 9a7c90e0de..da7fc7455d 100644 --- a/Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticOCLBeamformer.h +++ b/Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticOCLBeamformer.h @@ -1,129 +1,141 @@ /*=================================================================== 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 "mitkOclImageToImageFilter.h" +#include "mitkOclDataSetToDataSetFilter.h" #include namespace mitk { class OclImageToImageFilter; /** 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 PhotoacousticOCLBeamformer : public OclImageToImageFilter, public itk::Object +class PhotoacousticOCLBeamformer : public OclDataSetToDataSetFilter, public itk::Object { public: mitkClassMacroItkParent(PhotoacousticOCLBeamformer, 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. */ void SetInput(Image::Pointer image); + void SetInput(void* data, unsigned int* dimensions, unsigned int BpE); - mitk::Image::Pointer GetOutput(); + mitk::Image::Pointer GetOutputAsImage(); /** Update the filter */ void Update(); void SetOutputDim( unsigned int outputDim[3]) { m_OutputDim[0] = outputDim[0]; m_OutputDim[1] = outputDim[1]; m_OutputDim[2] = outputDim[2]; } void SetApodisation(float* apodisation, unsigned short apodArraySize) { m_ApodArraySize = apodArraySize; m_Apodisation = apodisation; } enum BeamformingAlgorithm { DASQuad, DMASQuad, DASSphe, DMASSphe }; void SetAlgorithm(BeamformingAlgorithm algorithm, bool PA) { m_Algorithm = algorithm; m_PAImage = PA; } void SetBeamformingParameters(float SpeedOfSound, float timeSpacing, float Pitch, float Angle, bool PAImage, unsigned short transducerElements) { m_SpeedOfSound = SpeedOfSound; m_TimeSpacing = timeSpacing; m_Pitch = Pitch; m_Angle = Angle; m_PAImage = PAImage; m_TransducerElements = transducerElements; } protected: /** Constructor */ PhotoacousticOCLBeamformer(); /** Destructor */ virtual ~PhotoacousticOCLBeamformer(); /** 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; BeamformingAlgorithm m_Algorithm; unsigned short m_PAImage; float m_SpeedOfSound; float m_TimeSpacing; float m_Pitch; float m_Angle; unsigned short m_TransducerElements; + + mitk::Image::Pointer m_InputImage; }; } -#endif \ No newline at end of file +#endif + +/* +float* data = new float[m_InputDim[0] * m_InputDim[1] * m_InputDim[2]]; +for (unsigned int i = 0; i < m_InputDim[0] * m_InputDim[1] * m_InputDim[2]; ++i) +{ + data[i] = i; +} +mitk::OclDataSetToDataSetFilter::SetInput(data, m_InputDim[0] * m_InputDim[1] * m_InputDim[2], sizeof(float));*/ \ No newline at end of file diff --git a/Modules/PhotoacousticsAlgorithms/Algorithms/mitkPhotoacousticBeamformingFilter.cpp b/Modules/PhotoacousticsAlgorithms/Algorithms/mitkPhotoacousticBeamformingFilter.cpp index ca03f169a8..c8acae5b99 100644 --- a/Modules/PhotoacousticsAlgorithms/Algorithms/mitkPhotoacousticBeamformingFilter.cpp +++ b/Modules/PhotoacousticsAlgorithms/Algorithms/mitkPhotoacousticBeamformingFilter.cpp @@ -1,613 +1,603 @@ /*=================================================================== 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 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; float inputDim[2] = { (float)input->GetDimension(0), (float)input->GetDimension(1) }; float outputDim[2] = { (float)output->GetDimension(0), (float)output->GetDimension(1) }; unsigned short chunkSize = 2; // TODO: make this slightly less arbitrary unsigned int oclOutputDim[3] = { output->GetDimension(0), output->GetDimension(1), output->GetDimension(2) }; unsigned int oclOutputDimChunk[3] = { output->GetDimension(0), output->GetDimension(1), chunkSize}; unsigned int oclInputDimChunk[3] = { input->GetDimension(0), input->GetDimension(1), chunkSize}; unsigned int oclOutputDimLastChunk[3] = { output->GetDimension(0), output->GetDimension(1), input->GetDimension(2) % chunkSize }; unsigned int oclInputDimLastChunk[3] = { input->GetDimension(0), input->GetDimension(1), input->GetDimension(2) % chunkSize }; const int apodArraySize = m_Conf.TransducerElements * 4; // 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; } int progInterval = output->GetDimension(2) / 20 > 1 ? output->GetDimension(2) / 20 : 1; // the interval at which we update the gui progress bar auto begin = std::chrono::high_resolution_clock::now(); // debbuging the performance... if (!m_Conf.UseGPU) { 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_OutputData = new float[m_Conf.ReconstructionLines*m_Conf.SamplesPerLine]; - m_InputDataPuffer = new float[input->GetDimension(0)*input->GetDimension(1)]; - - // first, we convert any data to float, which we use by default + // 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; - delete[] m_InputDataPuffer; m_OutputData = nullptr; m_InputData = nullptr; } } else { mitk::PhotoacousticOCLBeamformer::Pointer m_oclFilter = mitk::PhotoacousticOCLBeamformer::New(); try { if (m_Conf.Algorithm == beamformingSettings::BeamformingAlgorithm::DAS) { if (m_Conf.DelayCalculationMethod == beamformingSettings::DelayCalc::QuadApprox) m_oclFilter->SetAlgorithm(PhotoacousticOCLBeamformer::BeamformingAlgorithm::DASQuad, true); else if (m_Conf.DelayCalculationMethod == beamformingSettings::DelayCalc::Spherical) m_oclFilter->SetAlgorithm(PhotoacousticOCLBeamformer::BeamformingAlgorithm::DASSphe, true); } else if (m_Conf.Algorithm == beamformingSettings::BeamformingAlgorithm::DMAS) { if (m_Conf.DelayCalculationMethod == beamformingSettings::DelayCalc::QuadApprox) m_oclFilter->SetAlgorithm(PhotoacousticOCLBeamformer::BeamformingAlgorithm::DMASQuad, true); else if (m_Conf.DelayCalculationMethod == beamformingSettings::DelayCalc::Spherical) m_oclFilter->SetAlgorithm(PhotoacousticOCLBeamformer::BeamformingAlgorithm::DMASSphe, true); } m_oclFilter->SetApodisation(ApodWindow, apodArraySize); m_oclFilter->SetOutputDim(oclOutputDimChunk); m_oclFilter->SetBeamformingParameters(m_Conf.SpeedOfSound, m_Conf.TimeSpacing, m_Conf.Pitch, m_Conf.Angle, m_Conf.Photoacoustic, m_Conf.TransducerElements); + mitk::ImageReadAccessor inputReadAccessor(input); + + // first, we check whether the data 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; + } + if (chunkSize < oclOutputDim[2]) { - bool skip = false; - for (unsigned int i = 0; !skip && i < ceil((float)oclOutputDim[2] / (float)chunkSize); ++i) + unsigned short curChunkSize = chunkSize; + + for (unsigned int i = 0; i < ceil((float)oclOutputDim[2] / (float)chunkSize); ++i) { m_ProgressHandle(100 * ((float)(i * chunkSize) / (float)oclOutputDim[2]), "performing reconstruction"); - mitk::Image::Pointer chunk = mitk::Image::New(); - if ((int)((oclOutputDim[2]) - (i * chunkSize)) == (int)(1 + chunkSize)) + + if ((oclOutputDim[2]) - (i * chunkSize) >= chunkSize) { - // A 3d image of 3rd dimension == 1 can not be processed by openCL, make sure that this case never arises - oclInputDimLastChunk[2] = input->GetDimension(2) % chunkSize + chunkSize; - oclOutputDimLastChunk[2] = input->GetDimension(2) % chunkSize + chunkSize; - - chunk->Initialize(input->GetPixelType(), 3, oclInputDimLastChunk); - m_oclFilter->SetOutputDim(oclOutputDimLastChunk); - skip = true; //skip the last chunk + m_oclFilter->SetInput((void*)&m_InputData[i * chunkSize * input->GetDimension(0) * input->GetDimension(1)], oclInputDimChunk, sizeof(float)); } - else if ((oclOutputDim[2]) - (i * chunkSize) >= chunkSize) - chunk->Initialize(input->GetPixelType(), 3, oclInputDimChunk); else { - chunk->Initialize(input->GetPixelType(), 3, oclInputDimLastChunk); - m_oclFilter->SetOutputDim(oclOutputDimLastChunk); + m_oclFilter->SetInput((void*)&m_InputData[i * chunkSize * input->GetDimension(0) * input->GetDimension(1)], oclInputDimLastChunk, sizeof(float)); + curChunkSize = (unsigned short)oclOutputDimLastChunk[2]; // the last chunk might have a different size! } - - chunk->SetSpacing(input->GetGeometry()->GetSpacing()); - - mitk::ImageReadAccessor ChunkMove(input); - chunk->SetImportVolume((void*)&(((float*)const_cast(ChunkMove.GetData()))[i * chunkSize * input->GetDimension(0) * input->GetDimension(1)]), 0, 0, mitk::Image::ReferenceMemory); - - m_oclFilter->SetInput(chunk); - m_oclFilter->Update(); - auto out = m_oclFilter->GetOutput(); - unsigned int lastCopiedSlice = 0; - if(skip) - { - lastCopiedSlice = i * chunkSize + 1; - } + m_oclFilter->Update(); + float* out = (float*)m_oclFilter->GetOutput(); - for (unsigned int s = i * chunkSize; s < oclOutputDim[2]; ++s) // TODO: make the bounds here smaller... + for (unsigned short s = 0; s < curChunkSize; ++s) { - mitk::ImageReadAccessor copy(out, out->GetSliceData(s - i * chunkSize)); - output->SetImportSlice(const_cast(copy.GetData()), s, 0, 0, mitk::Image::ReferenceMemory); + output->SetImportSlice( (void*)&(out[s * oclOutputDim[0] * oclOutputDim[1]]), i * chunkSize + s, 0, 0, mitk::Image::ReferenceMemory); } } } else { m_ProgressHandle(50, "performing reconstruction"); m_oclFilter->SetOutputDim(oclOutputDim); m_oclFilter->SetInput(input); m_oclFilter->Update(); - auto out = m_oclFilter->GetOutput(); - mitk::ImageReadAccessor copy(out); - output->SetImportVolume(const_cast(copy.GetData()), 0, 0, mitk::Image::ReferenceMemory); + float* out = (float*)m_oclFilter->GetOutput(); + output->SetImportVolume((void*)out, 0, 0, mitk::Image::ReferenceMemory); } } catch (mitk::Exception &e) { std::string errorMessage = "Caught unexpected exception "; errorMessage.append(e.what()); MITK_ERROR << errorMessage; } } 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 * m_Conf.ReconstructionLines / 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); 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 * m_Conf.ReconstructionLines / 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); 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 * m_Conf.ReconstructionLines / 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); 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 * m_Conf.ReconstructionLines / 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); //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 8afda8ff42..b97208b5bb 100644 --- a/Modules/PhotoacousticsAlgorithms/Resources/DASQuadratic.cl +++ b/Modules/PhotoacousticsAlgorithms/Resources/DASQuadratic.cl @@ -1,76 +1,77 @@ /*=================================================================== 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( - __read_only image3d_t dSource, // input image + __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 // parameters + unsigned short TransducerElements, + unsigned int inputL, + unsigned int inputS, + unsigned int Slices // 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); unsigned short outputS = get_global_size(1); unsigned short outputL = get_global_size(0); - - // get image width and weight - const unsigned int inputL = get_image_width( dSource ); - const unsigned int inputS = get_image_height( dSource ); - const unsigned int Slices = get_image_depth( dSource ); // create an image sampler const sampler_t defaultSampler = CLK_NORMALIZED_COORDS_FALSE | CLK_ADDRESS_CLAMP_TO_EDGE | CLK_FILTER_NEAREST ; // 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) * TimeSpacing * SpeedOfSound / Pitch * outputL / TransducerElements) * s_i; 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); - - short AddSample = 0; - float output = 0; - float delayMultiplicator = pow(1 / (TimeSpacing*SpeedOfSound) * Pitch * TransducerElements / inputL, 2) / s_i / 2; + float apod_mult = apodArraySize / (maxLine - minLine); + + 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)] * read_imagef( dSource, defaultSampler, (int4)(l_s, AddSample, globalPosZ, 0 )).x; + 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/DASspherical.cl b/Modules/PhotoacousticsAlgorithms/Resources/DASspherical.cl index f9a11f7504..ab1be8f62a 100644 --- a/Modules/PhotoacousticsAlgorithms/Resources/DASspherical.cl +++ b/Modules/PhotoacousticsAlgorithms/Resources/DASspherical.cl @@ -1,79 +1,77 @@ /*=================================================================== 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( - __read_only image3d_t dSource, // input image + __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 // parameters + unsigned short TransducerElements, + unsigned int inputL, + unsigned int inputS, + unsigned int Slices // 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); unsigned short outputS = get_global_size(1); unsigned short outputL = get_global_size(0); - - // get image width and weight - const unsigned int inputL = get_image_width( dSource ); - const unsigned int inputS = get_image_height( dSource ); - const unsigned int Slices = get_image_depth( dSource ); // create an image sampler const sampler_t defaultSampler = CLK_NORMALIZED_COORDS_FALSE | CLK_ADDRESS_CLAMP_TO_EDGE | CLK_FILTER_NEAREST ; // 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) * TimeSpacing * SpeedOfSound / Pitch * outputL / TransducerElements) * s_i; 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); - - short AddSample = 0; - float output = 0; + float apod_mult = apodArraySize / (maxLine - minLine); + + 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)] * read_imagef( dSource, defaultSampler, (int4)(l_s, AddSample, globalPosZ, 0 )).x; + 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/DMASQuadratic.cl b/Modules/PhotoacousticsAlgorithms/Resources/DMASQuadratic.cl index 857b988746..e8ed09ae2a 100644 --- a/Modules/PhotoacousticsAlgorithms/Resources/DMASQuadratic.cl +++ b/Modules/PhotoacousticsAlgorithms/Resources/DMASQuadratic.cl @@ -1,92 +1,90 @@ /*=================================================================== 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 ckDMASQuad( - __read_only image3d_t dSource, // input image + __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 // parameters + unsigned short TransducerElements, + unsigned int inputL, + unsigned int inputS, + unsigned int Slices // 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); unsigned short outputS = get_global_size(1); unsigned short outputL = get_global_size(0); - - // get image width and weight - const unsigned int inputL = get_image_width( dSource ); - const unsigned int inputS = get_image_height( dSource ); - const unsigned int Slices = get_image_depth( dSource ); // create an image sampler const sampler_t defaultSampler = CLK_NORMALIZED_COORDS_FALSE | CLK_ADDRESS_CLAMP_TO_EDGE | CLK_FILTER_NEAREST ; // 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) * TimeSpacing * SpeedOfSound / Pitch * outputL / TransducerElements) * s_i; 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 = apodArraySize / (maxLine - minLine); float delayMultiplicator = pow((1 / (TimeSpacing*SpeedOfSound) * Pitch * TransducerElements / inputL), 2) / s_i / 2; - float mult = 0; - float output = 0; - float AddSample1 = 0; - float AddSample2 = 0; + float mult = 0; + float output = 0; + float AddSample1 = 0; + float AddSample2 = 0; for (short l_s1 = minLine; l_s1 < maxLine - 1; ++l_s1) { AddSample1 = delayMultiplicator * pow((l_s1 - l_i), 2) + s_i + (1-PAImage)*s_i; if (AddSample1 < inputS && AddSample1 >= 0) { for (short l_s2 = l_s1 + 1; l_s2 < maxLine; ++l_s2) { AddSample2 = delayMultiplicator * pow((l_s2 - l_i), 2) + s_i + (1-PAImage)*s_i; if (AddSample2 < inputS && AddSample2 >= 0) { - mult = read_imagef( dSource, defaultSampler, (int4)(l_s2, AddSample2, globalPosZ, 0 )).x - * apodArray[(short)((l_s2 - minLine)*apod_mult)] - * read_imagef( dSource, defaultSampler, (int4)(l_s1, AddSample1, globalPosZ, 0 )).x - * apodArray[(short)((l_s1 - minLine)*apod_mult)]; + mult = dSource[(unsigned int)(globalPosZ * inputL * inputS + AddSample2 * inputL + l_s2)] + * apodArray[(short)((l_s2 - minLine)*apod_mult)] + * dSource[(unsigned int)(globalPosZ * inputL * inputS + AddSample1 * inputL + l_s1)] + * apodArray[(short)((l_s1 - minLine)*apod_mult)]; output += sqrt(mult * ((float)(mult>0)-(float)(mult<0))) * ((mult > 0) - (mult < 0)); } } } else --usedLines; } dDest[ globalPosZ * outputL * outputS + globalPosY * outputL + globalPosX ] = 10 * output / (pow((float)usedLines, 2.0f) - (usedLines - 1)); } } \ No newline at end of file diff --git a/Modules/PhotoacousticsAlgorithms/Resources/DMASspherical.cl b/Modules/PhotoacousticsAlgorithms/Resources/DMASspherical.cl index fc480b0820..8d74488e57 100644 --- a/Modules/PhotoacousticsAlgorithms/Resources/DMASspherical.cl +++ b/Modules/PhotoacousticsAlgorithms/Resources/DMASspherical.cl @@ -1,98 +1,96 @@ /*=================================================================== 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 ckDMASSphe( - __read_only image3d_t dSource, // input image + __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 // parameters + unsigned short TransducerElements, + unsigned int inputL, + unsigned int inputS, + unsigned int Slices // 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); unsigned short outputS = get_global_size(1); unsigned short outputL = get_global_size(0); - - // get image width and weight - const unsigned int inputL = get_image_width( dSource ); - const unsigned int inputS = get_image_height( dSource ); - const unsigned int Slices = get_image_depth( dSource ); // create an image sampler const sampler_t defaultSampler = CLK_NORMALIZED_COORDS_FALSE | CLK_ADDRESS_CLAMP_TO_EDGE | CLK_FILTER_NEAREST ; // 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) * TimeSpacing * SpeedOfSound / Pitch * outputL / TransducerElements) * s_i; 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 mult = 0; - float output = 0; - float AddSample1 = 0; - float AddSample2 = 0; + float apod_mult = apodArraySize / (maxLine - minLine); + + float mult = 0; + float output = 0; + float AddSample1 = 0; + float AddSample2 = 0; for (short l_s1 = minLine; l_s1 < maxLine - 1; ++l_s1) { AddSample1 = sqrt( pow(s_i, 2) + pow((1 / (TimeSpacing*SpeedOfSound) * ((l_s1 - l_i)*Pitch*TransducerElements)/inputL), 2) ) + (1-PAImage)*s_i; if (AddSample1 < inputS && AddSample1 >= 0) { for (short l_s2 = l_s1 + 1; l_s2 < maxLine; ++l_s2) { AddSample2 = sqrt( pow(s_i, 2) + pow((1 / (TimeSpacing*SpeedOfSound) * ((l_s2 - l_i)*Pitch*TransducerElements)/inputL), 2) ) + (1-PAImage)*s_i; if (AddSample2 < inputS && AddSample2 >= 0) { - mult = read_imagef( dSource, defaultSampler, (int4)(l_s2, AddSample2, globalPosZ, 0 )).x - * apodArray[(short)((l_s2 - minLine)*apod_mult)] - * read_imagef( dSource, defaultSampler, (int4)(l_s1, AddSample1, globalPosZ, 0 )).x - * apodArray[(short)((l_s1 - minLine)*apod_mult)]; + mult = dSource[(unsigned int)(globalPosZ * inputL * inputS + AddSample2 * inputL + l_s2)] + * apodArray[(short)((l_s2 - minLine)*apod_mult)] + * dSource[(unsigned int)(globalPosZ * inputL * inputS + AddSample1 * inputL + l_s1)] + * apodArray[(short)((l_s1 - minLine)*apod_mult)]; output += sqrt(mult * ((float)(mult>0)-(float)(mult<0))) * ((mult > 0) - (mult < 0)); } } } else --usedLines; } dDest[ globalPosZ * outputL * outputS + globalPosY * outputL + globalPosX ] = 10 * output / (pow((float)usedLines, 2.0f) - (usedLines - 1)); } } \ No newline at end of file