diff --git a/Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticBModeFilter.cpp b/Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticBModeFilter.cpp index 7c32c388af..f1ed562d11 100644 --- a/Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticBModeFilter.cpp +++ b/Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticBModeFilter.cpp @@ -1,186 +1,218 @@ /*=================================================================== 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 "mitkPhotoacousticBModeFilter.h" #include "usServiceReference.h" #include #ifdef PHOTOACOUSTICS_USE_GPU mitk::PhotoacousticOCLBModeFilter::PhotoacousticOCLBModeFilter() : m_PixelCalculation(NULL) { this->AddSourceFile("BModeAbs.cl"); this->AddSourceFile("BModeAbsLog.cl"); this->m_FilterID = "PixelCalculation"; } mitk::PhotoacousticOCLBModeFilter::~PhotoacousticOCLBModeFilter() { if (this->m_PixelCalculation) { clReleaseKernel(m_PixelCalculation); } } void mitk::PhotoacousticOCLBModeFilter::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::PhotoacousticOCLBModeFilter::Execute() { try { - this->InitExec(this->m_PixelCalculation); + this->InitExec(this->m_PixelCalculation, m_InputDim, sizeof(float)); } catch (const mitk::Exception& e) { MITK_ERROR << "Catched exception while initializing filter: " << e.what(); return; } + cl_int clErr; + clErr = clSetKernelArg(this->m_PixelCalculation, 2, sizeof(cl_uint), &(this->m_Size)); + + 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::PhotoacousticOCLBModeFilter::GetModule() { return us::GetModuleContext()->GetModule(); } bool mitk::PhotoacousticOCLBModeFilter::Initialize() { bool buildErr = true; cl_int clErr = 0; if (OclFilter::Initialize()) { if(m_UseLogFilter) this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "ckBmodeAbsLog", &clErr); else this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "ckBmodeAbs", &clErr); buildErr |= CHECK_OCL_ERR(clErr); } return (OclFilter::IsInitialized() && buildErr); } void mitk::PhotoacousticOCLBModeFilter::SetInput(mitk::Image::Pointer image) { - if (image->GetDimension() != 3) + 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); + m_Size = m_InputDim[0] * m_InputDim[1] * m_InputDim[2]; +} + +mitk::Image::Pointer mitk::PhotoacousticOCLBModeFilter::GetOutput() +{ + mitk::Image::Pointer outputImage = mitk::Image::New(); + + if (m_Output->IsModified(GPU_DATA)) { - mitkThrowException(mitk::Exception) << "Input for " << this->GetNameOfClass() << - " is not 3D. The filter only supports 3D. Please change your input."; + void* pData = m_Output->TransferDataToCPU(m_CommandQue); + + const unsigned int dimension = 3; + unsigned int dimensions[3] = { m_InputDim[0], m_InputDim[1], m_InputDim[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); } - OclImageToImageFilter::SetInput(image); + + MITK_DEBUG << "Image Initialized."; + + return outputImage; } #endif mitk::PhotoacousticBModeFilter::PhotoacousticBModeFilter() : m_UseLogFilter(false) { this->SetNumberOfIndexedInputs(1); this->SetNumberOfRequiredInputs(1); } mitk::PhotoacousticBModeFilter::~PhotoacousticBModeFilter() { } void mitk::PhotoacousticBModeFilter::GenerateInputRequestedRegion() { Superclass::GenerateInputRequestedRegion(); mitk::Image* output = this->GetOutput(); mitk::Image* input = const_cast (this->GetInput()); if (!output->IsInitialized()) { return; } input->SetRequestedRegionToLargestPossibleRegion(); //GenerateTimeInInputRegion(output, input); } void mitk::PhotoacousticBModeFilter::GenerateOutputInformation() { mitk::Image::ConstPointer input = this->GetInput(); mitk::Image::Pointer output = this->GetOutput(); if ((output->IsInitialized()) && (this->GetMTime() <= m_TimeOfHeaderInitialization.GetMTime())) return; itkDebugMacro(<< "GenerateOutputInformation()"); output->Initialize(input); output->GetGeometry()->SetSpacing(input->GetGeometry()->GetSpacing()); output->GetGeometry()->Modified(); output->SetPropertyList(input->GetPropertyList()->Clone()); m_TimeOfHeaderInitialization.Modified(); } void mitk::PhotoacousticBModeFilter::GenerateData() { GenerateOutputInformation(); mitk::Image::Pointer input = this->GetInput(); mitk::Image::Pointer output = this->GetOutput(); if (!output->IsInitialized()) return; mitk::ImageReadAccessor reader(input); unsigned int size = output->GetDimension(0) * output->GetDimension(1) * output->GetDimension(2); float* InputData = (float*)const_cast(reader.GetData()); float* OutputData = new float[size]; if(!m_UseLogFilter) for (unsigned int i = 0; i < size; ++i) { OutputData[i] = abs(InputData[i]); } else { for (unsigned int i = 0; i < size; ++i) { OutputData[i] = log(abs(InputData[i])); } } output->SetImportVolume(OutputData, 0, 0, mitk::Image::ImportMemoryManagementType::ManageMemory); m_TimeOfHeaderInitialization.Modified(); } \ No newline at end of file diff --git a/Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticBModeFilter.h b/Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticBModeFilter.h index 00066cfff9..e37d9f23c9 100644 --- a/Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticBModeFilter.h +++ b/Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticBModeFilter.h @@ -1,129 +1,132 @@ /*=================================================================== 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 "mitkOclImageToImageFilter.h" #endif #include +#include "mitkOclDataSetToDataSetFilter.h" #include "mitkImageToImageFilter.h" namespace mitk { #ifdef PHOTOACOUSTICS_USE_GPU - class OclImageToImageFilter; /** Documentation * - * \brief The OclBinaryThresholdImageFilter computes a binary segmentation based on given - threshold values. - + * \brief The PhotoacousticOCLBModeFilter simply takes the absolute of all pixels in the image. * - * 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. + * The filter gives the option to use a log after taking the absolute. */ - class PhotoacousticOCLBModeFilter : public OclImageToImageFilter, public itk::Object + class PhotoacousticOCLBModeFilter : public OclDataSetToDataSetFilter, public itk::Object { public: mitkClassMacroItkParent(PhotoacousticOCLBModeFilter, 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. + * @brief SetInput Set the input image + * @param image a image. */ void SetInput(Image::Pointer image); /** Update the filter */ void Update(); void SetParameters(bool useLogFilter) { m_UseLogFilter = useLogFilter; - } + } + + /** + * @brief GetOutputAsImage Returns an mitk::Image constructed from the processed data + */ + mitk::Image::Pointer GetOutput(); protected: - /** Constructor */ PhotoacousticOCLBModeFilter(); - /** Destructor */ virtual ~PhotoacousticOCLBModeFilter(); - /** 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; bool m_UseLogFilter; + + mitk::Image::Pointer m_InputImage; + unsigned int m_InputDim[3]; + unsigned int m_Size; + }; #endif class PhotoacousticBModeFilter : public ImageToImageFilter { public: mitkClassMacro(PhotoacousticBModeFilter, ImageToImageFilter); itkFactorylessNewMacro(Self) itkCloneMacro(Self) void SetParameters(bool useLogFilter) { m_UseLogFilter = useLogFilter; } protected: PhotoacousticBModeFilter(); ~PhotoacousticBModeFilter(); virtual void GenerateInputRequestedRegion() override; virtual void GenerateOutputInformation() override; virtual void GenerateData() override; //##Description //## @brief Time when Header was last initialized itk::TimeStamp m_TimeOfHeaderInitialization; bool m_UseLogFilter; }; } #endif \ No newline at end of file diff --git a/Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticOCLBeamformer.h b/Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticOCLBeamformer.h index 8b77730c68..5d80900b6c 100644 --- a/Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticOCLBeamformer.h +++ b/Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticOCLBeamformer.h @@ -1,144 +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_ #ifdef PHOTOACOUSTICS_USE_GPU #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 OclDataSetToDataSetFilter, public itk::Object { public: mitkClassMacroItkParent(PhotoacousticOCLBeamformer, 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; } enum BeamformingAlgorithm { DASQuad, DMASQuad, DASSphe, DMASSphe }; /** Set which Algorithm should be used for beamforming */ void SetAlgorithm(BeamformingAlgorithm algorithm, bool PA) { m_Algorithm = algorithm; m_PAImage = PA; } /** Set various beamforming parameters */ 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 #endif \ No newline at end of file diff --git a/Modules/PhotoacousticsAlgorithms/Resources/BModeAbs.cl b/Modules/PhotoacousticsAlgorithms/Resources/BModeAbs.cl index ce692acbb6..5dfd66d02d 100644 --- a/Modules/PhotoacousticsAlgorithms/Resources/BModeAbs.cl +++ b/Modules/PhotoacousticsAlgorithms/Resources/BModeAbs.cl @@ -1,41 +1,38 @@ /*=================================================================== 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 ckBmodeAbs( - __read_only image3d_t dSource, // input image - __global float* dDest // output buffer + __global float* dSource, // input image + __global float* dDest, // output buffer + unsigned int size ) { // get thread identifier unsigned int globalPosX = get_global_id(0); unsigned int globalPosY = get_global_id(1); unsigned int globalPosZ = get_global_id(2); - // 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 ; + // get image width and height + unsigned short inputS = get_global_size(1); + unsigned short inputL = get_global_size(0); + unsigned short slices = get_global_size(2); // terminate non-valid threads - if ( globalPosX < inputL && globalPosY < inputS && globalPosZ < Slices ) + if ( globalPosX + inputL * globalPosY + inputL * inputS * globalPosZ < size ) { - - dDest[ globalPosZ * inputL * inputS + globalPosY * inputL + globalPosX ] = fabs(read_imagef( dSource, defaultSampler, (int4)(globalPosX, globalPosY, globalPosZ, 0 )).x); + dDest[ globalPosZ * inputL * inputS + globalPosY * inputL + globalPosX ] = fabs(dSource[ globalPosZ * inputL * inputS + globalPosY * inputL + globalPosX ]); } } \ No newline at end of file diff --git a/Modules/PhotoacousticsAlgorithms/Resources/BModeAbsLog.cl b/Modules/PhotoacousticsAlgorithms/Resources/BModeAbsLog.cl index 1daa1709e1..842e5a5b52 100644 --- a/Modules/PhotoacousticsAlgorithms/Resources/BModeAbsLog.cl +++ b/Modules/PhotoacousticsAlgorithms/Resources/BModeAbsLog.cl @@ -1,40 +1,38 @@ /*=================================================================== 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 ckBmodeAbsLog( - __read_only image3d_t dSource, // input image - __global float* dDest // output buffer + __global float* dSource, // input image + __global float* dDest, // output buffer + unsigned int size ) { // get thread identifier unsigned int globalPosX = get_global_id(0); unsigned int globalPosY = get_global_id(1); unsigned int globalPosZ = get_global_id(2); - // 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 ; + // get image width and height + unsigned short inputS = get_global_size(1); + unsigned short inputL = get_global_size(0); + unsigned short slices = get_global_size(2); // terminate non-valid threads - if ( globalPosX < inputL && globalPosY < inputS && globalPosZ < Slices ) + if ( globalPosX + inputL * globalPosY + inputL * inputS * globalPosZ < size ) { - dDest[ globalPosZ * inputL * inputS + globalPosY * inputL + globalPosX ] = log(fabs((float)read_imagef( dSource, defaultSampler, (int4)(globalPosX, globalPosY, globalPosZ, 0 )).x)); + dDest[ globalPosZ * inputL * inputS + globalPosY * inputL + globalPosX ] = log(fabs(dSource[ globalPosZ * inputL * inputS + globalPosY * inputL + globalPosX ])); } } \ No newline at end of file