diff --git a/CMake/BuildConfigurations/camiPhotoacousticsWorkstation.cmake b/CMake/BuildConfigurations/camiPhotoacousticsWorkstation.cmake new file mode 100644 index 0000000000..14a6e49fd8 --- /dev/null +++ b/CMake/BuildConfigurations/camiPhotoacousticsWorkstation.cmake @@ -0,0 +1,21 @@ +message(STATUS "Configuring MITK Photoacoustics Setup Build") + +# Enable open cv and open igt link, which is a necessary configuration +set(MITK_USE_OpenCV ON CACHE BOOL "MITK Use OpenCV Library" FORCE) +set(MITK_USE_OpenIGTLink ON CACHE BOOL "MITK Use OpenIGTLink Library" FORCE) + +set(MITK_USE_OpenCL ON CACHE BOOL "MITK Use OpenCL Library" FORCE) + +# Enable default plugins and the navigation modules +set(MITK_CONFIG_PLUGINS + org.mitk.gui.qt.datamanager + org.mitk.gui.qt.stdmultiwidgeteditor + org.mitk.gui.qt.imagenavigator + org.mitk.gui.qt.properties + org.mitk.gui.qt.viewnavigator + org.mitk.gui.qt.ultrasound + org.mitk.gui.qt.photoacoustics.imageprocessing + org.mitk.gui.qt.measurementtoolbox + org.mitk.gui.qt.pointsetinteraction +) + diff --git a/CMake/FindOpenCL.cmake b/CMake/FindOpenCL.cmake index e30fb7a484..848594e9e3 100644 --- a/CMake/FindOpenCL.cmake +++ b/CMake/FindOpenCL.cmake @@ -1,89 +1,88 @@ # - Try to find OpenCL # This module tries to find an OpenCL implementation on your system. It supports # AMD / ATI, Apple and NVIDIA implementations, but should work, too. # # To set manually the paths, define these environment variables: # OpenCL_INCPATH - Include path (e.g. OpenCL_INCPATH=/opt/cuda/4.0/cuda/include) # OpenCL_LIBPATH - Library path (e.h. OpenCL_LIBPATH=/usr/lib64/nvidia) # # Once done this will define # OPENCL_FOUND - system has OpenCL # OPENCL_INCLUDE_DIRS - the OpenCL include directory # OPENCL_LIBRARIES - link these to use OpenCL # # WIN32 should work, but is untested FIND_PACKAGE(PackageHandleStandardArgs) SET (OPENCL_VERSION_STRING "0.1.0") SET(OPENCL_VERSION_MAJOR 0) SET(OPENCL_VERSION_MINOR 1) SET(OPENCL_VERSION_PATCH 0) IF(APPLE) FIND_LIBRARY(OPENCL_LIBRARIES OpenCL DOC "OpenCL lib for OSX") FIND_PATH(OPENCL_INCLUDE_DIRS OpenCL/cl.h DOC "Include for OpenCL on OSX") FIND_PATH(_OPENCL_CPP_INCLUDE_DIRS OpenCL/cl.hpp DOC "Include for OpenCL CPP bindings on OSX") ELSE() IF (WIN32) FIND_PATH(OPENCL_INCLUDE_DIRS CL/cl.h) FIND_PATH(_OPENCL_CPP_INCLUDE_DIRS CL/cl.hpp) IF(CMAKE_SIZEOF_VOID_P MATCHES "8") SET(OPENCL_LIB_DIR "$ENV{ATISTREAMSDKROOT}/lib/x86_64") if(NOT IS_DIRECTORY ${OPENCL_LIB_DIR}) SET(OPENCL_LIB_DIR "$ENV{CUDA_PATH}/lib/x64") endif() ELSE() SET(OPENCL_LIB_DIR "$ENV{ATISTREAMSDKROOT}/lib/x86") if(NOT IS_DIRECTORY ${OPENCL_LIB_DIR}) # need to convert path in the cmake style ? SET(OPENCL_LIB_DIR "$ENV{CUDA_PATH}/lib/Win32") endif() ENDIF() file(TO_CMAKE_PATH ${OPENCL_LIB_DIR} OPENCL_LIB_DIR) GET_FILENAME_COMPONENT(OPENCL_LIB_DIR ${OPENCL_LIB_DIR} ABSOLUTE) FIND_LIBRARY(OPENCL_LIBRARIES OpenCL.lib PATHS ${OPENCL_LIB_DIR} ENV OpenCL_LIBPATH) GET_FILENAME_COMPONENT(_OPENCL_INC_CAND ${OPENCL_LIB_DIR}/../../include ABSOLUTE) # On Win32 search relative to the library FIND_PATH(OPENCL_INCLUDE_DIRS CL/cl.h PATHS "${_OPENCL_INC_CAND}" ENV OpenCL_INCPATH) FIND_PATH(_OPENCL_CPP_INCLUDE_DIRS CL/cl.hpp PATHS "${_OPENCL_INC_CAND}" ENV OpenCL_INCPATH) - ELSE (WIN32) + ELSE () # Unix style platforms - FIND_LIBRARY(OPENCL_LIBRARIES OpenCL - PATHS ENV LD_LIBRARY_PATH ENV OpenCL_LIBPATH + FIND_LIBRARY(OPENCL_LIBRARIES libOpenCL.so + PATHS "/usr/local/cuda/lib64" "/opt/AMDAPP/lib64" ENV LD_LIBRARY_PATH ENV OpenCL_LIBPATH ) - GET_FILENAME_COMPONENT(OPENCL_LIB_DIR ${OPENCL_LIBRARIES} PATH) GET_FILENAME_COMPONENT(_OPENCL_INC_CAND ${OPENCL_LIB_DIR}/../../include ABSOLUTE) # The AMD SDK currently does not place its headers # in /usr/include, therefore also search relative # to the library FIND_PATH(OPENCL_INCLUDE_DIRS CL/cl.h PATHS ${_OPENCL_INC_CAND} "/usr/local/cuda/include" "/opt/AMDAPP/include" ENV OpenCL_INCPATH) FIND_PATH(_OPENCL_CPP_INCLUDE_DIRS CL/cl.hpp PATHS ${_OPENCL_INC_CAND} "/usr/local/cuda/include" "/opt/AMDAPP/include" ENV OpenCL_INCPATH) - ENDIF (WIN32) + ENDIF () -ENDIF (APPLE) +ENDIF () FIND_PACKAGE_HANDLE_STANDARD_ARGS(OpenCL DEFAULT_MSG OPENCL_LIBRARIES OPENCL_INCLUDE_DIRS) IF(_OPENCL_CPP_INCLUDE_DIRS) SET( OPENCL_HAS_CPP_BINDINGS TRUE ) LIST( APPEND OPENCL_INCLUDE_DIRS ${_OPENCL_CPP_INCLUDE_DIRS} ) # This is often the same, so clean up LIST( REMOVE_DUPLICATES OPENCL_INCLUDE_DIRS ) ENDIF(_OPENCL_CPP_INCLUDE_DIRS) MARK_AS_ADVANCED( OPENCL_INCLUDE_DIRS ) diff --git a/Modules/OpenCL/files.cmake b/Modules/OpenCL/files.cmake index 799077854e..c1f158a404 100644 --- a/Modules/OpenCL/files.cmake +++ b/Modules/OpenCL/files.cmake @@ -1,23 +1,26 @@ set(CPP_FILES # helper classes mitkOclUtils.cpp mitkOclResourceServiceImpl_Private.cpp mitkOclImageFormats.cpp # module activator mitkOpenCLActivator.cpp # base data and filter objects mitkOclBaseData.cpp mitkOclImage.cpp + mitkOclDataSet.cpp mitkOclFilter.cpp mitkOclImageFilter.cpp + mitkOclDataSetFilter.cpp mitkOclImageToImageFilter.cpp + mitkOclDataSetToDataSetFilter.cpp # own filter implementations mitkOclBinaryThresholdImageFilter.cpp ) set(RESOURCE_FILES BinaryThresholdFilter.cl ) diff --git a/Modules/OpenCL/mitkOclDataSet.cpp b/Modules/OpenCL/mitkOclDataSet.cpp new file mode 100644 index 0000000000..514b2d100d --- /dev/null +++ b/Modules/OpenCL/mitkOclDataSet.cpp @@ -0,0 +1,160 @@ +/*=================================================================== + +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_BpE(1) +{ +} + +mitk::OclDataSet::~OclDataSet() +{ + MITK_INFO << "OclDataSet Destructor"; + + //release GMEM Image buffer + if (m_gpuBuffer) clReleaseMemObject(m_gpuBuffer); +} + + +cl_mem mitk::OclDataSet::CreateGPUBuffer() +{ + MITK_INFO << "InitializeGPUBuffer call with: BPE=" << m_BpE; + + us::ServiceReference ref = GetModuleContext()->GetServiceReference(); + OclResourceService* resources = GetModuleContext()->GetService(ref); + + m_context = resources->GetContext(); + + int clErr; + if (m_gpuBuffer) clReleaseMemObject(m_gpuBuffer); + + 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) + { + 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::SetBufferSize(unsigned int size) +{ + m_bufferSize = size; +} + +void mitk::OclDataSet::SetBpE(unsigned short BpE) +{ + m_BpE = BpE; +} \ No newline at end of file diff --git a/Modules/OpenCL/mitkOclDataSet.h b/Modules/OpenCL/mitkOclDataSet.h new file mode 100644 index 0000000000..1bc5dcadde --- /dev/null +++ b/Modules/OpenCL/mitkOclDataSet.h @@ -0,0 +1,134 @@ +/*=================================================================== + +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(); + + /** \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 amount of elements in the DataSet */ + int GetBufferSize() const + { + return this->m_bufferSize; + } + + short GetBytesPerElement() const + { + return this->m_BpE; + } + + /** @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; + +private: + /*! GMEM Buffer Size in elements*/ + unsigned int m_bufferSize; + + bool m_gpuModified; + bool m_cpuModified; + + /*! Reference to the data */ + void* m_Data; + + /*! Bytes per Element in Buffer*/ + unsigned short m_BpE; +}; + +} +#endif //__mitkOclDataSet_h diff --git a/Modules/OpenCL/mitkOclDataSetFilter.cpp b/Modules/OpenCL/mitkOclDataSetFilter.cpp new file mode 100644 index 0000000000..f30719b5ae --- /dev/null +++ b/Modules/OpenCL/mitkOclDataSetFilter.cpp @@ -0,0 +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->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 new file mode 100644 index 0000000000..881291bbd0 --- /dev/null +++ b/Modules/OpenCL/mitkOclDataSetFilter.h @@ -0,0 +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 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 new file mode 100644 index 0000000000..852f487a40 --- /dev/null +++ b/Modules/OpenCL/mitkOclDataSetToDataSetFilter.cpp @@ -0,0 +1,128 @@ +/*=================================================================== + +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, size_t outputDataSize, 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 or output buffer size changed + if (!clBuffOut || (size_t)m_Output->GetBufferSize() != outputDataSize) + { + MITK_DEBUG << "Create GPU DataSet call " << uiDataSetWidth << "x" << uiDataSetHeight << "x" << uiDataSetDepth; + m_Output->SetBpE(outputBpE); + m_Output->SetBufferSize(outputDataSize); + 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); +} + +bool mitk::OclDataSetToDataSetFilter::InitExecNoInput(cl_kernel ckKernel, unsigned int* dimensions, size_t outputDataSize, unsigned int outputBpE) +{ + cl_int clErr = 0; + + // 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 clBuffOut = m_Output->GetGPUBuffer(); + + // output DataSet not initialized + if (!clBuffOut) + { + MITK_DEBUG << "Create GPU DataSet call " << uiDataSetWidth << "x" << uiDataSetHeight << "x" << uiDataSetDepth; + m_Output->SetBpE(outputBpE); + m_Output->SetBufferSize(outputDataSize); + clBuffOut = m_Output->CreateGPUBuffer(); + m_CurrentSizeOutput = outputBpE; + } + + clErr = clSetKernelArg(ckKernel, 0, 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/OpenCL/mitkOclDataSetToDataSetFilter.h b/Modules/OpenCL/mitkOclDataSetToDataSetFilter.h new file mode 100644 index 0000000000..0891e8b876 --- /dev/null +++ b/Modules/OpenCL/mitkOclDataSetToDataSetFilter.h @@ -0,0 +1,81 @@ +/*=================================================================== + +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 __mitkOclDataSetToDataSetFilter_h +#define __mitkOclDataSetToDataSetFilter_h + +#include "mitkOclDataSetFilter.h" + +namespace mitk +{ +class OclDataSetFilter; +class OclDataSetToDataSetFilter; + +/** @class OclDataSetToDataSetFilter + * @brief The OclDataSetToDataSetFilter is the base class for all OpenCL DataSet filter generating DataSets. + */ +class MITKOPENCL_EXPORT OclDataSetToDataSetFilter: public OclDataSetFilter +{ +public: + /*! + * \brief Returns an pointer to the filtered data. + */ + void* GetOutput(); + + /*! + * \brief Returns a pointer to the graphics memory. + * + * Use this method when executing two and more filters on the GPU for fast access. + * This method does not copy the data to RAM. It returns only a pointer. + */ + mitk::OclDataSet::Pointer GetGPUOutput(); + +protected: + /** + * @brief OclDataSetToDataSetFilter Default constructor. + */ + OclDataSetToDataSetFilter(); + + /** @brief Destructor */ + virtual ~OclDataSetToDataSetFilter(); + + /** Output DataSet */ + mitk::OclDataSet::Pointer m_Output; + + /** @brief (Virtual) method Update() to be implemented in derived classes. */ + virtual void Update() = 0; + + /** + * @brief InitExec Initialize the execution + * @param ckKernel The GPU kernel. + * @throws mitk::Exception if something goes wrong. + * @return True for success. + */ + bool InitExec(cl_kernel ckKernel, unsigned int* dimensions, size_t outputDataSize, unsigned int outputBpE); + + bool InitExecNoInput(cl_kernel ckKernel, unsigned int* dimensions, size_t outputDataSize, unsigned int outputBpE); + + /** @brief Get the memory size needed for each element */ + virtual int GetBytesPerElem(); + + unsigned int m_CurrentSizeOutput; + +private: + /** Block dimensions */ + unsigned int m_BlockDims[3]; +}; +} +#endif // __mitkOclDataSetToDataSetFilter_h diff --git a/Modules/OpenCL/mitkOclFilter.cpp b/Modules/OpenCL/mitkOclFilter.cpp index 5312e59117..d1b1c23aab 100644 --- a/Modules/OpenCL/mitkOclFilter.cpp +++ b/Modules/OpenCL/mitkOclFilter.cpp @@ -1,248 +1,284 @@ /*=================================================================== 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. ===================================================================*/ //Ocl #include "mitkOclFilter.h" #include "mitkOclUtils.h" #include "mitkOpenCLActivator.h" //Mitk #include #include //usService #include "usServiceReference.h" #include #include #include #include #include #include mitk::OclFilter::OclFilter() : m_ClCompilerFlags(""), m_ClProgram(nullptr), m_CommandQue(nullptr), m_FilterID("mitkOclFilter"), m_Preambel(" "), m_Initialized(false) { } mitk::OclFilter::OclFilter(const char* filename) : m_ClCompilerFlags(""), m_ClProgram(nullptr), m_CommandQue(nullptr), m_FilterID(filename), m_Preambel(" "), m_Initialized(false) { m_ClFiles.push_back(filename); } mitk::OclFilter::~OclFilter() { MITK_DEBUG << "OclFilter Destructor"; // release program if (m_ClProgram) { us::ServiceReference ref = GetModuleContext()->GetServiceReference(); OclResourceService* resources = GetModuleContext()->GetService(ref); // remove program from storage resources->RemoveProgram(m_FilterID); } } bool mitk::OclFilter::ExecuteKernel( cl_kernel kernel, unsigned int workSizeDim ) { cl_int clErr = 0; clErr = clEnqueueNDRangeKernel( this->m_CommandQue, kernel, workSizeDim, nullptr, this->m_GlobalWorkSize, m_LocalWorkSize, 0, nullptr, nullptr); CHECK_OCL_ERR( clErr ); return ( clErr == CL_SUCCESS ); } +bool mitk::OclFilter::ExecuteKernelChunks( cl_kernel kernel, unsigned int workSizeDim, size_t* chunksDim ) +{ + size_t offset[3] ={0, 0, 0}; + cl_int clErr = 0; + + if(workSizeDim == 2) + { + for(offset[0] = 0; offset[0] < m_GlobalWorkSize[0]; offset[0] += chunksDim[0]) + { + for(offset[1] = 0; offset[1] < m_GlobalWorkSize[1]; offset[1] += chunksDim[1]) + { + clErr |= clEnqueueNDRangeKernel( this->m_CommandQue, kernel, workSizeDim, + offset, chunksDim, m_LocalWorkSize, 0, nullptr, nullptr); + } + } + } + else if(workSizeDim == 3) + { + for(offset[0] = 0; offset[0] < m_GlobalWorkSize[0]; offset[0] += chunksDim[0]) + { + for(offset[1] = 0; offset[1] < m_GlobalWorkSize[1]; offset[1] += chunksDim[1]) + { + for(offset[2] = 0; offset[2] < m_GlobalWorkSize[2]; offset[2] += chunksDim[2]) + { + clErr |= clEnqueueNDRangeKernel( this->m_CommandQue, kernel, workSizeDim, + offset, chunksDim, m_LocalWorkSize, 0, nullptr, nullptr); + } + } + } + } + + CHECK_OCL_ERR(clErr); + + return ( clErr == CL_SUCCESS ); +} + bool mitk::OclFilter::Initialize() { us::ServiceReference ref = GetModuleContext()->GetServiceReference(); OclResourceService* resources = GetModuleContext()->GetService(ref); m_CommandQue = resources->GetCommandQueue(); cl_int clErr = 0; m_Initialized = CHECK_OCL_ERR(clErr); if ( m_ClFiles.empty()) { MITK_ERROR<<"No OpenCL Source FILE specified"; return false; } if (m_ClProgram == nullptr) { try { this->m_ClProgram = resources->GetProgram( this->m_FilterID ); } catch(const mitk::Exception& e) { - MITK_INFO << "Program not stored in resource manager, compiling."; + MITK_INFO << "Program not stored in resource manager, compiling. " << e; this->CompileSource(); } } return m_Initialized; } void mitk::OclFilter::LoadSourceFiles(CStringList &sourceCode, ClSizeList &sourceCodeSize) { for( CStringList::iterator it = m_ClFiles.begin(); it != m_ClFiles.end(); ++it ) { MITK_DEBUG << "Load file :" << *it; us::ModuleResource mdr = GetModule()->GetResource(*it); if( !mdr.IsValid() ) MITK_WARN << "Could not load resource: " << mdr.GetName() << " is invalid!"; us::ModuleResourceStream rss(mdr); // read resource file to a string std::istreambuf_iterator eos; std::string source(std::istreambuf_iterator(rss), eos); // add preambel and build up string to compile std::string src(m_Preambel); src.append("\n"); src.append(source); // allocate new char buffer char* tmp = new char[src.size() + 1]; strcpy(tmp,src.c_str()); // add source to list sourceCode.push_back((const char*)tmp); sourceCodeSize.push_back(src.size()); } } void mitk::OclFilter::CompileSource() { // helper variable int clErr = 0; CStringList sourceCode; ClSizeList sourceCodeSize; if (m_ClFiles.empty()) { MITK_ERROR("ocl.filter") << "No shader source file was set"; return; } //get a valid opencl context us::ServiceReference ref = GetModuleContext()->GetServiceReference(); OclResourceService* resources = GetModuleContext()->GetService(ref); cl_context gpuContext = resources->GetContext(); // load the program source from file LoadSourceFiles(sourceCode, sourceCodeSize); if ( !sourceCode.empty() ) { // create program from all files in the file list m_ClProgram = clCreateProgramWithSource(gpuContext, sourceCode.size(), &sourceCode[0], &sourceCodeSize[0], &clErr); CHECK_OCL_ERR(clErr); // build the source code MITK_DEBUG << "Building Program Source"; std::string compilerOptions = ""; compilerOptions.append(m_ClCompilerFlags); MITK_DEBUG("ocl.filter") << "cl compiler flags: " << compilerOptions.c_str(); clErr = clBuildProgram(m_ClProgram, 0, nullptr, compilerOptions.c_str(), nullptr, nullptr); CHECK_OCL_ERR(clErr); // if OpenCL Source build failed if (clErr != CL_SUCCESS) { MITK_ERROR("ocl.filter") << "Failed to build source"; oclLogBuildInfo(m_ClProgram, resources->GetCurrentDevice() ); oclLogBinary(m_ClProgram, resources->GetCurrentDevice() ); m_Initialized = false; } // store the succesfully build program into the program storage provided by the resource service resources->InsertProgram(m_ClProgram, m_FilterID, true); // free the char buffers with the source code for( CStringList::iterator it = sourceCode.begin(); it != sourceCode.end(); ++it ) { delete[] *it; } } else { MITK_ERROR("ocl.filter") << "Could not load from source"; m_Initialized = false; } } void mitk::OclFilter::SetWorkingSize(unsigned int locx, unsigned int dimx, unsigned int locy, unsigned int dimy, unsigned int locz, unsigned int dimz) { // set the local work size this->m_LocalWorkSize[0] = locx; this->m_LocalWorkSize[1] = locy; this->m_LocalWorkSize[2] = locz; this->m_GlobalWorkSize[0] = dimx; this->m_GlobalWorkSize[1] = dimy; this->m_GlobalWorkSize[2] = dimz; // estimate the global work size this->m_GlobalWorkSize[0] = iDivUp( dimx, this->m_LocalWorkSize[0]) * this->m_LocalWorkSize[0]; if ( dimy > 1) this->m_GlobalWorkSize[1] = iDivUp( dimy, this->m_LocalWorkSize[1]) * this->m_LocalWorkSize[1]; if( dimz > 1 ) this->m_GlobalWorkSize[2] = iDivUp( dimz, this->m_LocalWorkSize[2]) * this->m_LocalWorkSize[2]; } void mitk::OclFilter::SetSourcePreambel(const char* preambel) { this->m_Preambel = preambel; } void mitk::OclFilter::AddSourceFile(const char* filename) { m_ClFiles.push_back(filename); } void mitk::OclFilter::SetCompilerFlags(const char* flags) { m_ClCompilerFlags = flags; } bool mitk::OclFilter::IsInitialized() { return m_Initialized; } diff --git a/Modules/OpenCL/mitkOclFilter.h b/Modules/OpenCL/mitkOclFilter.h index e4b1bb67fb..ba3705a466 100644 --- a/Modules/OpenCL/mitkOclFilter.h +++ b/Modules/OpenCL/mitkOclFilter.h @@ -1,149 +1,153 @@ /*=================================================================== 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 __mitkOclFilter_h #define __mitkOclFilter_h #include "mitkOclUtils.h" #include "mitkCommon.h" #include #include #include namespace mitk { /** @class OclFilter @brief Superclass for all OpenCL based filter. This class takes care of loading and compiling the external GPU program code. */ class MITKOPENCL_EXPORT OclFilter { public: /** * @brief Add a source file from the resource files to the * OpenCL shader file list. Multiple files can be added to the list. * * @param name of the file in the resource system */ void AddSourceFile(const char* filename); /** * @brief Set specific compilerflags to compile the CL source. Default is set to nullptr; * example: "-cl-fast-relaxed-math -cl-mad-enable -cl-strict-aliasing" * * @param flags to the modulefolder that contains the gpuSource */ void SetCompilerFlags(const char* flags); /** * @brief Returns true if the initialization was successfull */ virtual bool IsInitialized(); /** @brief Destructor */ virtual ~OclFilter(); protected: typedef std::vector CStringList; typedef std::vector ClSizeList; /** @brief Constructor */ OclFilter(); /** @brief Constructor ( overloaded ) */ OclFilter(const char* filename); /** @brief String that contains the compiler flags */ const char* m_ClCompilerFlags; /** @brief The compiled OpenCL program */ cl_program m_ClProgram; /** @brief Command queue for the filter */ cl_command_queue m_CommandQue; /** @brief Unique ID of the filter, needs to be specified in the constructor of the derived class */ std::string m_FilterID; /*! @brief source preambel for e.g. #define commands to be inserted into the OpenCL source */ const char* m_Preambel; /** @brief List of sourcefiles that will be compiled for this filter.*/ CStringList m_ClFiles; /** @brief status of the filter */ bool m_Initialized; /** @brief The local work size fo the filter */ size_t m_LocalWorkSize[3]; /** @brief The global work size of the filter */ size_t m_GlobalWorkSize[3]; /** @brief Set the working size for the following OpenCL kernel call */ void SetWorkingSize(unsigned int locx, unsigned int dimx, unsigned int locy = 1, unsigned int dimy = 1, unsigned int locz = 1, unsigned int dimz = 1); /** @brief Execute the given kernel on the OpenCL Index-Space defined by the local and global work sizes */ bool ExecuteKernel( cl_kernel kernel, unsigned int workSizeDim ); + /** @brief Execute the given kernel on the OpenCL Index-Space defined by the local and global work sizes, but divide it into chunks of dimension chunksDim + */ + bool ExecuteKernelChunks( cl_kernel kernel, unsigned int workSizeDim, size_t* chunksDim ); + /** * \brief Initialize all necessary parts of the filter * * The Initialize() method creates the command queue and the m_clProgram. * The program is either compiled from the given source or taken from the * OclResourceManager if the program was compiled already. */ bool Initialize(); /** * @brief Compile the program source * * @param preambel e.g. defines for the shader code */ void CompileSource(); /** * @brief Add some source code on the beginning of the loaded source * * In this way, some preprocessor flags for the CL compiler can at the beginning of the filter * @param preambel Source preambel for e.g. #define commands to be inserted into the OpenCL source */ void SetSourcePreambel(const char* preambel); /** * @brief Get the Module of the filter. Needs to be implemented by every subclass. * The filter will load the OpenCL sourcefiles from this module context. */ virtual us::Module* GetModule() = 0; /** * @brief Helper functions that load sourcefiles from the module context in the Initialize function. * @param SourceCodeList holds the sourcecode for every file as string, the SourceCodeSizeList holst the * size of every file in bytes. */ void LoadSourceFiles(CStringList &SourceCodeList, ClSizeList &SourceCodeSizeList); }; } #endif // __mitkOclFilter_h diff --git a/Modules/OpenCL/mitkOclImageToImageFilter.cpp b/Modules/OpenCL/mitkOclImageToImageFilter.cpp index f78d6f9def..10f88b06a4 100644 --- a/Modules/OpenCL/mitkOclImageToImageFilter.cpp +++ b/Modules/OpenCL/mitkOclImageToImageFilter.cpp @@ -1,193 +1,193 @@ /*=================================================================== 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 "mitkOclImageToImageFilter.h" #include "mitkOclImage.h" #include "mitkException.h" mitk::OclImageToImageFilter::OclImageToImageFilter() { m_Output = mitk::OclImage::New(); } mitk::OclImageToImageFilter::~OclImageToImageFilter() { } mitk::OclImage::Pointer mitk::OclImageToImageFilter::GetGPUOutput() { // initialize some variables m_Output->SetPixelType(m_Input->GetPixelType()); // create new image, for passing the essential information to the output m_Output->InitializeMITKImage(); const unsigned int dimension = m_Input->GetDimension(); unsigned int* dimensions = m_Input->GetDimensions(); m_Output->SetDimensions( dimensions ); m_Output->SetDimension( (unsigned short)dimension ); m_Output->GetMITKImage()->Initialize( this->GetOutputType(), dimension, dimensions); const mitk::SlicedGeometry3D::Pointer p_slg = m_Input->GetMITKImage()->GetSlicedGeometry(0); m_Output->GetMITKImage()->SetSpacing( p_slg->GetSpacing() ); m_Output->GetMITKImage()->SetGeometry( m_Input->GetMITKImage()->GetGeometry() ); return this->m_Output; } mitk::Image::Pointer mitk::OclImageToImageFilter::GetOutput() { if (m_Output->IsModified(GPU_DATA)) { void* pData = m_Output->TransferDataToCPU(m_CommandQue); const unsigned int dimension = m_Input->GetDimension(); unsigned int* dimensions = m_Input->GetDimensions(); const mitk::SlicedGeometry3D::Pointer p_slg = m_Input->GetMITKImage()->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); } MITK_DEBUG << "Image Initialized."; return m_Output->GetMITKImage(); } mitk::PixelType mitk::OclImageToImageFilter::GetOutputType() { // get the current image format from the input image const cl_image_format* currentImFormat = this->m_Input->GetPixelType(); // return the value according to the current channel type switch( currentImFormat->image_channel_data_type ) { case CL_UNORM_INT8: return mitk::MakeScalarPixelType(); case CL_UNSIGNED_INT8: return mitk::MakeScalarPixelType(); case CL_UNORM_INT16: return mitk::MakeScalarPixelType(); default: return mitk::MakeScalarPixelType(); } } int mitk::OclImageToImageFilter::GetBytesPerElem() { return (this->m_CurrentType + 1); } bool mitk::OclImageToImageFilter::InitExec(cl_kernel ckKernel) { cl_int clErr = 0; if( m_Input.IsNull() ) mitkThrow() << "Input image is null."; // get image size once const unsigned int uiImageWidth = m_Input->GetDimension(0); const unsigned int uiImageHeight = m_Input->GetDimension(1); const unsigned int uiImageDepth = m_Input->GetDimension(2); // compute work sizes this->SetWorkingSize( 8, uiImageWidth, 8, uiImageHeight , 8, uiImageDepth ); cl_mem clBuffIn = m_Input->GetGPUImage(this->m_CommandQue); cl_mem clBuffOut = m_Output->GetGPUBuffer(); if (!clBuffIn) { if ( m_Input->TransferDataToGPU(m_CommandQue) != CL_SUCCESS ) { mitkThrow()<< "Could not create / initialize gpu image."; } clBuffIn = m_Input->GetGPUImage(m_CommandQue); } // output image not initialized if (!clBuffOut) { //TODO bpp, or SetImageWidth/Height/... MITK_DEBUG << "Create GPU Image call " << uiImageWidth<< "x"<CreateGPUImage(uiImageWidth, uiImageHeight, uiImageDepth, this->m_CurrentType + 1); } 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 ); } bool mitk::OclImageToImageFilter::InitExec(cl_kernel ckKernel, unsigned int* dimensions) { cl_int clErr = 0; if( m_Input.IsNull() ) mitkThrow() << "Input image is null."; // get image size once const unsigned int uiImageWidth = dimensions[0]; const unsigned int uiImageHeight = dimensions[1]; - const unsigned int uiImageDepth = dimensions[2]; + const unsigned int uiImageDepth = dimensions[2]+1; // compute work sizes this->SetWorkingSize( 8, uiImageWidth, 8, uiImageHeight , 8, uiImageDepth ); cl_mem clBuffIn = m_Input->GetGPUImage(this->m_CommandQue); cl_mem clBuffOut = m_Output->GetGPUBuffer(); if (!clBuffIn) { if ( m_Input->TransferDataToGPU(m_CommandQue) != CL_SUCCESS ) { mitkThrow()<< "Could not create / initialize gpu image."; } clBuffIn = m_Input->GetGPUImage(m_CommandQue); } // output image not initialized //TODO bpp, or SetImageWidth/Height/... MITK_INFO << "Create GPU Image call " << uiImageWidth<< "x"<CreateGPUImage(uiImageWidth, uiImageHeight, uiImageDepth, this->m_CurrentType + 1); 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/OpenCL/mitkOclUtils.cpp b/Modules/OpenCL/mitkOclUtils.cpp index d743b3d494..a8fa4c0c53 100644 --- a/Modules/OpenCL/mitkOclUtils.cpp +++ b/Modules/OpenCL/mitkOclUtils.cpp @@ -1,572 +1,576 @@ /*=================================================================== 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 "mitkOclUtils.h" #include "mitkLogMacros.h" #include #include unsigned int iDivUp(unsigned int dividend, unsigned int divisor){ return (dividend % divisor == 0) ? (dividend / divisor) : (dividend / divisor + 1); } cl_int oclGetPlatformID(cl_platform_id* selectedPlatform) { cl_uint num_platforms = 0; cl_platform_id* clPlatformIDs; cl_int ciErrNum = 0; ciErrNum = clGetPlatformIDs( 0, nullptr, &num_platforms); if ( ciErrNum != CL_SUCCESS) { MITK_ERROR<<" Error " << ciErrNum << " in clGetPlatformIDs() \n"; throw std::bad_exception(); } else { clPlatformIDs = new cl_platform_id[num_platforms]; ciErrNum = clGetPlatformIDs( num_platforms, clPlatformIDs, nullptr); if(ciErrNum == CL_SUCCESS) { *selectedPlatform = clPlatformIDs[0]; } } return CL_SUCCESS; } void oclPrintMemObjectInfo(cl_mem memobj) { cl_int clErr = 0; MITK_INFO << "Examining cl_mem object: " << memobj << "\n------------------\n"; // CL_MEM_TYPE cl_mem_object_type objtype; clErr = clGetMemObjectInfo( memobj, CL_MEM_TYPE, sizeof(cl_mem_object_type),&objtype, nullptr); CHECK_OCL_ERR( clErr ); switch(objtype) { case CL_MEM_OBJECT_BUFFER: MITK_INFO << "CL_MEM_TYPE \t" << "BUFFER_OBJ" << "\n"; break; case CL_MEM_OBJECT_IMAGE2D: MITK_INFO << "CL_MEM_TYPE \t" << "2D IMAGE" << "\n"; break; case CL_MEM_OBJECT_IMAGE3D: MITK_INFO << "CL_MEM_TYPE \t" << "3D IMAGE" << "\n"; break; default: MITK_INFO << "CL_MEM_TYPE \t" << "[could not resolve]" << "\n"; break; } // CL_MEM_FLAGS cl_mem_flags flags; clErr = clGetMemObjectInfo( memobj, CL_MEM_FLAGS, sizeof(cl_mem_flags),&flags, nullptr); CHECK_OCL_ERR( clErr ); switch(flags) { case CL_MEM_READ_ONLY: MITK_INFO << "CL_MEM_FLAGS \t" << "CL_MEM_READ_ONLY" << "\n"; break; case CL_MEM_WRITE_ONLY: MITK_INFO << "CL_MEM_FLAGS \t" << "CL_MEM_WRITE_ONLY" << "\n"; break; case CL_MEM_READ_WRITE: MITK_INFO << "CL_MEM_FLAGS \t" << "CL_MEM_READ_WRITE" << "\n"; break; default: MITK_INFO << "CL_MEM_FLAGS \t" << "not resolved, " << flags << "\n"; break; } // get CL_MEM_SIZE size_t memsize; clErr = clGetMemObjectInfo( memobj, CL_MEM_SIZE, sizeof(memsize),&memsize, nullptr); CHECK_OCL_ERR( clErr ); MITK_INFO << "CL_MEM_SIZE \t" << memsize << "\n"; // get CL_MEM_HOST_PTR float *hostptr; clErr = clGetMemObjectInfo( memobj, CL_MEM_HOST_PTR, sizeof(void*), (void*) &hostptr, nullptr); CHECK_OCL_ERR( clErr ); MITK_INFO << "CL_MEM_HOST_PTR \t" << hostptr << "\n"; // get CL_CONTEXT cl_context gpuctxt; clErr = clGetMemObjectInfo( memobj, CL_MEM_CONTEXT, sizeof(cl_context), &gpuctxt, nullptr); CHECK_OCL_ERR( clErr ); MITK_INFO << "CL_CONTEXT \t\t" << gpuctxt << "\n"; // get CL_MEM_REFERENCE_COUNT cl_uint refs; clErr = clGetMemObjectInfo( memobj, CL_MEM_REFERENCE_COUNT, sizeof(cl_uint), &refs, nullptr); CHECK_OCL_ERR(clErr); MITK_INFO << "CL_REF_COUNT \t" << refs << "\n"; MITK_INFO << "================== \n" << std::endl; } void oclPrintDeviceInfo(cl_device_id device) { char device_string[1024]; clGetDeviceInfo(device, CL_DEVICE_NAME, sizeof(device_string), &device_string, nullptr); MITK_INFO("ocl.log")<< " Device : " << device_string; // CL_DEVICE_INFO cl_device_type type; clGetDeviceInfo(device, CL_DEVICE_TYPE, sizeof(type), &type, nullptr); if( type & CL_DEVICE_TYPE_CPU ) MITK_INFO("ocl.log")<<" CL_DEVICE_TYPE: CL_DEVICE_TYPE_CPU"; if( type & CL_DEVICE_TYPE_GPU ) MITK_INFO("ocl.log")<<" CL_DEVICE_TYPE: CL_DEVICE_TYPE_GPU"; if( type & CL_DEVICE_TYPE_ACCELERATOR ) MITK_INFO("ocl.log")<<" CL_DEVICE_TYPE: CL_DEVICE_TYPE_ACCELERATOR"; if( type & CL_DEVICE_TYPE_DEFAULT ) MITK_INFO("ocl.log")<<" CL_DEVICE_TYPE: CL_DEVICE_TYPE_DEFAULT"; // CL_DEVICE_MAX_COMPUTE_UNITS cl_uint compute_units; clGetDeviceInfo(device, CL_DEVICE_MAX_COMPUTE_UNITS, sizeof(compute_units), &compute_units, nullptr); MITK_INFO("ocl.log")<<" CL_DEVICE_MAX_COMPUTE_UNITS:" << compute_units; // CL_DEVICE_MAX_WORK_GROUP_SIZE size_t workitem_size[3]; clGetDeviceInfo(device, CL_DEVICE_MAX_WORK_ITEM_SIZES, sizeof(workitem_size), &workitem_size, nullptr); MITK_INFO("ocl.log")<<" CL_DEVICE_MAX_WORK_ITEM_SIZES:\t"<< workitem_size[0]<< workitem_size[1]<< workitem_size[2]; // CL_DEVICE_MAX_WORK_GROUP_SIZE size_t workgroup_size; clGetDeviceInfo(device, CL_DEVICE_MAX_WORK_GROUP_SIZE, sizeof(workgroup_size), &workgroup_size, nullptr); MITK_INFO("ocl.log")<<" CL_DEVICE_MAX_WORK_GROUP_SIZE:" << workgroup_size; // CL_DEVICE_MAX_CLOCK_FREQUENCY cl_uint clock_frequency; clGetDeviceInfo(device, CL_DEVICE_MAX_CLOCK_FREQUENCY, sizeof(clock_frequency), &clock_frequency, nullptr); MITK_INFO("ocl.log")<<" CL_DEVICE_MAX_CLOCK_FREQUENCY:"<< clock_frequency / 1000; // CL_DEVICE_IMAGE_SUPPORT cl_bool image_support; clGetDeviceInfo(device, CL_DEVICE_IMAGE_SUPPORT, sizeof(image_support), &image_support, nullptr); MITK_INFO("ocl.log")<<" CL_DEVICE_IMAGE_SUPPORT:\t" << image_support; // CL_DEVICE_GLOBAL_MEM_SIZE cl_ulong mem_size; clGetDeviceInfo(device, CL_DEVICE_GLOBAL_MEM_SIZE, sizeof(mem_size), &mem_size, nullptr); MITK_INFO("ocl.log")<<" CL_DEVICE_GLOBAL_MEM_SIZE:\t\t"<<(unsigned int)(mem_size / (1024 * 1024))<<"Mbytes"; // CL_DEVICE_LOCAL_MEM_SIZE clGetDeviceInfo(device, CL_DEVICE_LOCAL_MEM_SIZE, sizeof(mem_size), &mem_size, nullptr); MITK_INFO("ocl.log")<<" CL_DEVICE_LOCAL_MEM_SIZE:\t\t"<< (unsigned int)(mem_size / (1024)) <<"KByte\n"; + // CL_DEVICE_MAX_CONSTANT_BUFFER_SIZE + clGetDeviceInfo(device, CL_DEVICE_MAX_CONSTANT_BUFFER_SIZE, sizeof(mem_size), &mem_size, nullptr); + MITK_INFO("ocl.log") << " CL_DEVICE_MAX_CONSTANT_BUFFER_SIZE:\t\t" << (unsigned int)(mem_size / (1024)) << "KByte"; + //check for image support properties clGetDeviceInfo(device, CL_DEVICE_IMAGE2D_MAX_WIDTH, sizeof(workgroup_size), &workgroup_size, nullptr); MITK_INFO("ocl.log")<<" CL_DEVICE_IMAGE2D_MAX_WIDTH:\t" << workgroup_size; clGetDeviceInfo(device, CL_DEVICE_IMAGE2D_MAX_HEIGHT, sizeof(workgroup_size), &workgroup_size, nullptr); MITK_INFO("ocl.log")<<" CL_DEVICE_IMAGE2D_MAX_HEIGHT:\t" << workgroup_size; clGetDeviceInfo(device, CL_DEVICE_IMAGE3D_MAX_WIDTH, sizeof(workgroup_size), &workgroup_size, nullptr); MITK_INFO("ocl.log")<<" CL_DEVICE_IMAGE3D_MAX_WIDTH:\t" << workgroup_size; clGetDeviceInfo(device, CL_DEVICE_IMAGE3D_MAX_HEIGHT, sizeof(workgroup_size), &workgroup_size, nullptr); MITK_INFO("ocl.log")<<" CL_DEVICE_IMAGE3D_MAX_HEIGHT:\t" << workgroup_size; clGetDeviceInfo(device, CL_DEVICE_IMAGE3D_MAX_DEPTH, sizeof(workgroup_size), &workgroup_size, nullptr); MITK_INFO("ocl.log")<<" CL_DEVICE_IMAGE3D_MAX_DEPTH:\t" << workgroup_size; // CL_DEVICE_QUEUE_PROPERTIES cl_command_queue_properties queue_properties; clGetDeviceInfo(device, CL_DEVICE_QUEUE_PROPERTIES, sizeof(queue_properties), &queue_properties, nullptr); if( queue_properties & CL_QUEUE_OUT_OF_ORDER_EXEC_MODE_ENABLE ) MITK_INFO("ocl.log")<<" CL_DEVICE_QUEUE_PROPERTIES:\t\t"<< "CL_QUEUE_OUT_OF_ORDER_EXEC_MODE_ENABLE"; if( queue_properties & CL_QUEUE_PROFILING_ENABLE ) MITK_INFO("ocl.log")<<" CL_DEVICE_QUEUE_PROPERTIES:\t\t"<< "CL_QUEUE_PROFILING_ENABLE"; } std::string GetOclErrorAsString( int _clErr ) { std::string returnString("unkown error number: "+std::to_string(_clErr)+" \n"); switch(_clErr) { case CL_SUCCESS: returnString = "CL_SUCCESS\n"; break; case CL_DEVICE_NOT_FOUND: returnString = "CL_DEVICE_NOT_FOUND\n"; break; case CL_DEVICE_NOT_AVAILABLE: returnString = "CL_DEVICE_NOT_AVAILABLE\n"; break; /*case CL_DEVICE_COMPILER_NOT_AVAILABLE: returnString = "CL_DEVICE_COMPILER_NOT_AVAILABLE\n"; break; */ case CL_MEM_OBJECT_ALLOCATION_FAILURE : returnString = "CL_MEM_OBJECT_ALLOCATION_FAILURE\n"; break; case CL_OUT_OF_RESOURCES: returnString = "CL_OUT_OF_RESOURCES\n"; break; case CL_OUT_OF_HOST_MEMORY: returnString = "CL_OUT_OF_HOST_MEMORY\n"; break; case CL_PROFILING_INFO_NOT_AVAILABLE: returnString = "CL_PROFILING_INFO_NOT_AVAILABLE\n"; break; case CL_MEM_COPY_OVERLAP: returnString = "CL_MEM_COPY_OVERLAP\n"; break; case CL_IMAGE_FORMAT_MISMATCH: returnString = "CL_IMAGE_FORMAT_MISMATCH\n"; break; case CL_IMAGE_FORMAT_NOT_SUPPORTED: returnString = "CL_IMAGE_FORMAT_NOT_SUPPORTED\n"; break; case CL_BUILD_PROGRAM_FAILURE: returnString = "CL_BUILD_PROGRAM_FAILURE\n"; break; case CL_MAP_FAILURE: returnString = "CL_MAP_FAILURE\n"; break; case CL_INVALID_VALUE: returnString = "CL_INVALID_VALUE\n"; break; case CL_INVALID_DEVICE_TYPE: returnString = "CL_INVALID_DEVICE_TYPE\n"; break; case CL_INVALID_PLATFORM: returnString = "CL_INVALID_PLATFORM\n"; break; case CL_INVALID_DEVICE: returnString = "CL_INVALID_DEVICE\n"; break; case CL_INVALID_CONTEXT : returnString = "CL_INVALID_CONTEXT\n"; break; case CL_INVALID_QUEUE_PROPERTIES: returnString = "CL_INVALID_QUEUE_PROPERTIES\n"; break; case CL_INVALID_COMMAND_QUEUE: returnString = "CL_INVALID_COMMAND_QUEUE\n"; break; case CL_INVALID_HOST_PTR: returnString = "CL_INVALID_HOST_PTR\n"; break; case CL_INVALID_MEM_OBJECT: returnString = "CL_INVALID_MEM_OBJECT\n"; break; case CL_INVALID_IMAGE_FORMAT_DESCRIPTOR: returnString = "CL_INVALID_IMAGE_FORMAT_DESCRIPTOR\n"; break; case CL_INVALID_IMAGE_SIZE: returnString = "CL_INVALID_IMAGE_SIZE\n"; break; case CL_INVALID_SAMPLER : returnString = "CL_INVALID_SAMPLER\n"; break; case CL_INVALID_BINARY: returnString = "CL_INVALID_BINARY\n"; break; case CL_INVALID_BUILD_OPTIONS: returnString = "CL_INVALID_BUILD_OPTIONS\n"; break; case CL_INVALID_PROGRAM: returnString = "CL_INVALID_PROGRAM\n"; break; case CL_INVALID_PROGRAM_EXECUTABLE: returnString = "CL_INVALID_PROGRAM_EXECUTABLE\n"; break; case CL_INVALID_KERNEL_NAME: returnString = "CL_INVALID_KERNEL_NAME\n"; break; case CL_INVALID_KERNEL_DEFINITION: returnString = "CL_INVALID_KERNEL_DEFINITION\n"; break; case CL_INVALID_KERNEL : returnString = "CL_INVALID_KERNEL\n"; break; case CL_INVALID_ARG_INDEX : returnString = "CL_INVALID_ARG_INDEX\n"; break; case CL_INVALID_ARG_VALUE : returnString = "CL_INVALID_ARG_VALUE\n"; break; case CL_INVALID_ARG_SIZE : returnString = "CL_INVALID_ARG_SIZE\n"; break; case CL_INVALID_KERNEL_ARGS : returnString = "CL_INVALID_KERNEL_ARGS\n"; break; case CL_INVALID_WORK_DIMENSION: returnString = "CL_INVALID_WORK_DIMENSION\n"; break; case CL_INVALID_WORK_GROUP_SIZE: returnString = "CL_INVALID_WORK_GROUP_SIZE\n"; break; case CL_INVALID_WORK_ITEM_SIZE: returnString = "CL_INVALID_WORK_ITEM_SIZE\n"; break; case CL_INVALID_GLOBAL_OFFSET: returnString = "CL_INVALID_GLOBAL_OFFSET\n"; break; case CL_INVALID_EVENT_WAIT_LIST: returnString = "CL_INVALID_EVENT_WAIT_LIST\n"; break; case CL_INVALID_EVENT: returnString = "CL_INVALID_EVENT\n"; break; case CL_INVALID_OPERATION: returnString = "CL_INVALID_OPERATION\n"; break; case CL_INVALID_GL_OBJECT: returnString = "CL_INVALID_GL_OBJECT\n"; break; case CL_INVALID_BUFFER_SIZE : returnString = "CL_INVALID_BUFFER_SIZE\n"; break; case CL_INVALID_MIP_LEVEL : returnString = "CL_INVALID_MIP_LEVEL\n"; break; default: break; } return returnString; } void GetOclError(int _clErr) { if(_clErr == CL_SUCCESS) MITK_WARN << "Called GetOclErr() with no error value: [CL_SUCCESS]"; else MITK_ERROR << GetOclErrorAsString(_clErr); } bool oclCheckError(int _err, const char* filepath, int lineno) { if (_err) { MITK_ERROR<< "OpenCL Error at " << filepath <<":"<< lineno; GetOclError(_err); return 0; } return 1; } void GetSupportedImageFormats(cl_context _context, cl_mem_object_type _type) { const unsigned int entries = 500; cl_image_format* formats = new cl_image_format[entries]; cl_uint _written = 0; // OpenCL constant to catch error IDs cl_int ciErr1; // Get list of supported image formats for READ_ONLY access ciErr1 = clGetSupportedImageFormats( _context, CL_MEM_READ_ONLY, _type, entries, formats, &_written); CHECK_OCL_ERR(ciErr1); MITK_INFO << "Supported Image Formats, Image: CL_MEM_READ_ONLY \n"; for (unsigned int i=0; i<_written; i++) { MITK_INFO<< "ChannelType: " << GetImageTypeAsString(formats[i].image_channel_data_type) << "| ChannelOrder: "<< GetImageTypeAsString(formats[i].image_channel_order) <<"\n"; } _written = 0; // Get list of supported image formats for READ_WRITE access ciErr1 = clGetSupportedImageFormats( _context, CL_MEM_READ_WRITE, _type, entries, formats, &_written); CHECK_OCL_ERR(ciErr1); MITK_INFO << "Supported Image Formats, Image: CL_MEM_READ_WRITE (found: " << _written <<") \n"; for (unsigned int i=0; i<_written; i++) { MITK_INFO<< "ChannelType: " << GetImageTypeAsString(formats[i].image_channel_data_type) << "| ChannelOrder: "<< GetImageTypeAsString(formats[i].image_channel_order) <<"\n"; } _written = 0; // Get list of supported image formats for WRITE_ONLY access ciErr1 = clGetSupportedImageFormats( _context, CL_MEM_WRITE_ONLY, _type, entries, formats, &_written); CHECK_OCL_ERR(ciErr1); MITK_INFO << "Supported Image Formats, Image: CL_MEM_WRITE_ONLY (found: " << _written <<") \n"; for (unsigned int i=0; i<_written; i++) { MITK_INFO<< "ChannelType: " << GetImageTypeAsString(formats[i].image_channel_data_type) << "| ChannelOrder: "<< GetImageTypeAsString(formats[i].image_channel_order) <<"\n"; } } std::string GetImageTypeAsString( const unsigned int _in) { switch(_in) { case CL_R: return "CL_R "; break; case CL_A: return "CL_A "; break; case CL_RG: return "CL_RG "; break; case CL_RA: return "CL_RA "; break; case CL_RGB: return "CL_RGB "; break; case CL_RGBA: return "CL_RGBA "; break; case CL_BGRA: return "CL_BGRA "; break; case CL_ARGB: return "CL_ARGB "; break; case CL_INTENSITY: return "CL_INTENSITY "; break; case CL_LUMINANCE: return "CL_LUMINANCE "; break; case CL_SNORM_INT8: return "CL_SNORM_INT8 "; break; case CL_SNORM_INT16: return "CL_SNORM_INT16 "; break; case CL_UNORM_INT8: return "CL_UNORM_INT8 "; break; case CL_UNORM_INT16: return "CL_UNORM_INT16 "; break; case CL_UNORM_SHORT_565: return "CL_UNORM_SHORT_565 "; break; case CL_UNORM_SHORT_555: return "CL_UNORM_SHORT_555 "; break; case CL_UNORM_INT_101010: return "CL_UNORM_INT_101010 "; break; case CL_SIGNED_INT8: return "CL_SIGNED_INT8 "; break; case CL_SIGNED_INT16: return "CL_SIGNED_INT16 "; break; case CL_SIGNED_INT32: return "CL_SIGNED_INT32 "; break; case CL_UNSIGNED_INT8: return "CL_UNSIGNED_INT8 "; break; case CL_UNSIGNED_INT16: return "CL_UNSIGNED_INT16 "; break; case CL_UNSIGNED_INT32: return "CL_UNSIGNED_INT32 "; break; case CL_HALF_FLOAT: return "CL_HALF_FLOAT "; break; case CL_FLOAT: return "CL_FLOAT "; break; default: return "--"; break; } } void oclLogBinary(cl_program clProg, cl_device_id clDev) { // Grab the number of devices associated with the program cl_uint num_devices; clGetProgramInfo(clProg, CL_PROGRAM_NUM_DEVICES, sizeof(cl_uint), &num_devices, nullptr); // Grab the device ids cl_device_id* devices = (cl_device_id*) malloc(num_devices * sizeof(cl_device_id)); clGetProgramInfo(clProg, CL_PROGRAM_DEVICES, num_devices * sizeof(cl_device_id), devices, 0); // Grab the sizes of the binaries size_t* binary_sizes = (size_t*)malloc(num_devices * sizeof(size_t)); clGetProgramInfo(clProg, CL_PROGRAM_BINARY_SIZES, num_devices * sizeof(size_t), binary_sizes, nullptr); // Now get the binaries char** ptx_code = (char**)malloc(num_devices * sizeof(char*)); for( unsigned int i=0; iAddSourceFile("BModeAbs.cl"); - this->AddSourceFile("BModeAbsLog.cl"); - - this->m_FilterID = "PixelCalculation"; -} - -mitk::PhotoacousticBModeFilter::~PhotoacousticBModeFilter() -{ - if (this->m_PixelCalculation) - { - clReleaseKernel(m_PixelCalculation); - } -} - -void mitk::PhotoacousticBModeFilter::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::PhotoacousticBModeFilter::Execute() -{ - try - { - this->InitExec(this->m_PixelCalculation); - } - catch (const mitk::Exception& e) - { - MITK_ERROR << "Catched exception while initializing filter: " << e.what(); - return; - } - - // 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::PhotoacousticBModeFilter::GetModule() -{ - return us::GetModuleContext()->GetModule(); -} - -bool mitk::PhotoacousticBModeFilter::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::PhotoacousticBModeFilter::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); -} diff --git a/Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticOCLBeamformer.cpp b/Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticOCLBeamformer.cpp deleted file mode 100644 index 114810ddc9..0000000000 --- a/Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticOCLBeamformer.cpp +++ /dev/null @@ -1,177 +0,0 @@ -/*=================================================================== - -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 ) -{ - this->AddSourceFile("DASQuadratic.cl"); - this->AddSourceFile("DMASQuadratic.cl"); - this->AddSourceFile("DASspherical.cl"); - this->AddSourceFile("DMASspherical.cl"); - - 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); - } - 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)); - - 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::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); -} - -mitk::Image::Pointer mitk::PhotoacousticOCLBeamformer::GetOutput() -{ - if (m_Output->IsModified(GPU_DATA)) - { - void* pData = m_Output->TransferDataToCPU(m_CommandQue); - - const unsigned int dimension = 3; - unsigned int* dimensions = m_OutputDim; - - const mitk::SlicedGeometry3D::Pointer p_slg = m_Input->GetMITKImage()->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); - } - - MITK_DEBUG << "Image Initialized."; - - return m_Output->GetMITKImage(); -} diff --git a/Modules/PhotoacousticsAlgorithms/CMakeLists.txt b/Modules/PhotoacousticsAlgorithms/CMakeLists.txt index a1bd2a69cf..7b8587334b 100644 --- a/Modules/PhotoacousticsAlgorithms/CMakeLists.txt +++ b/Modules/PhotoacousticsAlgorithms/CMakeLists.txt @@ -1,8 +1,16 @@ +set(dependencies_list MitkCore MitkAlgorithmsExt) + +IF(MITK_USE_OpenCL) + add_definitions(-DPHOTOACOUSTICS_USE_GPU) + set(dependencies_list ${dependencies_list} MitkOpenCL) + message("Using OpenCL in PhotoacousticAlgorithms") +ENDIF(MITK_USE_OpenCL) + MITK_CREATE_MODULE( SUBPROJECTS - DEPENDS MitkCore MitkAlgorithmsExt MitkOpenCL + DEPENDS ${dependencies_list} #AUTOLOAD_WITH MitkCore INCLUDE_DIRS PUBLIC Algorithms/ITKUltrasound Algorithms Algorithms/OCL INTERNAL_INCLUDE_DIRS ${INCLUDE_DIRS_INTERNAL} PACKAGE_DEPENDS ITK|ITKFFT+ITKImageCompose+ITKImageIntensity -) +) \ No newline at end of file diff --git a/Modules/PhotoacousticsAlgorithms/Algorithms/ITKUltrasound/LICENSE b/Modules/PhotoacousticsAlgorithms/ITKFilter/ITKUltrasound/LICENSE similarity index 100% rename from Modules/PhotoacousticsAlgorithms/Algorithms/ITKUltrasound/LICENSE rename to Modules/PhotoacousticsAlgorithms/ITKFilter/ITKUltrasound/LICENSE diff --git a/Modules/PhotoacousticsAlgorithms/Algorithms/ITKUltrasound/README.rst b/Modules/PhotoacousticsAlgorithms/ITKFilter/ITKUltrasound/README.rst similarity index 100% rename from Modules/PhotoacousticsAlgorithms/Algorithms/ITKUltrasound/README.rst rename to Modules/PhotoacousticsAlgorithms/ITKFilter/ITKUltrasound/README.rst diff --git a/Modules/PhotoacousticsAlgorithms/Algorithms/ITKUltrasound/itkAnalyticSignalImageFilter.h b/Modules/PhotoacousticsAlgorithms/ITKFilter/ITKUltrasound/itkAnalyticSignalImageFilter.h similarity index 100% rename from Modules/PhotoacousticsAlgorithms/Algorithms/ITKUltrasound/itkAnalyticSignalImageFilter.h rename to Modules/PhotoacousticsAlgorithms/ITKFilter/ITKUltrasound/itkAnalyticSignalImageFilter.h diff --git a/Modules/PhotoacousticsAlgorithms/Algorithms/ITKUltrasound/itkAnalyticSignalImageFilter.hxx b/Modules/PhotoacousticsAlgorithms/ITKFilter/ITKUltrasound/itkAnalyticSignalImageFilter.hxx similarity index 100% rename from Modules/PhotoacousticsAlgorithms/Algorithms/ITKUltrasound/itkAnalyticSignalImageFilter.hxx rename to Modules/PhotoacousticsAlgorithms/ITKFilter/ITKUltrasound/itkAnalyticSignalImageFilter.hxx diff --git a/Modules/PhotoacousticsAlgorithms/Algorithms/ITKUltrasound/itkBModeImageFilter.h b/Modules/PhotoacousticsAlgorithms/ITKFilter/ITKUltrasound/itkBModeImageFilter.h similarity index 100% rename from Modules/PhotoacousticsAlgorithms/Algorithms/ITKUltrasound/itkBModeImageFilter.h rename to Modules/PhotoacousticsAlgorithms/ITKFilter/ITKUltrasound/itkBModeImageFilter.h diff --git a/Modules/PhotoacousticsAlgorithms/Algorithms/ITKUltrasound/itkBModeImageFilter.hxx b/Modules/PhotoacousticsAlgorithms/ITKFilter/ITKUltrasound/itkBModeImageFilter.hxx similarity index 100% rename from Modules/PhotoacousticsAlgorithms/Algorithms/ITKUltrasound/itkBModeImageFilter.hxx rename to Modules/PhotoacousticsAlgorithms/ITKFilter/ITKUltrasound/itkBModeImageFilter.hxx diff --git a/Modules/PhotoacousticsAlgorithms/Algorithms/ITKUltrasound/itkCurvilinearArraySpecialCoordinatesImage.h b/Modules/PhotoacousticsAlgorithms/ITKFilter/ITKUltrasound/itkCurvilinearArraySpecialCoordinatesImage.h similarity index 100% rename from Modules/PhotoacousticsAlgorithms/Algorithms/ITKUltrasound/itkCurvilinearArraySpecialCoordinatesImage.h rename to Modules/PhotoacousticsAlgorithms/ITKFilter/ITKUltrasound/itkCurvilinearArraySpecialCoordinatesImage.h diff --git a/Modules/PhotoacousticsAlgorithms/Algorithms/ITKUltrasound/itkCurvilinearArraySpecialCoordinatesImage.hxx b/Modules/PhotoacousticsAlgorithms/ITKFilter/ITKUltrasound/itkCurvilinearArraySpecialCoordinatesImage.hxx similarity index 100% rename from Modules/PhotoacousticsAlgorithms/Algorithms/ITKUltrasound/itkCurvilinearArraySpecialCoordinatesImage.hxx rename to Modules/PhotoacousticsAlgorithms/ITKFilter/ITKUltrasound/itkCurvilinearArraySpecialCoordinatesImage.hxx diff --git a/Modules/PhotoacousticsAlgorithms/Algorithms/ITKUltrasound/itkFFT1DComplexConjugateToRealImageFilter.h b/Modules/PhotoacousticsAlgorithms/ITKFilter/ITKUltrasound/itkFFT1DComplexConjugateToRealImageFilter.h similarity index 100% rename from Modules/PhotoacousticsAlgorithms/Algorithms/ITKUltrasound/itkFFT1DComplexConjugateToRealImageFilter.h rename to Modules/PhotoacousticsAlgorithms/ITKFilter/ITKUltrasound/itkFFT1DComplexConjugateToRealImageFilter.h diff --git a/Modules/PhotoacousticsAlgorithms/Algorithms/ITKUltrasound/itkFFT1DComplexConjugateToRealImageFilter.hxx b/Modules/PhotoacousticsAlgorithms/ITKFilter/ITKUltrasound/itkFFT1DComplexConjugateToRealImageFilter.hxx similarity index 100% rename from Modules/PhotoacousticsAlgorithms/Algorithms/ITKUltrasound/itkFFT1DComplexConjugateToRealImageFilter.hxx rename to Modules/PhotoacousticsAlgorithms/ITKFilter/ITKUltrasound/itkFFT1DComplexConjugateToRealImageFilter.hxx diff --git a/Modules/PhotoacousticsAlgorithms/Algorithms/ITKUltrasound/itkFFT1DComplexToComplexImageFilter.h b/Modules/PhotoacousticsAlgorithms/ITKFilter/ITKUltrasound/itkFFT1DComplexToComplexImageFilter.h similarity index 100% rename from Modules/PhotoacousticsAlgorithms/Algorithms/ITKUltrasound/itkFFT1DComplexToComplexImageFilter.h rename to Modules/PhotoacousticsAlgorithms/ITKFilter/ITKUltrasound/itkFFT1DComplexToComplexImageFilter.h diff --git a/Modules/PhotoacousticsAlgorithms/Algorithms/ITKUltrasound/itkFFT1DComplexToComplexImageFilter.hxx b/Modules/PhotoacousticsAlgorithms/ITKFilter/ITKUltrasound/itkFFT1DComplexToComplexImageFilter.hxx similarity index 100% rename from Modules/PhotoacousticsAlgorithms/Algorithms/ITKUltrasound/itkFFT1DComplexToComplexImageFilter.hxx rename to Modules/PhotoacousticsAlgorithms/ITKFilter/ITKUltrasound/itkFFT1DComplexToComplexImageFilter.hxx diff --git a/Modules/PhotoacousticsAlgorithms/Algorithms/ITKUltrasound/itkFFT1DRealToComplexConjugateImageFilter.h b/Modules/PhotoacousticsAlgorithms/ITKFilter/ITKUltrasound/itkFFT1DRealToComplexConjugateImageFilter.h similarity index 100% rename from Modules/PhotoacousticsAlgorithms/Algorithms/ITKUltrasound/itkFFT1DRealToComplexConjugateImageFilter.h rename to Modules/PhotoacousticsAlgorithms/ITKFilter/ITKUltrasound/itkFFT1DRealToComplexConjugateImageFilter.h diff --git a/Modules/PhotoacousticsAlgorithms/Algorithms/ITKUltrasound/itkFFT1DRealToComplexConjugateImageFilter.hxx b/Modules/PhotoacousticsAlgorithms/ITKFilter/ITKUltrasound/itkFFT1DRealToComplexConjugateImageFilter.hxx similarity index 100% rename from Modules/PhotoacousticsAlgorithms/Algorithms/ITKUltrasound/itkFFT1DRealToComplexConjugateImageFilter.hxx rename to Modules/PhotoacousticsAlgorithms/ITKFilter/ITKUltrasound/itkFFT1DRealToComplexConjugateImageFilter.hxx diff --git a/Modules/PhotoacousticsAlgorithms/Algorithms/ITKUltrasound/itkRegionFromReferenceImageFilter.h b/Modules/PhotoacousticsAlgorithms/ITKFilter/ITKUltrasound/itkRegionFromReferenceImageFilter.h similarity index 100% rename from Modules/PhotoacousticsAlgorithms/Algorithms/ITKUltrasound/itkRegionFromReferenceImageFilter.h rename to Modules/PhotoacousticsAlgorithms/ITKFilter/ITKUltrasound/itkRegionFromReferenceImageFilter.h diff --git a/Modules/PhotoacousticsAlgorithms/Algorithms/ITKUltrasound/itkRegionFromReferenceImageFilter.hxx b/Modules/PhotoacousticsAlgorithms/ITKFilter/ITKUltrasound/itkRegionFromReferenceImageFilter.hxx similarity index 100% rename from Modules/PhotoacousticsAlgorithms/Algorithms/ITKUltrasound/itkRegionFromReferenceImageFilter.hxx rename to Modules/PhotoacousticsAlgorithms/ITKFilter/ITKUltrasound/itkRegionFromReferenceImageFilter.hxx diff --git a/Modules/PhotoacousticsAlgorithms/Algorithms/ITKUltrasound/itkSpectra1DImageFilter.h b/Modules/PhotoacousticsAlgorithms/ITKFilter/ITKUltrasound/itkSpectra1DImageFilter.h similarity index 100% rename from Modules/PhotoacousticsAlgorithms/Algorithms/ITKUltrasound/itkSpectra1DImageFilter.h rename to Modules/PhotoacousticsAlgorithms/ITKFilter/ITKUltrasound/itkSpectra1DImageFilter.h diff --git a/Modules/PhotoacousticsAlgorithms/Algorithms/ITKUltrasound/itkSpectra1DImageFilter.hxx b/Modules/PhotoacousticsAlgorithms/ITKFilter/ITKUltrasound/itkSpectra1DImageFilter.hxx similarity index 100% rename from Modules/PhotoacousticsAlgorithms/Algorithms/ITKUltrasound/itkSpectra1DImageFilter.hxx rename to Modules/PhotoacousticsAlgorithms/ITKFilter/ITKUltrasound/itkSpectra1DImageFilter.hxx diff --git a/Modules/PhotoacousticsAlgorithms/Algorithms/ITKUltrasound/itkSpectra1DSupportWindowImageFilter.h b/Modules/PhotoacousticsAlgorithms/ITKFilter/ITKUltrasound/itkSpectra1DSupportWindowImageFilter.h similarity index 100% rename from Modules/PhotoacousticsAlgorithms/Algorithms/ITKUltrasound/itkSpectra1DSupportWindowImageFilter.h rename to Modules/PhotoacousticsAlgorithms/ITKFilter/ITKUltrasound/itkSpectra1DSupportWindowImageFilter.h diff --git a/Modules/PhotoacousticsAlgorithms/Algorithms/ITKUltrasound/itkSpectra1DSupportWindowImageFilter.hxx b/Modules/PhotoacousticsAlgorithms/ITKFilter/ITKUltrasound/itkSpectra1DSupportWindowImageFilter.hxx similarity index 100% rename from Modules/PhotoacousticsAlgorithms/Algorithms/ITKUltrasound/itkSpectra1DSupportWindowImageFilter.hxx rename to Modules/PhotoacousticsAlgorithms/ITKFilter/ITKUltrasound/itkSpectra1DSupportWindowImageFilter.hxx diff --git a/Modules/PhotoacousticsAlgorithms/Algorithms/ITKUltrasound/itkSpectra1DSupportWindowToMaskImageFilter.h b/Modules/PhotoacousticsAlgorithms/ITKFilter/ITKUltrasound/itkSpectra1DSupportWindowToMaskImageFilter.h similarity index 100% rename from Modules/PhotoacousticsAlgorithms/Algorithms/ITKUltrasound/itkSpectra1DSupportWindowToMaskImageFilter.h rename to Modules/PhotoacousticsAlgorithms/ITKFilter/ITKUltrasound/itkSpectra1DSupportWindowToMaskImageFilter.h diff --git a/Modules/PhotoacousticsAlgorithms/Algorithms/ITKUltrasound/itkSpectra1DSupportWindowToMaskImageFilter.hxx b/Modules/PhotoacousticsAlgorithms/ITKFilter/ITKUltrasound/itkSpectra1DSupportWindowToMaskImageFilter.hxx similarity index 100% rename from Modules/PhotoacousticsAlgorithms/Algorithms/ITKUltrasound/itkSpectra1DSupportWindowToMaskImageFilter.hxx rename to Modules/PhotoacousticsAlgorithms/ITKFilter/ITKUltrasound/itkSpectra1DSupportWindowToMaskImageFilter.hxx diff --git a/Modules/PhotoacousticsAlgorithms/Algorithms/ITKUltrasound/itkTimeGainCompensationImageFilter.h b/Modules/PhotoacousticsAlgorithms/ITKFilter/ITKUltrasound/itkTimeGainCompensationImageFilter.h similarity index 100% rename from Modules/PhotoacousticsAlgorithms/Algorithms/ITKUltrasound/itkTimeGainCompensationImageFilter.h rename to Modules/PhotoacousticsAlgorithms/ITKFilter/ITKUltrasound/itkTimeGainCompensationImageFilter.h diff --git a/Modules/PhotoacousticsAlgorithms/Algorithms/ITKUltrasound/itkTimeGainCompensationImageFilter.hxx b/Modules/PhotoacousticsAlgorithms/ITKFilter/ITKUltrasound/itkTimeGainCompensationImageFilter.hxx similarity index 100% rename from Modules/PhotoacousticsAlgorithms/Algorithms/ITKUltrasound/itkTimeGainCompensationImageFilter.hxx rename to Modules/PhotoacousticsAlgorithms/ITKFilter/ITKUltrasound/itkTimeGainCompensationImageFilter.hxx diff --git a/Modules/PhotoacousticsAlgorithms/Algorithms/ITKUltrasound/itkVnlFFT1DComplexConjugateToRealImageFilter.h b/Modules/PhotoacousticsAlgorithms/ITKFilter/ITKUltrasound/itkVnlFFT1DComplexConjugateToRealImageFilter.h similarity index 100% rename from Modules/PhotoacousticsAlgorithms/Algorithms/ITKUltrasound/itkVnlFFT1DComplexConjugateToRealImageFilter.h rename to Modules/PhotoacousticsAlgorithms/ITKFilter/ITKUltrasound/itkVnlFFT1DComplexConjugateToRealImageFilter.h diff --git a/Modules/PhotoacousticsAlgorithms/Algorithms/ITKUltrasound/itkVnlFFT1DComplexConjugateToRealImageFilter.hxx b/Modules/PhotoacousticsAlgorithms/ITKFilter/ITKUltrasound/itkVnlFFT1DComplexConjugateToRealImageFilter.hxx similarity index 100% rename from Modules/PhotoacousticsAlgorithms/Algorithms/ITKUltrasound/itkVnlFFT1DComplexConjugateToRealImageFilter.hxx rename to Modules/PhotoacousticsAlgorithms/ITKFilter/ITKUltrasound/itkVnlFFT1DComplexConjugateToRealImageFilter.hxx diff --git a/Modules/PhotoacousticsAlgorithms/Algorithms/ITKUltrasound/itkVnlFFT1DComplexToComplexImageFilter.h b/Modules/PhotoacousticsAlgorithms/ITKFilter/ITKUltrasound/itkVnlFFT1DComplexToComplexImageFilter.h similarity index 100% rename from Modules/PhotoacousticsAlgorithms/Algorithms/ITKUltrasound/itkVnlFFT1DComplexToComplexImageFilter.h rename to Modules/PhotoacousticsAlgorithms/ITKFilter/ITKUltrasound/itkVnlFFT1DComplexToComplexImageFilter.h diff --git a/Modules/PhotoacousticsAlgorithms/Algorithms/ITKUltrasound/itkVnlFFT1DComplexToComplexImageFilter.hxx b/Modules/PhotoacousticsAlgorithms/ITKFilter/ITKUltrasound/itkVnlFFT1DComplexToComplexImageFilter.hxx similarity index 100% rename from Modules/PhotoacousticsAlgorithms/Algorithms/ITKUltrasound/itkVnlFFT1DComplexToComplexImageFilter.hxx rename to Modules/PhotoacousticsAlgorithms/ITKFilter/ITKUltrasound/itkVnlFFT1DComplexToComplexImageFilter.hxx diff --git a/Modules/PhotoacousticsAlgorithms/Algorithms/ITKUltrasound/itkVnlFFT1DRealToComplexConjugateImageFilter.h b/Modules/PhotoacousticsAlgorithms/ITKFilter/ITKUltrasound/itkVnlFFT1DRealToComplexConjugateImageFilter.h similarity index 100% rename from Modules/PhotoacousticsAlgorithms/Algorithms/ITKUltrasound/itkVnlFFT1DRealToComplexConjugateImageFilter.h rename to Modules/PhotoacousticsAlgorithms/ITKFilter/ITKUltrasound/itkVnlFFT1DRealToComplexConjugateImageFilter.h diff --git a/Modules/PhotoacousticsAlgorithms/Algorithms/ITKUltrasound/itkVnlFFT1DRealToComplexConjugateImageFilter.hxx b/Modules/PhotoacousticsAlgorithms/ITKFilter/ITKUltrasound/itkVnlFFT1DRealToComplexConjugateImageFilter.hxx similarity index 100% rename from Modules/PhotoacousticsAlgorithms/Algorithms/ITKUltrasound/itkVnlFFT1DRealToComplexConjugateImageFilter.hxx rename to Modules/PhotoacousticsAlgorithms/ITKFilter/ITKUltrasound/itkVnlFFT1DRealToComplexConjugateImageFilter.hxx diff --git a/Modules/PhotoacousticsAlgorithms/Algorithms/itkPhotoacousticBModeImageFilter.h b/Modules/PhotoacousticsAlgorithms/ITKFilter/itkPhotoacousticBModeImageFilter.h similarity index 100% rename from Modules/PhotoacousticsAlgorithms/Algorithms/itkPhotoacousticBModeImageFilter.h rename to Modules/PhotoacousticsAlgorithms/ITKFilter/itkPhotoacousticBModeImageFilter.h diff --git a/Modules/PhotoacousticsAlgorithms/Algorithms/itkPhotoacousticBModeImageFilter.hxx b/Modules/PhotoacousticsAlgorithms/ITKFilter/itkPhotoacousticBModeImageFilter.hxx similarity index 100% rename from Modules/PhotoacousticsAlgorithms/Algorithms/itkPhotoacousticBModeImageFilter.hxx rename to Modules/PhotoacousticsAlgorithms/ITKFilter/itkPhotoacousticBModeImageFilter.hxx 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 diff --git a/Modules/PhotoacousticsAlgorithms/Resources/DAS.cl b/Modules/PhotoacousticsAlgorithms/Resources/DAS.cl new file mode 100644 index 0000000000..2e8812e0dd --- /dev/null +++ b/Modules/PhotoacousticsAlgorithms/Resources/DAS.cl @@ -0,0 +1,65 @@ +/*=================================================================== + +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 ckDAS( + __global float* dSource, // input image + __global float* dDest, // output buffer + __global unsigned short* usedLines, + __global unsigned int* memoryLocations, + __global unsigned short* AddSamples, + __constant 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 curUsedLines = usedLines[globalPosY * 3 * outputL + 3 * globalPosX]; + unsigned short minLine = usedLines[globalPosY * 3 * outputL + 3 * globalPosX + 1]; + unsigned short maxLine = usedLines[globalPosY * 3 *outputL + 3 * globalPosX + 2]; + + float apod_mult = (float)apodArraySize / (float)curUsedLines; + + unsigned short AddSample = 0; + + float output = 0; + float mult = 0; + + unsigned int MemoryStartAccessPoint = memoryLocations[globalPosY * outputL + globalPosX]; + + for (short l_s = minLine; l_s < maxLine; ++l_s) + { + AddSample = AddSamples[MemoryStartAccessPoint + l_s - minLine]; + if (AddSample < inputS && AddSample >= 0) { + output += apodArray[(int)((l_s - minLine)*apod_mult)] * dSource[(int)(globalPosZ * inputL * inputS + AddSample * inputL + l_s)]; + } + else + --curUsedLines; + } + + dDest[ globalPosZ * outputL * outputS + globalPosY * outputL + globalPosX ] = output / (float)curUsedLines; + } +} \ No newline at end of file diff --git a/Modules/PhotoacousticsAlgorithms/Resources/DASQuadratic.cl b/Modules/PhotoacousticsAlgorithms/Resources/DASQuadratic.cl index 8afda8ff42..12a1456807 100644 --- a/Modules/PhotoacousticsAlgorithms/Resources/DASQuadratic.cl +++ b/Modules/PhotoacousticsAlgorithms/Resources/DASQuadratic.cl @@ -1,76 +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( - __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, + 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); - - 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; + 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); - - short AddSample = 0; - float output = 0; - float delayMultiplicator = pow(1 / (TimeSpacing*SpeedOfSound) * Pitch * TransducerElements / inputL, 2) / s_i / 2; + 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)] * read_imagef( dSource, defaultSampler, (int4)(l_s, AddSample, globalPosZ, 0 )).x; + if (AddSample < inputS && AddSample >= 0) { + 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 similarity index 67% rename from Modules/PhotoacousticsAlgorithms/Resources/DASspherical.cl rename to Modules/PhotoacousticsAlgorithms/Resources/DASSpherical.cl index f9a11f7504..7cedb732a0 100644 --- a/Modules/PhotoacousticsAlgorithms/Resources/DASspherical.cl +++ b/Modules/PhotoacousticsAlgorithms/Resources/DASSpherical.cl @@ -1,79 +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( - __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, + 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); - - 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; + 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); - - short AddSample = 0; - float output = 0; + 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)] * 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/DMAS.cl b/Modules/PhotoacousticsAlgorithms/Resources/DMAS.cl new file mode 100644 index 0000000000..f2ddcd5088 --- /dev/null +++ b/Modules/PhotoacousticsAlgorithms/Resources/DMAS.cl @@ -0,0 +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 ckDMAS( + __global float* dSource, // input image + __global float* dDest, // output buffer + __global unsigned short* usedLines, + __global unsigned int* memoryLocations, + __global unsigned short* AddSamples, + __constant 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 curUsedLines = usedLines[globalPosY * 3 * outputL + 3 * globalPosX]; + unsigned short minLine = usedLines[globalPosY * 3 * outputL + 3 * globalPosX + 1]; + unsigned short maxLine = usedLines[globalPosY * 3 *outputL + 3 * globalPosX + 2]; + + float apod_mult = (float)apodArraySize / (float)curUsedLines; + + unsigned short AddSample1 = 0; + unsigned short AddSample2 = 0; + + float output = 0; + float mult = 0; + + unsigned int MemoryStartAccessPoint = memoryLocations[globalPosY * outputL + globalPosX]; + + for (short l_s1 = minLine; l_s1 < maxLine; ++l_s1) + { + AddSample1 = AddSamples[MemoryStartAccessPoint + l_s1 - minLine]; + if (AddSample1 < inputS && AddSample1 >= 0) { + for (short l_s2 = l_s1 + 1; l_s2 < maxLine; ++l_s2) + { + AddSample2 = AddSamples[MemoryStartAccessPoint + 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/DMASQuadratic.cl b/Modules/PhotoacousticsAlgorithms/Resources/DMASQuadratic.cl index 857b988746..ec7ec8a868 100644 --- a/Modules/PhotoacousticsAlgorithms/Resources/DMASQuadratic.cl +++ b/Modules/PhotoacousticsAlgorithms/Resources/DMASQuadratic.cl @@ -1,92 +1,84 @@ /*=================================================================== 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, + 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); - - 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; + + 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 delayMultiplicator = pow((1 / (TimeSpacing*SpeedOfSound) * Pitch * TransducerElements / inputL), 2) / s_i / 2; - - float mult = 0; - float output = 0; - float AddSample1 = 0; - float AddSample2 = 0; - - for (short l_s1 = minLine; l_s1 < maxLine - 1; ++l_s1) + float apod_mult = apodArraySize / (maxLine - minLine); + + short AddSample1 = 0; + short AddSample2 = 0; + + float output = 0; + float mult = 0; + + float delayMultiplicator = pow(1 / (TimeSpacing*SpeedOfSound) * Pitch * TransducerElements / inputL, 2) / s_i / 2; + + for (short l_s1 = minLine; l_s1 < maxLine; ++l_s1) { AddSample1 = delayMultiplicator * pow((l_s1 - l_i), 2) + s_i + (1-PAImage)*s_i; - if (AddSample1 < inputS && AddSample1 >= 0) - { + 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)]; + if (AddSample1 < inputS && AddSample1 >= 0) { + mult = apodArray[(short)((l_s2 - minLine)*apod_mult)] * dSource[(unsigned int)(globalPosZ * inputL * inputS + AddSample2 * inputL + l_s2)] + * apodArray[(short)((l_s1 - minLine)*apod_mult)] * dSource[(unsigned int)(globalPosZ * inputL * inputS + AddSample1 * inputL + l_s1)]; 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)); + + dDest[ globalPosZ * outputL * outputS + globalPosY * outputL + globalPosX ] = 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 similarity index 52% rename from Modules/PhotoacousticsAlgorithms/Resources/DMASspherical.cl rename to Modules/PhotoacousticsAlgorithms/Resources/DMASSpherical.cl index fc480b0820..a883160ead 100644 --- a/Modules/PhotoacousticsAlgorithms/Resources/DMASspherical.cl +++ b/Modules/PhotoacousticsAlgorithms/Resources/DMASSpherical.cl @@ -1,98 +1,94 @@ /*=================================================================== 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, + 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); - - 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; + + 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 mult = 0; - float output = 0; - float AddSample1 = 0; - float AddSample2 = 0; + float apod_mult = apodArraySize / (maxLine - minLine); + + short AddSample1 = 0; + short AddSample2 = 0; + + float output = 0; + float mult = 0; + + float delayMultiplicator = pow(1 / (TimeSpacing*SpeedOfSound) * Pitch * TransducerElements / inputL, 2) / s_i / 2; - for (short l_s1 = minLine; l_s1 < maxLine - 1; ++l_s1) + for (short l_s1 = minLine; l_s1 < maxLine; ++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) - { + 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( + 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)]; + ) + (1-PAImage)*s_i; + if (AddSample1 < inputS && AddSample1 >= 0) { + mult = apodArray[(short)((l_s2 - minLine)*apod_mult)] * dSource[(unsigned int)(globalPosZ * inputL * inputS + AddSample2 * inputL + l_s2)] + * apodArray[(short)((l_s1 - minLine)*apod_mult)] * dSource[(unsigned int)(globalPosZ * inputL * inputS + AddSample1 * inputL + l_s1)]; 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)); + + dDest[ globalPosZ * outputL * outputS + globalPosY * outputL + globalPosX ] = output / (pow((float)usedLines, 2.0f) - (usedLines - 1)); } -} \ No newline at end of file +} + + \ 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..4f96c2f836 --- /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 * 3 * outputL + 3 * globalPosX]) + { + float l_i = (float)globalPosX / outputL * inputL; + float s_i = (float)globalPosY / outputS * inputS / 2; + + float l_s = usedLines[globalPosY * 3 * outputL + 3 * globalPosX + 1] + 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 * 3 * outputL + 3 * globalPosX]) + { + float l_i = (float)globalPosX / outputL * inputL; + float s_i = (float)globalPosY / outputS * inputS / 2; + + float l_s = usedLines[globalPosY * 3 * outputL + 3 * globalPosX + 1] + 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..4bc9b3bfa3 --- /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[3 * pos]; + } + + sums[id] = sum; + + if (id == (get_global_size(0) - 1) + (get_global_size(1) - 1) * get_global_size(0)) + finalSum[0] = sum + input[3 * id]; + } \ No newline at end of file diff --git a/Modules/PhotoacousticsAlgorithms/Resources/UsedLinesCalculation.cl b/Modules/PhotoacousticsAlgorithms/Resources/UsedLinesCalculation.cl new file mode 100644 index 0000000000..a0d4af6fb7 --- /dev/null +++ b/Modules/PhotoacousticsAlgorithms/Resources/UsedLinesCalculation.cl @@ -0,0 +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 * 3 * outputL + 3 * globalPosX] = (maxLine - minLine); //usedLines + dDest[globalPosY * 3 * outputL + 3 * globalPosX + 1] = minLine; //minLine + dDest[globalPosY * 3 * outputL + 3 * globalPosX + 2] = maxLine; //maxLine + } +} \ No newline at end of file diff --git a/Modules/PhotoacousticsAlgorithms/files.cmake b/Modules/PhotoacousticsAlgorithms/files.cmake index 08f72362af..fc025f27fb 100644 --- a/Modules/PhotoacousticsAlgorithms/files.cmake +++ b/Modules/PhotoacousticsAlgorithms/files.cmake @@ -1,18 +1,23 @@ +file(GLOB_RECURSE H_FILES RELATIVE "${CMAKE_CURRENT_SOURCE_DIR}" "${CMAKE_CURRENT_SOURCE_DIR}/include/*") + set(CPP_FILES - mitkPhotoacousticImage.cpp - - Algorithms/mitkPhotoacousticBeamformingFilter.cpp - - Algorithms/OCL/mitkPhotoacousticOCLBeamformer.cpp - - Algorithms/OCL/mitkPhotoacousticBModeFilter.cpp + source/mitkPhotoacousticImage.cpp + source/mitkPhotoacousticBeamformingFilter.cpp + source/OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.cpp + source/OpenCLFilter/mitkPhotoacousticBModeFilter.cpp + source/OpenCLFilter/mitkPhotoacousticOCLUsedLinesCalculation.cpp + source/OpenCLFilter/mitkPhotoacousticOCLDelayCalculation.cpp + source/OpenCLFilter/mitkPhotoacousticOCLMemoryLocSum.cpp ) set(RESOURCE_FILES DASQuadratic.cl - DMASQuadratic.cl - DASspherical.cl - DMASspherical.cl + DASSpherical.cl BModeAbs.cl BModeAbsLog.cl -) \ No newline at end of file + UsedLinesCalculation.cl + MemoryLocSum.cl + DelayCalculation.cl + DMAS.cl + DAS.cl +) diff --git a/Modules/PhotoacousticsAlgorithms/include/OpenCLFilter/mitkPhotoacousticBModeFilter.h b/Modules/PhotoacousticsAlgorithms/include/OpenCLFilter/mitkPhotoacousticBModeFilter.h new file mode 100644 index 0000000000..af464ed599 --- /dev/null +++ b/Modules/PhotoacousticsAlgorithms/include/OpenCLFilter/mitkPhotoacousticBModeFilter.h @@ -0,0 +1,130 @@ +/*=================================================================== + +The Medical Imaging Interaction Toolkit (MITK) + +Copyright (c) German Cancer Research Center, +Division of Medical and Biological Informatics. +All rights reserved. + +This software is distributed WITHOUT ANY WARRANTY; without +even the implied warranty of MERCHANTABILITY or FITNESS FOR +A PARTICULAR PURPOSE. + +See LICENSE.txt or http://www.mitk.org for details. + +===================================================================*/ + +#ifndef _MITKPHOTOACOUSTICSBMODEFILTER_H_ +#define _MITKPHOTOACOUSTICSBMODEFILTER_H_ + +#ifdef PHOTOACOUSTICS_USE_GPU +#include "mitkOclDataSetToDataSetFilter.h" +#endif + +#include +#include "mitkImageToImageFilter.h" + +namespace mitk +{ + #ifdef PHOTOACOUSTICS_USE_GPU + + /** Documentation + * + * \brief The PhotoacousticOCLBModeFilter simply takes the absolute of all pixels in the image. + * + * The filter gives the option to use a log after taking the absolute. + */ + + class PhotoacousticOCLBModeFilter : public OclDataSetToDataSetFilter, public itk::Object + { + + public: + mitkClassMacroItkParent(PhotoacousticOCLBModeFilter, itk::Object); + itkNewMacro(Self); + + /** + * @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: + + PhotoacousticOCLBModeFilter(); + + virtual ~PhotoacousticOCLBModeFilter(); + + 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 diff --git a/Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticOCLBeamformer.h b/Modules/PhotoacousticsAlgorithms/include/OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.h similarity index 53% rename from Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticOCLBeamformer.h rename to Modules/PhotoacousticsAlgorithms/include/OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.h index 9a7c90e0de..060ee5c0e3 100644 --- a/Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticOCLBeamformer.h +++ b/Modules/PhotoacousticsAlgorithms/include/OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.h @@ -1,129 +1,140 @@ /*=================================================================== 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" +#include "mitkOclDataSetToDataSetFilter.h" #include +#include "mitkPhotoacousticOCLDelayCalculation.h" +#include "mitkPhotoacousticOCLMemoryLocSum.h" +#include "mitkPhotoacousticOCLUsedLinesCalculation.h" + +#include "mitkPhotoacousticBeamformingSettings.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 PhotoacousticOCLBeamformingFilter : public OclDataSetToDataSetFilter, public itk::Object { public: - mitkClassMacroItkParent(PhotoacousticOCLBeamformer, itk::Object); + mitkClassMacroItkParent(PhotoacousticOCLBeamformingFilter, 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 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); - mitk::Image::Pointer GetOutput(); + void* GetOutput(); + + /** + * @brief GetOutputAsImage Returns an mitk::Image constructed from the processed data + */ + 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]; - } + /** 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 }; - - 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) + void SetConfig(BeamformingSettings settings) { - m_SpeedOfSound = SpeedOfSound; - m_TimeSpacing = timeSpacing; - m_Pitch = Pitch; - m_Angle = Angle; - m_PAImage = PAImage; - m_TransducerElements = transducerElements; + m_ConfOld = m_Conf; + m_Conf = settings; } protected: /** Constructor */ - PhotoacousticOCLBeamformer(); + PhotoacousticOCLBeamformingFilter(); /** Destructor */ - virtual ~PhotoacousticOCLBeamformer(); + virtual ~PhotoacousticOCLBeamformingFilter(); /** Initialize the filter */ bool Initialize(); + void UpdateDataBuffers(); + 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]; + 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; + + BeamformingSettings m_Conf; + BeamformingSettings m_ConfOld; + + 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; + + cl_mem m_ApodizationBuffer; + cl_mem m_MemoryLocationsBuffer; + cl_mem m_DelaysBuffer; + cl_mem m_UsedLinesBuffer; + + std::chrono::steady_clock::time_point m_Begin; }; } -#endif \ No newline at end of file +#endif +#endif diff --git a/Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticBModeFilter.h b/Modules/PhotoacousticsAlgorithms/include/OpenCLFilter/mitkPhotoacousticOCLDelayCalculation.h similarity index 62% copy from Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticBModeFilter.h copy to Modules/PhotoacousticsAlgorithms/include/OpenCLFilter/mitkPhotoacousticOCLDelayCalculation.h index adc24614df..38bee14c68 100644 --- a/Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticBModeFilter.h +++ b/Modules/PhotoacousticsAlgorithms/include/OpenCLFilter/mitkPhotoacousticOCLDelayCalculation.h @@ -1,91 +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 _MITKPHOTOACOUSTICSOCLBEAMFORMER_H_ -#define _MITKPHOTOACOUSTICSOCLBEAMFORMER_H_ +#ifndef _MITKPHOTOACOUSTICSDELAYCALC_H_ +#define _MITKPHOTOACOUSTICSDELAYCALC_H_ -#include "mitkOclImageToImageFilter.h" +#ifdef PHOTOACOUSTICS_USE_GPU + +#include "mitkOclDataSetToDataSetFilter.h" #include +#include "mitkPhotoacousticBeamformingSettings.h" 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 PhotoacousticBModeFilter : public OclImageToImageFilter, public itk::Object + class OCLDelayCalculation : public OclDataSetToDataSetFilter, public itk::Object { public: - mitkClassMacroItkParent(PhotoacousticBModeFilter, 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. */ - void SetInput(Image::Pointer image); /** Update the filter */ void Update(); - void SetParameters(bool useLogFilter) + void SetConfig(BeamformingSettings conf) { - m_UseLogFilter = useLogFilter; + 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 */ - PhotoacousticBModeFilter(); + OCLDelayCalculation(); /** Destructor */ - virtual ~PhotoacousticBModeFilter(); + virtual ~OCLDelayCalculation(); /** Initialize the filter */ bool Initialize(); void Execute(); mitk::PixelType GetOutputType() { - return mitk::MakeScalarPixelType(); + return mitk::MakeScalarPixelType(); } int GetBytesPerElem() { - return sizeof(float); + return sizeof(unsigned short); } virtual us::Module* GetModule(); + int m_sizeThis; private: /** The OpenCL kernel for the filter */ cl_kernel m_PixelCalculation; - bool m_UseLogFilter; + + BeamformingSettings m_Conf; + 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/mitkPhotoacousticBModeFilter.h b/Modules/PhotoacousticsAlgorithms/include/OpenCLFilter/mitkPhotoacousticOCLMemoryLocSum.h similarity index 69% copy from Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticBModeFilter.h copy to Modules/PhotoacousticsAlgorithms/include/OpenCLFilter/mitkPhotoacousticOCLMemoryLocSum.h index adc24614df..bce56edbc2 100644 --- a/Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticBModeFilter.h +++ b/Modules/PhotoacousticsAlgorithms/include/OpenCLFilter/mitkPhotoacousticOCLMemoryLocSum.h @@ -1,91 +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 _MITKPHOTOACOUSTICSOCLBEAMFORMER_H_ -#define _MITKPHOTOACOUSTICSOCLBEAMFORMER_H_ +#ifndef _MITKPHOTOACOUSTICSPARALLELSUM_H_ +#define _MITKPHOTOACOUSTICSPARALLELSUM_H_ -#include "mitkOclImageToImageFilter.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 PhotoacousticBModeFilter : public OclImageToImageFilter, public itk::Object + class OCLMemoryLocSum : public OclDataSetToDataSetFilter, public itk::Object { public: - mitkClassMacroItkParent(PhotoacousticBModeFilter, 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. */ - void SetInput(Image::Pointer image); /** Update the filter */ void Update(); - void SetParameters(bool useLogFilter) + void SetInputDimensions(unsigned int* dimensions) { - m_UseLogFilter = useLogFilter; + m_Dim[0] = dimensions[0]; + m_Dim[1] = dimensions[1]; + } + + unsigned int GetSum() + { + return m_Sum; } protected: /** Constructor */ - PhotoacousticBModeFilter(); + OCLMemoryLocSum(); /** Destructor */ - virtual ~PhotoacousticBModeFilter(); + virtual ~OCLMemoryLocSum(); /** Initialize the filter */ bool Initialize(); void Execute(); mitk::PixelType GetOutputType() { - return mitk::MakeScalarPixelType(); + return mitk::MakeScalarPixelType(); } int GetBytesPerElem() { - return sizeof(float); + return sizeof(unsigned short); } virtual us::Module* GetModule(); + int m_sizeThis; private: /** The OpenCL kernel for the filter */ cl_kernel m_PixelCalculation; - bool m_UseLogFilter; + + unsigned int m_Sum; + unsigned int m_Dim[3]; }; } +#endif #endif \ No newline at end of file diff --git a/Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticBModeFilter.h b/Modules/PhotoacousticsAlgorithms/include/OpenCLFilter/mitkPhotoacousticOCLUsedLinesCalculation.h similarity index 69% rename from Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticBModeFilter.h rename to Modules/PhotoacousticsAlgorithms/include/OpenCLFilter/mitkPhotoacousticOCLUsedLinesCalculation.h index adc24614df..1203feca2a 100644 --- a/Modules/PhotoacousticsAlgorithms/Algorithms/OCL/mitkPhotoacousticBModeFilter.h +++ b/Modules/PhotoacousticsAlgorithms/include/OpenCLFilter/mitkPhotoacousticOCLUsedLinesCalculation.h @@ -1,91 +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 _MITKPHOTOACOUSTICSOCLBEAMFORMER_H_ -#define _MITKPHOTOACOUSTICSOCLBEAMFORMER_H_ +#ifndef _MITKPHOTOACOUSTICSOCLUSEDLINESCALCULATION_H_ +#define _MITKPHOTOACOUSTICSOCLUSEDLINESCALCULATION_H_ -#include "mitkOclImageToImageFilter.h" +#ifdef PHOTOACOUSTICS_USE_GPU + +#include "mitkOclDataSetToDataSetFilter.h" #include +#include "mitkPhotoacousticBeamformingSettings.h" 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 PhotoacousticBModeFilter : public OclImageToImageFilter, public itk::Object + class OCLUsedLinesCalculation : public OclDataSetToDataSetFilter, public itk::Object { public: - mitkClassMacroItkParent(PhotoacousticBModeFilter, itk::Object); + 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. */ - void SetInput(Image::Pointer image); /** Update the filter */ void Update(); - void SetParameters(bool useLogFilter) + void SetConfig(BeamformingSettings conf) { - m_UseLogFilter = useLogFilter; + m_Conf = conf; } protected: /** Constructor */ - PhotoacousticBModeFilter(); + OCLUsedLinesCalculation(); /** Destructor */ - virtual ~PhotoacousticBModeFilter(); + virtual ~OCLUsedLinesCalculation(); /** Initialize the filter */ bool Initialize(); void Execute(); mitk::PixelType GetOutputType() { - return mitk::MakeScalarPixelType(); + return mitk::MakeScalarPixelType(); } int GetBytesPerElem() { - return sizeof(float); + return sizeof(unsigned short); } virtual us::Module* GetModule(); + int m_sizeThis; + private: /** The OpenCL kernel for the filter */ cl_kernel m_PixelCalculation; - bool m_UseLogFilter; + + BeamformingSettings m_Conf; + float m_part; }; } -#endif \ No newline at end of file +#endif +#endif diff --git a/Modules/PhotoacousticsAlgorithms/Algorithms/mitkPhotoacousticBeamformingFilter.h b/Modules/PhotoacousticsAlgorithms/include/mitkPhotoacousticBeamformingFilter.h similarity index 71% rename from Modules/PhotoacousticsAlgorithms/Algorithms/mitkPhotoacousticBeamformingFilter.h rename to Modules/PhotoacousticsAlgorithms/include/mitkPhotoacousticBeamformingFilter.h index e2654fc7fd..f4956b6859 100644 --- a/Modules/PhotoacousticsAlgorithms/Algorithms/mitkPhotoacousticBeamformingFilter.h +++ b/Modules/PhotoacousticsAlgorithms/include/mitkPhotoacousticBeamformingFilter.h @@ -1,108 +1,80 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ - #ifndef MITK_PHOTOACOUSTICS_BEAMFORMING_FILTER #define MITK_PHOTOACOUSTICS_BEAMFORMING_FILTER #include "mitkImageToImageFilter.h" #include +#include "./OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.h" +#include "mitkPhotoacousticBeamformingSettings.h" namespace mitk { - //##Documentation //## @brief //## @ingroup Process + class BeamformingFilter : public ImageToImageFilter { public: mitkClassMacro(BeamformingFilter, ImageToImageFilter); itkFactorylessNewMacro(Self) itkCloneMacro(Self) - struct beamformingSettings - { - float Pitch = 0.0003; // [m] - float SpeedOfSound = 1540; // [m/s] - unsigned int SamplesPerLine = 2048; - unsigned int ReconstructionLines = 128; - float RecordTime = 0.00006; // [s] - float TimeSpacing = 0.0000000000001; // [s] - unsigned int TransducerElements = 128; - bool partial = false; - unsigned int CropBounds[2] = { 0,0 }; - - bool UseGPU = true; - - enum DelayCalc {QuadApprox, Spherical}; - DelayCalc DelayCalculationMethod = QuadApprox; - - enum Apodization {Hamm, Hann, Box}; - Apodization Apod = Hann; - - enum BeamformingAlgorithm {DMAS, DAS}; - BeamformingAlgorithm Algorithm = DAS; - - float Angle = 10; - bool Photoacoustic = true; - float BPHighPass = 50; - float BPLowPass = 50; - bool UseBP = false; - }; - - void Configure(beamformingSettings settings); + void Configure(BeamformingSettings settings); void SetProgressHandle(std::function progressHandle); protected: BeamformingFilter(); ~BeamformingFilter(); virtual void GenerateInputRequestedRegion() override; virtual void GenerateOutputInformation() override; virtual void GenerateData() override; //##Description //## @brief Time when Header was last initialized itk::TimeStamp m_TimeOfHeaderInitialization; std::function m_ProgressHandle; float* VonHannFunction(int samples); float* HammFunction(int samples); float* BoxFunction(int samples); void DASQuadraticLine(float* input, float* output, float inputDim[2], float outputDim[2], const short& line, float* apodisation, const short& apodArraySize); void DASSphericalLine(float* input, float* output, float inputDim[2], float outputDim[2], const short& line, float* apodisation, const short& apodArraySize); void DMASQuadraticLine(float* input, float* output, float inputDim[2], float outputDim[2], const short& line, float* apodisation, const short& apodArraySize); void DMASSphericalLine(float* input, float* output, float inputDim[2], float outputDim[2], const short& line, float* apodisation, const short& apodArraySize); float* m_OutputData; float* m_InputData; float* m_InputDataPuffer; - beamformingSettings m_Conf; - }; + BeamformingSettings m_Conf; + mitk::PhotoacousticOCLBeamformingFilter::Pointer m_BeamformingOclFilter; + }; } // namespace mitk #endif //MITK_PHOTOACOUSTICS_BEAMFORMING_FILTER diff --git a/Modules/PhotoacousticsAlgorithms/include/mitkPhotoacousticBeamformingSettings.h b/Modules/PhotoacousticsAlgorithms/include/mitkPhotoacousticBeamformingSettings.h new file mode 100644 index 0000000000..117f2919c8 --- /dev/null +++ b/Modules/PhotoacousticsAlgorithms/include/mitkPhotoacousticBeamformingSettings.h @@ -0,0 +1,76 @@ +/*=================================================================== + +The Medical Imaging Interaction Toolkit (MITK) + +Copyright (c) German Cancer Research Center, +Division of Medical and Biological Informatics. +All rights reserved. + +This software is distributed WITHOUT ANY WARRANTY; without +even the implied warranty of MERCHANTABILITY or FITNESS FOR +A PARTICULAR PURPOSE. + +See LICENSE.txt or http://www.mitk.org for details. + +===================================================================*/ + +#ifndef MITK_PHOTOACOUSTICS_BEAMFORMING_SETTINGS +#define MITK_PHOTOACOUSTICS_BEAMFORMING_SETTINGS + + + +namespace mitk { + //##Documentation + //## @brief + //## @ingroup Process + + class BeamformingSettings + { + public: + float Pitch = 0.0003; // [m] + float SpeedOfSound = 1540; // [m/s] + unsigned int SamplesPerLine = 2048; + unsigned int ReconstructionLines = 128; + float RecordTime = 0.00006; // [s] + float TimeSpacing = 0.0000000000001; // [s] + unsigned short TransducerElements = 128; + bool partial = false; + unsigned int CropBounds[2] = { 0,0 }; + unsigned int inputDim[3] = { 1,1,1 }; + unsigned int upperCutoff = 0; + + bool UseGPU = true; + + enum DelayCalc { QuadApprox, Spherical }; + DelayCalc DelayCalculationMethod = QuadApprox; + + enum Apodization { Hamm, Hann, Box }; + Apodization Apod = Hann; + + enum BeamformingAlgorithm { DMAS, DAS }; + BeamformingAlgorithm Algorithm = DAS; + + float Angle = 10; + bool isPhotoacousticImage = true; + float BPHighPass = 50; + float BPLowPass = 50; + bool UseBP = false; + + //this method ignores changes in BPLow/BPHigh/cropBounds/Algorithm/some more, as those are insignifiant in all current situations + static bool SettingsChangedOpenCL(const BeamformingSettings& lhs, const BeamformingSettings& rhs) + { + return !(((lhs.Angle - rhs.Angle) < 0.001f) && + (lhs.Apod == rhs.Apod) && + (lhs.DelayCalculationMethod == rhs.DelayCalculationMethod) && + (lhs.isPhotoacousticImage == rhs.isPhotoacousticImage) && + ((lhs.Pitch - rhs.Pitch) < 0.0000001f) && + (lhs.ReconstructionLines == rhs.ReconstructionLines) && + (lhs.RecordTime - rhs.RecordTime) < 0.00000001f && + (lhs.SamplesPerLine == rhs.SamplesPerLine) && + ((lhs.SpeedOfSound - rhs.SpeedOfSound) < 0.01f) && + ((lhs.TimeSpacing - rhs.TimeSpacing) < 0.0000000001f) && + (lhs.TransducerElements == rhs.TransducerElements)); + } + }; +} +#endif \ No newline at end of file diff --git a/Modules/PhotoacousticsAlgorithms/mitkPhotoacousticImage.h b/Modules/PhotoacousticsAlgorithms/include/mitkPhotoacousticImage.h similarity index 84% copy from Modules/PhotoacousticsAlgorithms/mitkPhotoacousticImage.h copy to Modules/PhotoacousticsAlgorithms/include/mitkPhotoacousticImage.h index daab4a43cf..9f00543a22 100644 --- a/Modules/PhotoacousticsAlgorithms/mitkPhotoacousticImage.h +++ b/Modules/PhotoacousticsAlgorithms/include/mitkPhotoacousticImage.h @@ -1,53 +1,53 @@ /*=================================================================== 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 mitkPhotoacousticImage_H_HEADER_INCLUDED #define mitkPhotoacousticImage_H_HEADER_INCLUDED #include "itkObject.h" #include "mitkCommon.h" #include "mitkImage.h" #include -#include "mitkOclImageToImageFilter.h" +#include "mitkPhotoacousticBeamformingSettings.h" #include "mitkPhotoacousticBeamformingFilter.h" - #include "MitkPhotoacousticsAlgorithmsExports.h" namespace mitk { class MITKPHOTOACOUSTICSALGORITHMS_EXPORT PhotoacousticImage : public itk::Object { public: mitkClassMacroItkParent(mitk::PhotoacousticImage, itk::Object); itkFactorylessNewMacro(Self); enum BModeMethod { ShapeDetection, Abs }; - mitk::Image::Pointer ApplyBmodeFilter(mitk::Image::Pointer inputImage, BModeMethod method = BModeMethod::Abs, bool UseLogFilter = false, float resampleSpacing = 0.15); + mitk::Image::Pointer ApplyBmodeFilter(mitk::Image::Pointer inputImage, BModeMethod method = BModeMethod::Abs, bool UseGPU = false, bool UseLogFilter = false, float resampleSpacing = 0.15); // mitk::Image::Pointer ApplyScatteringCompensation(mitk::Image::Pointer inputImage, int scatteringCoefficient); mitk::Image::Pointer ApplyResampling(mitk::Image::Pointer inputImage, unsigned int outputSize[2]); - mitk::Image::Pointer ApplyBeamforming(mitk::Image::Pointer inputImage, BeamformingFilter::beamformingSettings config, int cutoff, std::function progressHandle = [](int, std::string) {}); + mitk::Image::Pointer ApplyBeamforming(mitk::Image::Pointer inputImage, BeamformingSettings config, std::function progressHandle = [](int, std::string) {}); mitk::Image::Pointer ApplyCropping(mitk::Image::Pointer inputImage, int above, int below, int right, int left, int minSlice, int maxSlice); mitk::Image::Pointer BandpassFilter(mitk::Image::Pointer data, float recordTime, float BPHighPass, float BPLowPass, float alpha); protected: PhotoacousticImage(); virtual ~PhotoacousticImage(); + mitk::BeamformingFilter::Pointer m_BeamformingFilter; + itk::Image::Pointer BPFunction(mitk::Image::Pointer reference, int cutoffFrequencyPixelHighPass, int cutoffFrequencyPixelLowPass, float alpha); }; } // namespace mitk #endif /* mitkPhotoacousticImage_H_HEADER_INCLUDED */ diff --git a/Modules/PhotoacousticsAlgorithms/mitkPhotoacousticImage.cpp b/Modules/PhotoacousticsAlgorithms/mitkPhotoacousticImage.cpp index 9ba451e9cd..08f92552bc 100644 --- a/Modules/PhotoacousticsAlgorithms/mitkPhotoacousticImage.cpp +++ b/Modules/PhotoacousticsAlgorithms/mitkPhotoacousticImage.cpp @@ -1,641 +1,522 @@ /*=================================================================== 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 "ITKFilter/ITKUltrasound/itkBModeImageFilter.h" +#include "ITKFilter/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" +#include "ITKFilter/ITKUltrasound/itkFFT1DComplexConjugateToRealImageFilter.h" +#include "ITKFilter/ITKUltrasound/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 UseLogFilter, float resampleSpacing) +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) { - auto input = ApplyCropping(inputImage, 0, 0, 0, 0, 0, inputImage->GetDimension(2) - 1); - PhotoacousticBModeFilter::Pointer filter = PhotoacousticBModeFilter::New(); - filter->SetParameters(UseLogFilter); - filter->SetInput(input); - filter->Update(); + 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(); - auto out = filter->GetOutput(); + out = filter->GetOutput(); - if(resampleSpacing == 0) - return out; + 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) +mitk::Image::Pointer mitk::PhotoacousticImage::ApplyBeamforming(mitk::Image::Pointer inputImage, BeamformingSettings config, int cutoff, std::function progressHandle) { - // crop the image - // set the Maximum Size of the image to 4096 - unsigned short lowerCutoff = 0; - if (((unsigned int)(4096 + cutoff)) < inputImage->GetDimension(1)) - { - lowerCutoff = (unsigned short)(inputImage->GetDimension(1) - 4096); - } - - config.RecordTime = config.RecordTime - (cutoff + lowerCutoff) / inputImage->GetDimension(1) * config.RecordTime; // adjust the recorded time lost by cropping + 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, lowerCutoff, 0, 0, config.CropBounds[0], config.CropBounds[1]); - - // resample the image in horizontal direction - if (inputImage->GetDimension(0) != config.ReconstructionLines) - { - progressHandle(0, "resampling image"); - auto begin = std::chrono::high_resolution_clock::now(); - unsigned int dim[2] = { config.ReconstructionLines, processedImage->GetDimension(1) }; - processedImage = ApplyResampling(processedImage, dim); - auto end = std::chrono::high_resolution_clock::now(); - MITK_DEBUG << "Upsampling from " << inputImage->GetDimension(0) << " to " << config.ReconstructionLines << " lines completed in " << ((double)std::chrono::duration_cast(end - begin).count()) / 1000000 << "ms" << std::endl; - } + Image::Pointer processedImage = ApplyCropping(inputImage, cutoff, 0, 0, 0, config.CropBounds[0], config.CropBounds[1]); + config.inputDim[0] = processedImage->GetDimension(0); + config.inputDim[1] = processedImage->GetDimension(1); + config.inputDim[2] = processedImage->GetDimension(2); // 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; -} - - -/* -mitk::CropFilter::CropFilter() - : m_PixelCalculation(NULL) -{ - this->AddSourceFile("CropFilter.cl"); - - this->m_FilterID = "CropFilter"; -} - -mitk::CropFilter::~CropFilter() -{ - if (this->m_PixelCalculation) - { - clReleaseKernel(m_PixelCalculation); - } -} - -void mitk::CropFilter::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::CropFilter::Execute() -{ - cl_int clErr = 0; - - try - { - this->InitExec(this->m_PixelCalculation); - } - catch (const mitk::Exception& e) - { - MITK_ERROR << "Catched exception while initializing filter: " << e.what(); - return; - } - - us::ServiceReference ref = GetModuleContext()->GetServiceReference(); - OclResourceService* resources = GetModuleContext()->GetService(ref); - cl_context gpuContext = resources->GetContext(); - - unsigned int size = m_OutputDim[0] * m_OutputDim[1] * m_OutputDim[2] * m_Images; - - cl_mem cl_input = clCreateBuffer(gpuContext, CL_MEM_READ_ONLY | CL_MEM_USE_HOST_PTR, sizeof(float) * size, m_Data, &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_Images)); - - 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::CropFilter::GetModule() -{ - return us::GetModuleContext()->GetModule(); -} - -bool mitk::CropFilter::Initialize() -{ - bool buildErr = true; - cl_int clErr = 0; - - if (OclFilter::Initialize()) - { - this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "ckCrop", &clErr); - buildErr |= CHECK_OCL_ERR(clErr); - } - return (OclFilter::IsInitialized() && buildErr); -} - -void mitk::CropFilter::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); -} - -mitk::Image::Pointer mitk::CropFilter::GetOutput() -{ - if (m_Output->IsModified(GPU_DATA)) - { - void* pData = m_Output->TransferDataToCPU(m_CommandQue); - - const unsigned int dimension = 3; - unsigned int* dimensions = m_OutputDim; - - const mitk::SlicedGeometry3D::Pointer p_slg = m_Input->GetMITKImage()->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); - } - - MITK_DEBUG << "Image Initialized."; - - return m_Output->GetMITKImage(); -} - -*/ +} \ No newline at end of file diff --git a/Modules/PhotoacousticsAlgorithms/mitkPhotoacousticImage.h b/Modules/PhotoacousticsAlgorithms/mitkPhotoacousticImage.h index daab4a43cf..2200124e8d 100644 --- a/Modules/PhotoacousticsAlgorithms/mitkPhotoacousticImage.h +++ b/Modules/PhotoacousticsAlgorithms/mitkPhotoacousticImage.h @@ -1,53 +1,50 @@ /*=================================================================== 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 mitkPhotoacousticImage_H_HEADER_INCLUDED #define mitkPhotoacousticImage_H_HEADER_INCLUDED #include "itkObject.h" #include "mitkCommon.h" #include "mitkImage.h" #include -#include "mitkOclImageToImageFilter.h" #include "mitkPhotoacousticBeamformingFilter.h" - #include "MitkPhotoacousticsAlgorithmsExports.h" namespace mitk { class MITKPHOTOACOUSTICSALGORITHMS_EXPORT PhotoacousticImage : public itk::Object { public: mitkClassMacroItkParent(mitk::PhotoacousticImage, itk::Object); itkFactorylessNewMacro(Self); enum BModeMethod { ShapeDetection, Abs }; - mitk::Image::Pointer ApplyBmodeFilter(mitk::Image::Pointer inputImage, BModeMethod method = BModeMethod::Abs, bool UseLogFilter = false, float resampleSpacing = 0.15); + mitk::Image::Pointer ApplyBmodeFilter(mitk::Image::Pointer inputImage, BModeMethod method = BModeMethod::Abs, bool UseGPU = false, bool UseLogFilter = false, float resampleSpacing = 0.15); // mitk::Image::Pointer ApplyScatteringCompensation(mitk::Image::Pointer inputImage, int scatteringCoefficient); mitk::Image::Pointer ApplyResampling(mitk::Image::Pointer inputImage, unsigned int outputSize[2]); - mitk::Image::Pointer ApplyBeamforming(mitk::Image::Pointer inputImage, BeamformingFilter::beamformingSettings config, int cutoff, std::function progressHandle = [](int, std::string) {}); + mitk::Image::Pointer ApplyBeamforming(mitk::Image::Pointer inputImage, BeamformingSettings config, int cutoff, std::function progressHandle = [](int, std::string) {}); mitk::Image::Pointer ApplyCropping(mitk::Image::Pointer inputImage, int above, int below, int right, int left, int minSlice, int maxSlice); mitk::Image::Pointer BandpassFilter(mitk::Image::Pointer data, float recordTime, float BPHighPass, float BPLowPass, float alpha); protected: PhotoacousticImage(); virtual ~PhotoacousticImage(); itk::Image::Pointer BPFunction(mitk::Image::Pointer reference, int cutoffFrequencyPixelHighPass, int cutoffFrequencyPixelLowPass, float alpha); }; } // namespace mitk #endif /* mitkPhotoacousticImage_H_HEADER_INCLUDED */ diff --git a/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticBModeFilter.cpp b/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticBModeFilter.cpp new file mode 100644 index 0000000000..d19bd50198 --- /dev/null +++ b/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticBModeFilter.cpp @@ -0,0 +1,221 @@ +/*=================================================================== + +The Medical Imaging Interaction Toolkit (MITK) + +Copyright (c) German Cancer Research Center, +Division of Medical and Biological Informatics. +All rights reserved. + +This software is distributed WITHOUT ANY WARRANTY; without +even the implied warranty of MERCHANTABILITY or FITNESS FOR +A PARTICULAR PURPOSE. + +See LICENSE.txt or http://www.mitk.org for details. + +===================================================================*/ + +#include "./OpenCLFilter/mitkPhotoacousticBModeFilter.h" +#include "usServiceReference.h" +#include + +#ifdef PHOTOACOUSTICS_USE_GPU + +mitk::PhotoacousticOCLBModeFilter::PhotoacousticOCLBModeFilter() + : m_PixelCalculation(NULL) +{ + this->AddSourceFile("BModeAbs.cl"); + this->AddSourceFile("BModeAbsLog.cl"); + + this->m_FilterID = "BModeFilter"; + + this->Initialize(); +} + +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 + { + size_t outputSize = m_InputDim[0] * m_InputDim[1] * m_InputDim[2]; + this->InitExec(this->m_PixelCalculation, m_InputDim, outputSize, 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) +{ + 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)) + { + 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); + } + + 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/source/OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.cpp b/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.cpp new file mode 100644 index 0000000000..d37f4222bc --- /dev/null +++ b/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.cpp @@ -0,0 +1,261 @@ +/*=================================================================== + +The Medical Imaging Interaction Toolkit (MITK) + +Copyright (c) German Cancer Research Center, +Division of Medical and Biological Informatics. +All rights reserved. + +This software is distributed WITHOUT ANY WARRANTY; without +even the implied warranty of MERCHANTABILITY or FITNESS FOR +A PARTICULAR PURPOSE. + +See LICENSE.txt or http://www.mitk.org for details. + +===================================================================*/ + +#ifdef PHOTOACOUSTICS_USE_GPU + +#include "./OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.h" +#include "usServiceReference.h" + +mitk::PhotoacousticOCLBeamformingFilter::PhotoacousticOCLBeamformingFilter() +: m_PixelCalculation( NULL ), m_InputImage(mitk::Image::New()), m_ApodizationBuffer(nullptr), m_MemoryLocationsBuffer(nullptr), m_DelaysBuffer(nullptr), m_UsedLinesBuffer(nullptr) +{ + this->AddSourceFile("DASQuadratic.cl"); + this->AddSourceFile("DAS.cl"); + this->AddSourceFile("DMAS.cl"); + this->m_FilterID = "OpenCLBeamformingFilter"; + + this->Initialize(); + + 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); + + 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(); +} + +mitk::PhotoacousticOCLBeamformingFilter::~PhotoacousticOCLBeamformingFilter() +{ + if ( this->m_PixelCalculation ) + { + clReleaseKernel( m_PixelCalculation ); + } + + if (m_ApodizationBuffer) clReleaseMemObject(m_ApodizationBuffer); +} + +void mitk::PhotoacousticOCLBeamformingFilter::Update() +{ + //Check if context & program available + if (!this->Initialize()) + { + us::ServiceReference ref = GetModuleContext()->GetServiceReference(); + OclResourceService* resources = GetModuleContext()->GetService(ref); + + // clean-up also the resources + resources->InvalidateStorage(); + mitkThrow() <<"Filter is not initialized. Cannot update."; + } + else{ + // Execute + this->Execute(); + } +} + +void mitk::PhotoacousticOCLBeamformingFilter::UpdateDataBuffers() +{ + if (BeamformingSettings::SettingsChangedOpenCL(m_Conf, m_ConfOld)) + { + cl_int clErr = 0; + MITK_INFO << "Updating Data for new configuration"; + //Initialize the Output + try + { + size_t outputSize = m_Conf.ReconstructionLines * m_Conf.SamplesPerLine * m_Conf.inputDim[2]; + m_OutputDim[0] = m_Conf.ReconstructionLines; + m_OutputDim[1] = m_Conf.SamplesPerLine; + m_OutputDim[2] = m_Conf.inputDim[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; + } + + // create the apodisation buffer + + if (m_Apodisation == nullptr) + { + MITK_INFO << "No apodisation function set; Beamforming will be done without any apodisation."; + m_Apodisation = new float[1]; + m_Apodisation[0] = 1; + m_ApodArraySize = 1; + } + + us::ServiceReference ref = GetModuleContext()->GetServiceReference(); + OclResourceService* resources = GetModuleContext()->GetService(ref); + cl_context gpuContext = resources->GetContext(); + + if (m_ApodizationBuffer) clReleaseMemObject(m_ApodizationBuffer); + + this->m_ApodizationBuffer = clCreateBuffer(gpuContext, CL_MEM_READ_ONLY | CL_MEM_USE_HOST_PTR, sizeof(float) * m_ApodArraySize, m_Apodisation, &clErr); + CHECK_OCL_ERR(clErr); + + // calculate used lines + + m_UsedLinesCalculation->SetConfig(m_Conf); + m_UsedLinesCalculation->Update(); + + // calculate memory locations of the Delays + + 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(); + + m_UsedLinesBuffer = m_UsedLinesCalculation->GetGPUOutput()->GetGPUBuffer(); + m_MemoryLocationsBuffer = m_SumFilter->GetGPUOutput()->GetGPUBuffer(); + + // calculate the Delays + + m_DelayCalculation->SetConfig(m_Conf); + m_DelayCalculation->SetInputs(m_UsedLinesBuffer, m_MemoryLocationsBuffer, size); + m_DelayCalculation->Update(); + + m_DelaysBuffer = m_DelayCalculation->GetGPUOutput()->GetGPUBuffer(); + + m_ConfOld = m_Conf; + } +} + +void mitk::PhotoacousticOCLBeamformingFilter::Execute() +{ + cl_int clErr = 0; + UpdateDataBuffers(); + + clErr = clSetKernelArg(this->m_PixelCalculation, 2, sizeof(cl_mem), &(this->m_UsedLinesBuffer)); + clErr |= clSetKernelArg(this->m_PixelCalculation, 3, sizeof(cl_mem), &(this->m_MemoryLocationsBuffer)); + clErr |= clSetKernelArg(this->m_PixelCalculation, 4, sizeof(cl_mem), &(this->m_DelaysBuffer)); + clErr |= clSetKernelArg(this->m_PixelCalculation, 5, sizeof(cl_mem), &(this->m_ApodizationBuffer)); + clErr |= clSetKernelArg(this->m_PixelCalculation, 6, sizeof(cl_ushort), &(this->m_ApodArraySize)); + clErr |= clSetKernelArg(this->m_PixelCalculation, 7, sizeof(cl_uint), &(this->m_Conf.inputDim[0])); + clErr |= clSetKernelArg(this->m_PixelCalculation, 8, sizeof(cl_uint), &(this->m_Conf.inputDim[1])); + clErr |= clSetKernelArg(this->m_PixelCalculation, 9, sizeof(cl_uint), &(this->m_Conf.inputDim[2])); + clErr |= clSetKernelArg(this->m_PixelCalculation, 10, sizeof(cl_uint), &(this->m_Conf.ReconstructionLines)); + clErr |= clSetKernelArg(this->m_PixelCalculation, 11, sizeof(cl_uint), &(this->m_Conf.SamplesPerLine)); + + CHECK_OCL_ERR(clErr); + + // execute the filter on a 3D NDRange + if (m_OutputDim[2] == 1 || m_ChunkSize[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 BeamformingSettings::BeamformingAlgorithm::DAS: + { + this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "ckDAS", &clErr); + break; + } + case BeamformingSettings::BeamformingAlgorithm::DMAS: + { + this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "ckDMAS", &clErr); + break; + } + default: + { + MITK_INFO << "No beamforming algorithm specified, setting to DAS"; + this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "ckDAS", &clErr); + break; + } + } + buildErr |= CHECK_OCL_ERR( clErr ); + } + return (OclFilter::IsInitialized() && buildErr ); +} + +void mitk::PhotoacousticOCLBeamformingFilter::SetInput(mitk::Image::Pointer image) +{ + OclDataSetToDataSetFilter::SetInput(image); + + m_InputImage = image; + m_Conf.inputDim[0] = m_InputImage->GetDimension(0); + m_Conf.inputDim[1] = m_InputImage->GetDimension(1); + m_Conf.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_Conf.inputDim[0] = dimensions[0]; + m_Conf.inputDim[1] = dimensions[1]; + m_Conf.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->SetImportVolume(pData, 0, 0, mitk::Image::ReferenceMemory); + } + + MITK_DEBUG << "Image Initialized."; + + return outputImage; +} + +void* mitk::PhotoacousticOCLBeamformingFilter::GetOutput() +{ + return OclDataSetToDataSetFilter::GetOutput(); +} +#endif diff --git a/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLDelayCalculation.cpp b/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLDelayCalculation.cpp new file mode 100644 index 0000000000..13b216d8d8 --- /dev/null +++ b/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLDelayCalculation.cpp @@ -0,0 +1,119 @@ +/*=================================================================== + +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 "./OpenCLFilter/mitkPhotoacousticOCLDelayCalculation.h" +#include "usServiceReference.h" +#include "mitkImageReadAccessor.h" + +mitk::OCLDelayCalculation::OCLDelayCalculation() + : m_PixelCalculation(NULL) +{ + this->AddSourceFile("DelayCalculation.cl"); + this->m_FilterID = "DelayCalculation"; + + this->Initialize(); +} + +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 == 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 == BeamformingSettings::DelayCalc::Spherical) + m_DelayMultiplicatorRaw = 1 / (m_Conf.TimeSpacing*m_Conf.SpeedOfSound) * (m_Conf.Pitch*m_Conf.TransducerElements); + + m_IsPAImage = m_Conf.isPhotoacousticImage; + + 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 == BeamformingSettings::DelayCalc::QuadApprox) + this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "ckDelayCalculationQuad", &clErr); + if (m_Conf.DelayCalculationMethod == 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/source/OpenCLFilter/mitkPhotoacousticOCLMemoryLocSum.cpp b/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLMemoryLocSum.cpp new file mode 100644 index 0000000000..ebfac4ec3f --- /dev/null +++ b/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLMemoryLocSum.cpp @@ -0,0 +1,119 @@ +/*=================================================================== + +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 "./OpenCLFilter/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"; + + this->Initialize(); +} + +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/source/OpenCLFilter/mitkPhotoacousticOCLUsedLinesCalculation.cpp b/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLUsedLinesCalculation.cpp new file mode 100644 index 0000000000..7f5f1eee1b --- /dev/null +++ b/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLUsedLinesCalculation.cpp @@ -0,0 +1,108 @@ +/*=================================================================== + +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 +#define _USE_MATH_DEFINES + +#include +#include "./OpenCLFilter/mitkPhotoacousticOCLUsedLinesCalculation.h" +#include "usServiceReference.h" +#include "mitkImageReadAccessor.h" + +mitk::OCLUsedLinesCalculation::OCLUsedLinesCalculation() + : m_PixelCalculation(NULL) +{ + this->AddSourceFile("UsedLinesCalculation.cl"); + this->m_FilterID = "UsedLinesCalculation"; + + this->Initialize(); +} + +mitk::OCLUsedLinesCalculation::~OCLUsedLinesCalculation() +{ + if (this->m_PixelCalculation) + { + clReleaseKernel(m_PixelCalculation); + } +} + +void mitk::OCLUsedLinesCalculation::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::OCLUsedLinesCalculation::Execute() +{ + cl_int clErr = 0; + + unsigned int gridDim[3] = { m_Conf.ReconstructionLines, m_Conf.SamplesPerLine, 1 }; + size_t outputSize = gridDim[0] * gridDim[1] * 3; + + try + { + this->InitExecNoInput(this->m_PixelCalculation, gridDim, outputSize, sizeof(unsigned short)); + } + catch (const mitk::Exception& e) + { + MITK_ERROR << "Caught exception while initializing filter: " << e.what(); + return; + } + + m_part = (tan(m_Conf.Angle / 360 * 2 * M_PI) * ((m_Conf.SpeedOfSound * m_Conf.TimeSpacing)) / (m_Conf.Pitch * m_Conf.TransducerElements)) * m_Conf.inputDim[0]; + + clErr = clSetKernelArg(this->m_PixelCalculation, 1, sizeof(cl_float), &(this->m_part)); + clErr |= clSetKernelArg(this->m_PixelCalculation, 2, sizeof(cl_uint), &(this->m_Conf.inputDim[0])); + clErr |= clSetKernelArg(this->m_PixelCalculation, 3, sizeof(cl_uint), &(this->m_Conf.inputDim[1])); + + 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); +} + +us::Module *mitk::OCLUsedLinesCalculation::GetModule() +{ + return us::GetModuleContext()->GetModule(); +} + +bool mitk::OCLUsedLinesCalculation::Initialize() +{ + bool buildErr = true; + cl_int clErr = 0; + + if (OclFilter::Initialize()) + { + this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "ckUsedLines", &clErr); + buildErr |= CHECK_OCL_ERR(clErr); + } + return (OclFilter::IsInitialized() && buildErr); +} +#endif diff --git a/Modules/PhotoacousticsAlgorithms/Algorithms/mitkPhotoacousticBeamformingFilter.cpp b/Modules/PhotoacousticsAlgorithms/source/mitkPhotoacousticBeamformingFilter.cpp similarity index 67% rename from Modules/PhotoacousticsAlgorithms/Algorithms/mitkPhotoacousticBeamformingFilter.cpp rename to Modules/PhotoacousticsAlgorithms/source/mitkPhotoacousticBeamformingFilter.cpp index 47587fb719..9edb570cc2 100644 --- a/Modules/PhotoacousticsAlgorithms/Algorithms/mitkPhotoacousticBeamformingFilter.cpp +++ b/Modules/PhotoacousticsAlgorithms/source/mitkPhotoacousticBeamformingFilter.cpp @@ -1,607 +1,541 @@ /*=================================================================== 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 "mitkPhotoacousticBeamformingFilter.h" mitk::BeamformingFilter::BeamformingFilter() : m_OutputData(nullptr), m_InputData(nullptr) { this->SetNumberOfIndexedInputs(1); this->SetNumberOfRequiredInputs(1); m_ProgressHandle = [](int, std::string) {}; + + m_BeamformingOclFilter = mitk::PhotoacousticOCLBeamformingFilter::New(); } 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 + 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; + 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) { + int progInterval = output->GetDimension(2) / 20 > 1 ? output->GetDimension(2) / 20 : 1; + // the interval at which we update the gui progress bar + + float inputDim[2] = { (float)input->GetDimension(0), (float)input->GetDimension(1) }; + float outputDim[2] = { (float)output->GetDimension(0), (float)output->GetDimension(1) }; + for (unsigned int i = 0; i < output->GetDimension(2); ++i) // seperate Slices should get Beamforming seperately applied { mitk::ImageReadAccessor inputReadAccessor(input, input->GetSliceData(i)); - m_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.Algorithm == BeamformingSettings::BeamformingAlgorithm::DAS) { - if (m_Conf.DelayCalculationMethod == beamformingSettings::DelayCalc::QuadApprox) + 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) + 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) + else if (m_Conf.Algorithm == BeamformingSettings::BeamformingAlgorithm::DMAS) { - if (m_Conf.DelayCalculationMethod == beamformingSettings::DelayCalc::QuadApprox) + 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) + 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; } } + #ifdef PHOTOACOUSTICS_USE_GPU else { - mitk::PhotoacousticOCLBeamformer::Pointer m_oclFilter = mitk::PhotoacousticOCLBeamformer::New(); - try { - if (m_Conf.Algorithm == beamformingSettings::BeamformingAlgorithm::DAS) + // first, we check whether the data is float, other formats are unsupported + if (!(input->GetPixelType().GetTypeAsString() == "scalar (float)" || input->GetPixelType().GetTypeAsString() == " (float)")) { - 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); + MITK_ERROR << "Pixel type is not float, abort"; + return; } - 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); - - if (chunkSize < oclOutputDim[2]) - { - bool skip = false; - for (unsigned int i = 0; !skip && 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)) - { - // 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 - } - else if ((oclOutputDim[2]) - (i * chunkSize) >= chunkSize) - chunk->Initialize(input->GetPixelType(), 3, oclInputDimChunk); - else - { - chunk->Initialize(input->GetPixelType(), 3, oclInputDimLastChunk); - m_oclFilter->SetOutputDim(oclOutputDimLastChunk); - } - chunk->SetSpacing(input->GetGeometry()->GetSpacing()); + m_ProgressHandle(50, "performing reconstruction"); - 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_BeamformingOclFilter->SetApodisation(ApodWindow, apodArraySize); + m_BeamformingOclFilter->SetConfig(m_Conf); + m_BeamformingOclFilter->SetInput(input); + m_BeamformingOclFilter->Update(); - m_oclFilter->SetInput(chunk); - m_oclFilter->Update(); - auto out = m_oclFilter->GetOutput(); - - for (unsigned int s = i * chunkSize; s < oclOutputDim[2]; ++s) // TODO: make the bounds here smaller... - { - mitk::ImageReadAccessor copy(out, out->GetSliceData(s - i * chunkSize)); - output->SetImportSlice(const_cast(copy.GetData()), 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); - } + void* out = m_BeamformingOclFilter->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) +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 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; + AddSample = delayMultiplicator * pow((l_s - l_i), 2) + s_i + (1 - m_Conf.isPhotoacousticImage)*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 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; + ) + (1 - m_Conf.isPhotoacousticImage)*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 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; + AddSample[l_s] = (short)(delayMultiplicator * pow((minLine + l_s - l_i), 2) + s_i) + (1 - m_Conf.isPhotoacousticImage)*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)); + output[sample*(short)outputL + line] = 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 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; + ) + (1 - m_Conf.isPhotoacousticImage)*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)); + output[sample*(short)outputL + line] = output[sample*(short)outputL + line] / (float)(pow(usedLines, 2) - (usedLines - 1)); delete[] AddSample; } } diff --git a/Modules/PhotoacousticsAlgorithms/mitkPhotoacousticImage.cpp b/Modules/PhotoacousticsAlgorithms/source/mitkPhotoacousticImage.cpp similarity index 75% copy from Modules/PhotoacousticsAlgorithms/mitkPhotoacousticImage.cpp copy to Modules/PhotoacousticsAlgorithms/source/mitkPhotoacousticImage.cpp index 9ba451e9cd..a838dad9fc 100644 --- a/Modules/PhotoacousticsAlgorithms/mitkPhotoacousticImage.cpp +++ b/Modules/PhotoacousticsAlgorithms/source/mitkPhotoacousticImage.cpp @@ -1,641 +1,521 @@ /*=================================================================== 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 "../ITKFilter/ITKUltrasound/itkBModeImageFilter.h" +#include "../ITKFilter/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" +#include "../ITKFilter/ITKUltrasound/itkFFT1DComplexConjugateToRealImageFilter.h" +#include "../ITKFilter/ITKUltrasound/itkFFT1DRealToComplexConjugateImageFilter.h" -mitk::PhotoacousticImage::PhotoacousticImage() +mitk::PhotoacousticImage::PhotoacousticImage() : m_BeamformingFilter(BeamformingFilter::New()) { 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 UseLogFilter, float resampleSpacing) +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) { - auto input = ApplyCropping(inputImage, 0, 0, 0, 0, 0, inputImage->GetDimension(2) - 1); - PhotoacousticBModeFilter::Pointer filter = PhotoacousticBModeFilter::New(); - filter->SetParameters(UseLogFilter); - filter->SetInput(input); - filter->Update(); + 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(); - auto out = filter->GetOutput(); + out = filter->GetOutput(); - if(resampleSpacing == 0) - return out; + 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) +mitk::Image::Pointer mitk::PhotoacousticImage::ApplyBeamforming(mitk::Image::Pointer inputImage, BeamformingSettings config, std::function progressHandle) { - // crop the image - // set the Maximum Size of the image to 4096 - unsigned short lowerCutoff = 0; - if (((unsigned int)(4096 + cutoff)) < inputImage->GetDimension(1)) - { - lowerCutoff = (unsigned short)(inputImage->GetDimension(1) - 4096); - } - - config.RecordTime = config.RecordTime - (cutoff + lowerCutoff) / inputImage->GetDimension(1) * config.RecordTime; // adjust the recorded time lost by cropping + config.RecordTime = config.RecordTime - (float)(config.upperCutoff) / (float)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, lowerCutoff, 0, 0, config.CropBounds[0], config.CropBounds[1]); - - // resample the image in horizontal direction - if (inputImage->GetDimension(0) != config.ReconstructionLines) - { - progressHandle(0, "resampling image"); - auto begin = std::chrono::high_resolution_clock::now(); - unsigned int dim[2] = { config.ReconstructionLines, processedImage->GetDimension(1) }; - processedImage = ApplyResampling(processedImage, dim); - auto end = std::chrono::high_resolution_clock::now(); - MITK_DEBUG << "Upsampling from " << inputImage->GetDimension(0) << " to " << config.ReconstructionLines << " lines completed in " << ((double)std::chrono::duration_cast(end - begin).count()) / 1000000 << "ms" << std::endl; - } + Image::Pointer processedImage = ApplyCropping(inputImage, config.upperCutoff, 0, 0, 0, config.CropBounds[0], config.CropBounds[1]); + config.inputDim[0] = processedImage->GetDimension(0); + config.inputDim[1] = processedImage->GetDimension(1); + config.inputDim[2] = processedImage->GetDimension(2); // perform the beamforming - BeamformingFilter::Pointer Beamformer = BeamformingFilter::New(); - Beamformer->SetInput(processedImage); - Beamformer->Configure(config); - Beamformer->SetProgressHandle(progressHandle); - Beamformer->UpdateLargestPossibleRegion(); + m_BeamformingFilter->SetInput(processedImage); + m_BeamformingFilter->Configure(config); + m_BeamformingFilter->SetProgressHandle(progressHandle); + m_BeamformingFilter->UpdateLargestPossibleRegion(); - processedImage = Beamformer->GetOutput(); + processedImage = m_BeamformingFilter->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; -} - - -/* -mitk::CropFilter::CropFilter() - : m_PixelCalculation(NULL) -{ - this->AddSourceFile("CropFilter.cl"); - - this->m_FilterID = "CropFilter"; -} - -mitk::CropFilter::~CropFilter() -{ - if (this->m_PixelCalculation) - { - clReleaseKernel(m_PixelCalculation); - } -} - -void mitk::CropFilter::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::CropFilter::Execute() -{ - cl_int clErr = 0; - - try - { - this->InitExec(this->m_PixelCalculation); - } - catch (const mitk::Exception& e) - { - MITK_ERROR << "Catched exception while initializing filter: " << e.what(); - return; - } - - us::ServiceReference ref = GetModuleContext()->GetServiceReference(); - OclResourceService* resources = GetModuleContext()->GetService(ref); - cl_context gpuContext = resources->GetContext(); - - unsigned int size = m_OutputDim[0] * m_OutputDim[1] * m_OutputDim[2] * m_Images; - - cl_mem cl_input = clCreateBuffer(gpuContext, CL_MEM_READ_ONLY | CL_MEM_USE_HOST_PTR, sizeof(float) * size, m_Data, &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_Images)); - - 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::CropFilter::GetModule() -{ - return us::GetModuleContext()->GetModule(); -} - -bool mitk::CropFilter::Initialize() -{ - bool buildErr = true; - cl_int clErr = 0; - - if (OclFilter::Initialize()) - { - this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "ckCrop", &clErr); - buildErr |= CHECK_OCL_ERR(clErr); - } - return (OclFilter::IsInitialized() && buildErr); -} - -void mitk::CropFilter::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); -} - -mitk::Image::Pointer mitk::CropFilter::GetOutput() -{ - if (m_Output->IsModified(GPU_DATA)) - { - void* pData = m_Output->TransferDataToCPU(m_CommandQue); - - const unsigned int dimension = 3; - unsigned int* dimensions = m_OutputDim; - - const mitk::SlicedGeometry3D::Pointer p_slg = m_Input->GetMITKImage()->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); - } - - MITK_DEBUG << "Image Initialized."; - - return m_Output->GetMITKImage(); -} - -*/ +} \ No newline at end of file diff --git a/Plugins/org.mitk.gui.qt.photoacoustics.imageprocessing/CMakeLists.txt b/Plugins/org.mitk.gui.qt.photoacoustics.imageprocessing/CMakeLists.txt index ce750d05d7..0b44d554e6 100644 --- a/Plugins/org.mitk.gui.qt.photoacoustics.imageprocessing/CMakeLists.txt +++ b/Plugins/org.mitk.gui.qt.photoacoustics.imageprocessing/CMakeLists.txt @@ -1,7 +1,11 @@ project(org_mitk_gui_qt_photoacoustics_imageprocessing) mitk_create_plugin( EXPORT_DIRECTIVE IMAGEPROCESSING_EXPORT EXPORTED_INCLUDE_SUFFIXES src MODULE_DEPENDS MitkQtWidgetsExt MitkPhotoacousticsAlgorithms ) + +IF(MITK_USE_OpenCL) + add_definitions(-DPHOTOACOUSTICS_USE_GPU) +ENDIF(MITK_USE_OpenCL) \ No newline at end of file diff --git a/Plugins/org.mitk.gui.qt.photoacoustics.imageprocessing/src/internal/PAImageProcessing.cpp b/Plugins/org.mitk.gui.qt.photoacoustics.imageprocessing/src/internal/PAImageProcessing.cpp index 8c59dee114..33d9af950b 100644 --- a/Plugins/org.mitk.gui.qt.photoacoustics.imageprocessing/src/internal/PAImageProcessing.cpp +++ b/Plugins/org.mitk.gui.qt.photoacoustics.imageprocessing/src/internal/PAImageProcessing.cpp @@ -1,1080 +1,1115 @@ /*=================================================================== 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. ===================================================================*/ // Blueberry #include #include // Qmitk #include "PAImageProcessing.h" // Qt #include #include #include #include //mitk image #include #include "mitkPhotoacousticImage.h" #include "mitkPhotoacousticBeamformingFilter.h" //other #include #include #include const std::string PAImageProcessing::VIEW_ID = "org.mitk.views.paimageprocessing"; -PAImageProcessing::PAImageProcessing() : m_ResampleSpacing(0), m_UseLogfilter(false) +PAImageProcessing::PAImageProcessing() : m_ResampleSpacing(0), m_UseLogfilter(false), m_FilterBank(mitk::PhotoacousticImage::New()) { qRegisterMetaType(); qRegisterMetaType(); } void PAImageProcessing::SetFocus() { m_Controls.buttonApplyBModeFilter->setFocus(); } void PAImageProcessing::CreateQtPartControl(QWidget *parent) { // create GUI widgets from the Qt Designer's .ui file m_Controls.setupUi(parent); connect(m_Controls.buttonApplyBModeFilter, SIGNAL(clicked()), this, SLOT(StartBmodeThread())); connect(m_Controls.DoResampling, SIGNAL(clicked()), this, SLOT(UseResampling())); connect(m_Controls.Logfilter, SIGNAL(clicked()), this, SLOT(UseLogfilter())); connect(m_Controls.ResamplingValue, SIGNAL(valueChanged(double)), this, SLOT(SetResampling())); connect(m_Controls.buttonApplyBeamforming, SIGNAL(clicked()), this, SLOT(StartBeamformingThread())); connect(m_Controls.buttonApplyCropFilter, SIGNAL(clicked()), this, SLOT(StartCropThread())); connect(m_Controls.buttonApplyBandpass, SIGNAL(clicked()), this, SLOT(StartBandpassThread())); connect(m_Controls.UseImageSpacing, SIGNAL(clicked()), this, SLOT(UseImageSpacing())); connect(m_Controls.ScanDepth, SIGNAL(valueChanged(double)), this, SLOT(UpdateImageInfo())); connect(m_Controls.SpeedOfSound, SIGNAL(valueChanged(double)), this, SLOT(UpdateImageInfo())); + connect(m_Controls.SpeedOfSound, SIGNAL(valueChanged(double)), this, SLOT(ChangedSOSBeamforming())); + connect(m_Controls.BPSpeedOfSound, SIGNAL(valueChanged(double)), this, SLOT(ChangedSOSBandpass())); connect(m_Controls.Samples, SIGNAL(valueChanged(int)), this, SLOT(UpdateImageInfo())); connect(m_Controls.UseImageSpacing, SIGNAL(clicked()), this, SLOT(UpdateImageInfo())); connect(m_Controls.boundLow, SIGNAL(valueChanged(int)), this, SLOT(UpdateBounds())); connect(m_Controls.boundHigh, SIGNAL(valueChanged(int)), this, SLOT(UpdateBounds())); connect(m_Controls.Partial, SIGNAL(clicked()), this, SLOT(UpdateBounds())); connect(m_Controls.BatchProcessing, SIGNAL(clicked()), this, SLOT(BatchProcessing())); connect(m_Controls.StepBeamforming, SIGNAL(clicked()), this, SLOT(UpdateSaveBoxes())); connect(m_Controls.StepCropping, SIGNAL(clicked()), this, SLOT(UpdateSaveBoxes())); connect(m_Controls.StepBandpass, SIGNAL(clicked()), this, SLOT(UpdateSaveBoxes())); connect(m_Controls.StepBMode, SIGNAL(clicked()), this, SLOT(UpdateSaveBoxes())); UpdateSaveBoxes(); m_Controls.DoResampling->setChecked(false); m_Controls.ResamplingValue->setEnabled(false); m_Controls.progressBar->setMinimum(0); m_Controls.progressBar->setMaximum(100); m_Controls.progressBar->setVisible(false); m_Controls.UseImageSpacing->setToolTip("Image spacing of y-Axis must be in us, x-Axis in mm."); m_Controls.UseImageSpacing->setToolTipDuration(5000); m_Controls.ProgressInfo->setVisible(false); m_Controls.UseBP->hide(); + #ifndef PHOTOACOUSTICS_USE_GPU + m_Controls.UseGPUBf->setEnabled(false); + m_Controls.UseGPUBf->setChecked(false); + m_Controls.UseGPUBmode->setEnabled(false); + m_Controls.UseGPUBmode->setChecked(false); + #endif + UseImageSpacing(); } +void PAImageProcessing::ChangedSOSBandpass() +{ + m_Controls.SpeedOfSound->setValue(m_Controls.BPSpeedOfSound->value()); +} + +void PAImageProcessing::ChangedSOSBeamforming() +{ + m_Controls.BPSpeedOfSound->setValue(m_Controls.SpeedOfSound->value()); +} + std::vector splitpath( const std::string& str , const std::set delimiters) { std::vector result; char const* pch = str.c_str(); char const* start = pch; for (; *pch; ++pch) { if (delimiters.find(*pch) != delimiters.end()) { if (start != pch) { std::string str(start, pch); result.push_back(str); } else { result.push_back(""); } start = pch + 1; } } result.push_back(start); return result; } void PAImageProcessing::UpdateSaveBoxes() { if (m_Controls.StepBeamforming->isChecked()) m_Controls.SaveBeamforming->setEnabled(true); else m_Controls.SaveBeamforming->setEnabled(false); if (m_Controls.StepCropping->isChecked()) m_Controls.SaveCropping->setEnabled(true); else m_Controls.SaveCropping->setEnabled(false); if (m_Controls.StepBandpass->isChecked()) m_Controls.SaveBandpass->setEnabled(true); else m_Controls.SaveBandpass->setEnabled(false); if (m_Controls.StepBMode->isChecked()) m_Controls.SaveBMode->setEnabled(true); else m_Controls.SaveBMode->setEnabled(false); } void PAImageProcessing::BatchProcessing() { QFileDialog LoadDialog(nullptr, "Select Files to be processed"); LoadDialog.setFileMode(QFileDialog::FileMode::ExistingFiles); LoadDialog.setNameFilter(tr("Images (*.nrrd)")); LoadDialog.setViewMode(QFileDialog::Detail); QStringList fileNames; if (LoadDialog.exec()) fileNames = LoadDialog.selectedFiles(); QString saveDir = QFileDialog::getExistingDirectory(nullptr, tr("Select Directory To Save To"), "", QFileDialog::ShowDirsOnly | QFileDialog::DontResolveSymlinks); DisableControls(); std::set delims{'/'}; - mitk::PhotoacousticImage::Pointer filterbank = mitk::PhotoacousticImage::New(); - bool doSteps[] = { m_Controls.StepBeamforming->isChecked(), m_Controls.StepCropping->isChecked() , m_Controls.StepBandpass->isChecked(), m_Controls.StepBMode->isChecked() }; bool saveSteps[] = { m_Controls.SaveBeamforming->isChecked(), m_Controls.SaveCropping->isChecked() , m_Controls.SaveBandpass->isChecked(), m_Controls.SaveBMode->isChecked() }; for (int fileNumber = 0; fileNumber < fileNames.size(); ++fileNumber) { m_Controls.progressBar->setValue(0); m_Controls.progressBar->setVisible(true); m_Controls.ProgressInfo->setVisible(true); m_Controls.ProgressInfo->setText("loading file"); QString filename = fileNames.at(fileNumber); auto split = splitpath(filename.toStdString(), delims); std::string imageName = split.at(split.size()-1); // remove ".nrrd" imageName = imageName.substr(0, imageName.size()-5); mitk::Image::Pointer image = mitk::IOUtil::LoadImage(filename.toStdString().c_str()); UpdateBFSettings(image); // Beamforming if (doSteps[0]) { std::function progressHandle = [this](int progress, std::string progressInfo) { this->UpdateProgress(progress, progressInfo); }; m_Controls.progressBar->setValue(100); - image = filterbank->ApplyBeamforming(image, BFconfig, m_Controls.Cutoff->value(), progressHandle); + image = m_FilterBank->ApplyBeamforming(image, BFconfig, progressHandle); if (saveSteps[0]) { std::string saveFileName = saveDir.toStdString() + "/" + imageName + " beamformed" + ".nrrd"; mitk::IOUtil::Save(image, saveFileName); } } // Cropping if (doSteps[1]) { m_Controls.ProgressInfo->setText("cropping image"); - image = filterbank->ApplyCropping(image, m_Controls.CutoffAbove->value(), m_Controls.CutoffBelow->value(), 0, 0, 0, image->GetDimension(2) - 1); + image = m_FilterBank->ApplyCropping(image, m_Controls.CutoffAbove->value(), m_Controls.CutoffBelow->value(), 0, 0, 0, image->GetDimension(2) - 1); if (saveSteps[1]) { std::string saveFileName = saveDir.toStdString() + "/" + imageName + " cropped" + ".nrrd"; mitk::IOUtil::Save(image, saveFileName); } } // Bandpass if (doSteps[2]) { m_Controls.ProgressInfo->setText("applying bandpass"); float recordTime = image->GetDimension(1)*image->GetGeometry()->GetSpacing()[1] / 1000 / m_Controls.BPSpeedOfSound->value(); // add a safeguard so the program does not chrash when applying a Bandpass that reaches out of the bounds of the image float maxFrequency = 1 / (recordTime / image->GetDimension(1)) * image->GetDimension(1) / 2 / 2 / 1000; float BPHighPass = 1000000 * m_Controls.BPhigh->value(); // [Hz] float BPLowPass = maxFrequency - 1000000 * m_Controls.BPlow->value(); // [Hz] if (BPLowPass > maxFrequency && m_Controls.UseBP->isChecked()) { QMessageBox Msgbox; Msgbox.setText("LowPass too low, disabled it."); Msgbox.exec(); BPLowPass = 0; } if (BPLowPass < 0 && m_Controls.UseBP->isChecked()) { QMessageBox Msgbox; Msgbox.setText("LowPass too high, disabled it."); Msgbox.exec(); BPLowPass = 0; } if (BPHighPass > maxFrequency && m_Controls.UseBP->isChecked()) { QMessageBox Msgbox; Msgbox.setText("HighPass too high, disabled it."); Msgbox.exec(); BPHighPass = 0; } if (BPHighPass > maxFrequency - BFconfig.BPLowPass) { QMessageBox Msgbox; Msgbox.setText("HighPass higher than LowPass, disabled both."); Msgbox.exec(); BPHighPass = 0; BPLowPass = 0; } - image = filterbank->BandpassFilter(image, recordTime, BPHighPass, BPLowPass, m_Controls.BPFalloff->value()); + image = m_FilterBank->BandpassFilter(image, recordTime, BPHighPass, BPLowPass, m_Controls.BPFalloff->value()); if (saveSteps[2]) { std::string saveFileName = saveDir.toStdString() + "/" + imageName + " bandpassed" + ".nrrd"; mitk::IOUtil::Save(image, saveFileName); } } // Bmode if (doSteps[3]) { m_Controls.ProgressInfo->setText("applying bmode filter"); - - if (m_Controls.BModeMethod->currentText() == "Simple Abs Filter") - image = filterbank->ApplyBmodeFilter(image, mitk::PhotoacousticImage::BModeMethod::Abs, m_UseLogfilter, m_ResampleSpacing); - else if (m_Controls.BModeMethod->currentText() == "Shape Detection") - image = filterbank->ApplyBmodeFilter(image, mitk::PhotoacousticImage::BModeMethod::ShapeDetection, m_UseLogfilter, m_ResampleSpacing); + bool useGPU = m_Controls.UseGPUBmode->isChecked(); + + if (m_Controls.BModeMethod->currentText() == "Absolute Filter") + image = m_FilterBank->ApplyBmodeFilter(image, mitk::PhotoacousticImage::BModeMethod::Abs, useGPU, m_UseLogfilter, m_ResampleSpacing); + else if (m_Controls.BModeMethod->currentText() == "Envelope Detection") + image = m_FilterBank->ApplyBmodeFilter(image, mitk::PhotoacousticImage::BModeMethod::ShapeDetection, useGPU, m_UseLogfilter, m_ResampleSpacing); + if (saveSteps[3]) { std::string saveFileName = saveDir.toStdString() + "/" + imageName + " bmode" + ".nrrd"; mitk::IOUtil::Save(image, saveFileName); } } m_Controls.progressBar->setVisible(false); } EnableControls(); } void PAImageProcessing::StartBeamformingThread() { QList nodes = this->GetDataManagerSelection(); if (nodes.empty()) return; mitk::DataStorage::Pointer storage = this->GetDataStorage(); mitk::DataNode::Pointer node = nodes.front(); if (!node) { // Nothing selected. Inform the user and return QMessageBox::information(NULL, "Template", "Please load and select an image before starting image processing."); return; } mitk::BaseData* data = node->GetData(); if (data) { // test if this data item is an image or not (could also be a surface or something totally different) mitk::Image* image = dynamic_cast(data); if (image) { UpdateBFSettings(image); std::stringstream message; std::string name; message << "Performing beamforming for image "; if (node->GetName(name)) { // a property called "name" was found for this DataNode message << "'" << name << "'"; m_OldNodeName = name; } else m_OldNodeName = " "; message << "."; MITK_INFO << message.str(); m_Controls.progressBar->setValue(0); m_Controls.progressBar->setVisible(true); m_Controls.ProgressInfo->setVisible(true); m_Controls.ProgressInfo->setText("started"); m_Controls.buttonApplyBeamforming->setText("working..."); DisableControls(); BeamformingThread *thread = new BeamformingThread(); connect(thread, &BeamformingThread::result, this, &PAImageProcessing::HandleBeamformingResults); connect(thread, &BeamformingThread::updateProgress, this, &PAImageProcessing::UpdateProgress); connect(thread, &BeamformingThread::finished, thread, &QObject::deleteLater); thread->setConfig(BFconfig); - thread->setCutoff(m_Controls.Cutoff->value()); thread->setInputImage(image); + thread->setFilterBank(m_FilterBank); MITK_INFO << "Started new thread for Beamforming"; thread->start(); } } } void PAImageProcessing::HandleBeamformingResults(mitk::Image::Pointer image) { auto newNode = mitk::DataNode::New(); newNode->SetData(image); // name the new Data node std::stringstream newNodeName; newNodeName << m_OldNodeName << " "; - if (BFconfig.Algorithm == mitk::BeamformingFilter::beamformingSettings::BeamformingAlgorithm::DAS) + if (BFconfig.Algorithm == mitk::BeamformingSettings::BeamformingAlgorithm::DAS) newNodeName << "DAS bf, "; - else if (BFconfig.Algorithm == mitk::BeamformingFilter::beamformingSettings::BeamformingAlgorithm::DMAS) + else if (BFconfig.Algorithm == mitk::BeamformingSettings::BeamformingAlgorithm::DMAS) newNodeName << "DMAS bf, "; - if (BFconfig.DelayCalculationMethod == mitk::BeamformingFilter::beamformingSettings::DelayCalc::QuadApprox) + if (BFconfig.DelayCalculationMethod == mitk::BeamformingSettings::DelayCalc::QuadApprox) newNodeName << "q. delay"; - if (BFconfig.DelayCalculationMethod == mitk::BeamformingFilter::beamformingSettings::DelayCalc::Spherical) + if (BFconfig.DelayCalculationMethod == mitk::BeamformingSettings::DelayCalc::Spherical) newNodeName << "s. delay"; newNode->SetName(newNodeName.str()); // update level window for the current dynamic range mitk::LevelWindow levelWindow; newNode->GetLevelWindow(levelWindow); auto data = newNode->GetData(); levelWindow.SetAuto(dynamic_cast(data), true, true); newNode->SetLevelWindow(levelWindow); // add new node to data storage this->GetDataStorage()->Add(newNode); // disable progress bar m_Controls.progressBar->setVisible(false); m_Controls.ProgressInfo->setVisible(false); m_Controls.buttonApplyBeamforming->setText("Apply Beamforming"); EnableControls(); // update rendering mitk::RenderingManager::GetInstance()->InitializeViews( dynamic_cast(data)->GetGeometry(), mitk::RenderingManager::REQUEST_UPDATE_ALL, true); mitk::RenderingManager::GetInstance()->RequestUpdateAll(); } void PAImageProcessing::StartBmodeThread() { QList nodes = this->GetDataManagerSelection(); if (nodes.empty()) return; mitk::DataStorage::Pointer storage = this->GetDataStorage(); mitk::DataNode::Pointer node = nodes.front(); if (!node) { // Nothing selected. Inform the user and return QMessageBox::information(NULL, "Template", "Please load and select an image before starting image processing."); return; } mitk::BaseData* data = node->GetData(); if (data) { // test if this data item is an image or not (could also be a surface or something totally different) mitk::Image* image = dynamic_cast(data); if (image) { UpdateBFSettings(image); std::stringstream message; std::string name; message << "Performing image processing for image "; if (node->GetName(name)) { // a property called "name" was found for this DataNode message << "'" << name << "'"; m_OldNodeName = name; } else m_OldNodeName = " "; message << "."; MITK_INFO << message.str(); m_Controls.buttonApplyBModeFilter->setText("working..."); DisableControls(); BmodeThread *thread = new BmodeThread(); connect(thread, &BmodeThread::result, this, &PAImageProcessing::HandleBmodeResults); connect(thread, &BmodeThread::finished, thread, &QObject::deleteLater); - if(m_Controls.BModeMethod->currentText() == "Simple Abs Filter") - thread->setConfig(m_UseLogfilter, m_ResampleSpacing, mitk::PhotoacousticImage::BModeMethod::Abs); - else if(m_Controls.BModeMethod->currentText() == "Shape Detection") - thread->setConfig(m_UseLogfilter, m_ResampleSpacing, mitk::PhotoacousticImage::BModeMethod::ShapeDetection); + bool useGPU = m_Controls.UseGPUBmode->isChecked(); + if(m_Controls.BModeMethod->currentText() == "Absolute Filter") + thread->setConfig(m_UseLogfilter, m_ResampleSpacing, mitk::PhotoacousticImage::BModeMethod::Abs, useGPU); + else if(m_Controls.BModeMethod->currentText() == "Envelope Detection") + thread->setConfig(m_UseLogfilter, m_ResampleSpacing, mitk::PhotoacousticImage::BModeMethod::ShapeDetection, useGPU); thread->setInputImage(image); + thread->setFilterBank(m_FilterBank); MITK_INFO << "Started new thread for Image Processing"; thread->start(); } } } void PAImageProcessing::HandleBmodeResults(mitk::Image::Pointer image) { auto newNode = mitk::DataNode::New(); newNode->SetData(image); // name the new Data node std::stringstream newNodeName; newNodeName << m_OldNodeName << " "; newNodeName << "B-Mode"; newNode->SetName(newNodeName.str()); // update level window for the current dynamic range mitk::LevelWindow levelWindow; newNode->GetLevelWindow(levelWindow); auto data = newNode->GetData(); levelWindow.SetAuto(dynamic_cast(data), true, true); newNode->SetLevelWindow(levelWindow); // add new node to data storage this->GetDataStorage()->Add(newNode); // disable progress bar m_Controls.progressBar->setVisible(false); m_Controls.buttonApplyBModeFilter->setText("Apply B-mode Filter"); EnableControls(); // update rendering mitk::RenderingManager::GetInstance()->InitializeViews( dynamic_cast(data)->GetGeometry(), mitk::RenderingManager::REQUEST_UPDATE_ALL, true); mitk::RenderingManager::GetInstance()->RequestUpdateAll(); } void PAImageProcessing::StartCropThread() { QList nodes = this->GetDataManagerSelection(); if (nodes.empty()) return; mitk::DataStorage::Pointer storage = this->GetDataStorage(); mitk::DataNode::Pointer node = nodes.front(); if (!node) { // Nothing selected. Inform the user and return QMessageBox::information(NULL, "Template", "Please load and select an image before starting image cropping."); return; } mitk::BaseData* data = node->GetData(); if (data) { // test if this data item is an image or not (could also be a surface or something totally different) mitk::Image* image = dynamic_cast(data); if (image) { UpdateBFSettings(image); std::stringstream message; std::string name; message << "Performing image cropping for image "; if (node->GetName(name)) { // a property called "name" was found for this DataNode message << "'" << name << "'"; m_OldNodeName = name; } else m_OldNodeName = " "; message << "."; MITK_INFO << message.str(); m_Controls.buttonApplyCropFilter->setText("working..."); DisableControls(); CropThread *thread = new CropThread(); connect(thread, &CropThread::result, this, &PAImageProcessing::HandleCropResults); connect(thread, &CropThread::finished, thread, &QObject::deleteLater); thread->setConfig(m_Controls.CutoffAbove->value(), m_Controls.CutoffBelow->value()); thread->setInputImage(image); + thread->setFilterBank(m_FilterBank); MITK_INFO << "Started new thread for Image Cropping"; thread->start(); } } } void PAImageProcessing::HandleCropResults(mitk::Image::Pointer image) { auto newNode = mitk::DataNode::New(); newNode->SetData(image); // name the new Data node std::stringstream newNodeName; newNodeName << m_OldNodeName << " "; newNodeName << "Cropped"; newNode->SetName(newNodeName.str()); // update level window for the current dynamic range mitk::LevelWindow levelWindow; newNode->GetLevelWindow(levelWindow); auto data = newNode->GetData(); levelWindow.SetAuto(dynamic_cast(data), true, true); newNode->SetLevelWindow(levelWindow); // add new node to data storage this->GetDataStorage()->Add(newNode); m_Controls.buttonApplyCropFilter->setText("Apply Crop Filter"); EnableControls(); // update rendering mitk::RenderingManager::GetInstance()->InitializeViews( dynamic_cast(data)->GetGeometry(), mitk::RenderingManager::REQUEST_UPDATE_ALL, true); mitk::RenderingManager::GetInstance()->RequestUpdateAll(); } void PAImageProcessing::StartBandpassThread() { QList nodes = this->GetDataManagerSelection(); if (nodes.empty()) return; mitk::DataStorage::Pointer storage = this->GetDataStorage(); mitk::DataNode::Pointer node = nodes.front(); if (!node) { // Nothing selected. Inform the user and return QMessageBox::information(NULL, "Template", "Please load and select an image before starting image cropping."); return; } mitk::BaseData* data = node->GetData(); if (data) { // test if this data item is an image or not (could also be a surface or something totally different) mitk::Image* image = dynamic_cast(data); if (image) { UpdateBFSettings(image); std::stringstream message; std::string name; message << "Performing Bandpass filter on image "; if (node->GetName(name)) { // a property called "name" was found for this DataNode message << "'" << name << "'"; m_OldNodeName = name; } else m_OldNodeName = " "; message << "."; MITK_INFO << message.str(); m_Controls.buttonApplyBandpass->setText("working..."); DisableControls(); BandpassThread *thread = new BandpassThread(); connect(thread, &BandpassThread::result, this, &PAImageProcessing::HandleBandpassResults); connect(thread, &BandpassThread::finished, thread, &QObject::deleteLater); float recordTime = image->GetDimension(1)*image->GetGeometry()->GetSpacing()[1] / 1000 / m_Controls.BPSpeedOfSound->value(); // add a safeguard so the program does not chrash when applying a Bandpass that reaches out of the bounds of the image float maxFrequency = 1 / (recordTime / image->GetDimension(1)) * image->GetDimension(1) / 2 / 2 / 1000; float BPHighPass = 1000000 * m_Controls.BPhigh->value(); // [Hz] float BPLowPass = maxFrequency - 1000000 * m_Controls.BPlow->value(); // [Hz] if (BPLowPass > maxFrequency && m_Controls.UseBP->isChecked()) { QMessageBox Msgbox; Msgbox.setText("LowPass too low, disabled it."); Msgbox.exec(); BPLowPass = 0; } if (BPLowPass < 0 && m_Controls.UseBP->isChecked()) { QMessageBox Msgbox; Msgbox.setText("LowPass too high, disabled it."); Msgbox.exec(); BPLowPass = 0; } if (BPHighPass > maxFrequency && m_Controls.UseBP->isChecked()) { QMessageBox Msgbox; Msgbox.setText("HighPass too high, disabled it."); Msgbox.exec(); BPHighPass = 0; } if (BPHighPass > maxFrequency - BFconfig.BPLowPass) { QMessageBox Msgbox; Msgbox.setText("HighPass higher than LowPass, disabled both."); Msgbox.exec(); BPHighPass = 0; BPLowPass = 0; } thread->setConfig(BPHighPass, BPLowPass, m_Controls.BPFalloff->value(), recordTime); thread->setInputImage(image); + thread->setFilterBank(m_FilterBank); MITK_INFO << "Started new thread for Bandpass filter"; thread->start(); } } } void PAImageProcessing::HandleBandpassResults(mitk::Image::Pointer image) { auto newNode = mitk::DataNode::New(); newNode->SetData(image); // name the new Data node std::stringstream newNodeName; newNodeName << m_OldNodeName << " "; newNodeName << "Bandpassed"; newNode->SetName(newNodeName.str()); // update level window for the current dynamic range mitk::LevelWindow levelWindow; newNode->GetLevelWindow(levelWindow); auto data = newNode->GetData(); levelWindow.SetAuto(dynamic_cast(data), true, true); newNode->SetLevelWindow(levelWindow); // add new node to data storage this->GetDataStorage()->Add(newNode); m_Controls.buttonApplyBandpass->setText("Apply Bandpass"); EnableControls(); // update rendering mitk::RenderingManager::GetInstance()->InitializeViews( dynamic_cast(data)->GetGeometry(), mitk::RenderingManager::REQUEST_UPDATE_ALL, true); mitk::RenderingManager::GetInstance()->RequestUpdateAll(); } void PAImageProcessing::UpdateBounds() { if (!m_Controls.Partial->isChecked()) { m_Controls.boundLow->setEnabled(false); m_Controls.boundHigh->setEnabled(false); BFconfig.partial = false; return; } else { m_Controls.boundLow->setEnabled(true); m_Controls.boundHigh->setEnabled(true); BFconfig.partial = true; } - if(m_Controls.boundLow->value()>m_Controls.boundHigh->value()) + if(m_Controls.boundLow->value() > m_Controls.boundHigh->value()) { - MITK_INFO << "high bound < low bound -> setting both to beamform only first slice"; - m_Controls.boundLow->setValue(0); - m_Controls.boundHigh->setValue(0); - BFconfig.CropBounds[0] = 0; - BFconfig.CropBounds[1] = 0; + m_Controls.boundHigh->setValue(m_Controls.boundLow->value()); + BFconfig.CropBounds[0] = m_Controls.boundLow->value(); + BFconfig.CropBounds[1] = m_Controls.boundHigh->value(); } else { BFconfig.CropBounds[0] = m_Controls.boundLow->value(); BFconfig.CropBounds[1] = m_Controls.boundHigh->value(); } } void PAImageProcessing::UpdateProgress(int progress, std::string progressInfo) { if (progress < 100) m_Controls.progressBar->setValue(progress); else m_Controls.progressBar->setValue(100); m_Controls.ProgressInfo->setText(progressInfo.c_str()); qApp->processEvents(); } void PAImageProcessing::UpdateImageInfo() { QList nodes = this->GetDataManagerSelection(); if (nodes.empty()) return; mitk::DataNode::Pointer node = nodes.front(); if (!node) { // Nothing selected return; } mitk::BaseData* data = node->GetData(); if (data) { // test if this data item is an image or not (could also be a surface or something totally different) mitk::Image* image = dynamic_cast(data); if (image) { // beamforming configs if (m_Controls.UseImageSpacing->isChecked()) { m_Controls.ElementCount->setValue(image->GetDimension(0)); m_Controls.Pitch->setValue(image->GetGeometry()->GetSpacing()[0]); m_Controls.boundLow->setMaximum(image->GetDimension(2) - 1); + m_Controls.boundHigh->setMaximum(image->GetDimension(2) - 1); } UpdateRecordTime(image); // bandpass configs float speedOfSound = m_Controls.BPSpeedOfSound->value(); // [m/s] float recordTime = image->GetDimension(1)*image->GetGeometry()->GetSpacing()[1] / 1000 / speedOfSound; std::stringstream frequency; float maxFrequency = 1 / (recordTime / image->GetDimension(1)) * image->GetDimension(1) / 2 / 2 / 1000; frequency << maxFrequency / 1000000; //[MHz] frequency << "MHz"; m_Controls.BPhigh->setMaximum(maxFrequency / 1000000); m_Controls.BPlow->setMaximum(maxFrequency / 1000000); frequency << " is the maximal allowed frequency for the selected image."; m_Controls.BPhigh->setToolTip(frequency.str().c_str()); m_Controls.BPlow->setToolTip(frequency.str().c_str()); m_Controls.BPhigh->setToolTipDuration(5000); m_Controls.BPlow->setToolTipDuration(5000); } } } void PAImageProcessing::OnSelectionChanged( berry::IWorkbenchPart::Pointer /*source*/, const QList& nodes ) { // iterate all selected objects, adjust warning visibility foreach( mitk::DataNode::Pointer node, nodes ) { if( node.IsNotNull() && dynamic_cast(node->GetData()) ) { m_Controls.labelWarning->setVisible( false ); m_Controls.buttonApplyBModeFilter->setEnabled( true ); m_Controls.labelWarning2->setVisible(false); m_Controls.buttonApplyCropFilter->setEnabled(true); m_Controls.labelWarning3->setVisible(false); m_Controls.buttonApplyBandpass->setEnabled(true); + m_Controls.labelWarning4->setVisible(false); + m_Controls.buttonApplyBeamforming->setEnabled(true); UpdateImageInfo(); return; } } m_Controls.labelWarning->setVisible( true ); m_Controls.buttonApplyBModeFilter->setEnabled( false ); m_Controls.labelWarning2->setVisible(true); m_Controls.buttonApplyCropFilter->setEnabled(false); m_Controls.labelWarning3->setVisible(true); m_Controls.buttonApplyBandpass->setEnabled(false); + m_Controls.labelWarning4->setVisible(true); + m_Controls.buttonApplyBeamforming->setEnabled(false); } void PAImageProcessing::UseResampling() { if (m_Controls.DoResampling->isChecked()) { m_Controls.ResamplingValue->setEnabled(true); m_ResampleSpacing = m_Controls.ResamplingValue->value(); } else { m_Controls.ResamplingValue->setEnabled(false); m_ResampleSpacing = 0; } } void PAImageProcessing::UseLogfilter() { m_UseLogfilter = m_Controls.Logfilter->isChecked(); } void PAImageProcessing::SetResampling() { m_ResampleSpacing = m_Controls.ResamplingValue->value(); } void PAImageProcessing::UpdateBFSettings(mitk::Image::Pointer image) { if ("DAS" == m_Controls.BFAlgorithm->currentText()) - BFconfig.Algorithm = mitk::BeamformingFilter::beamformingSettings::BeamformingAlgorithm::DAS; + BFconfig.Algorithm = mitk::BeamformingSettings::BeamformingAlgorithm::DAS; else if ("DMAS" == m_Controls.BFAlgorithm->currentText()) - BFconfig.Algorithm = mitk::BeamformingFilter::beamformingSettings::BeamformingAlgorithm::DMAS; + BFconfig.Algorithm = mitk::BeamformingSettings::BeamformingAlgorithm::DMAS; if ("Quad. Approx." == m_Controls.DelayCalculation->currentText()) { - BFconfig.DelayCalculationMethod = mitk::BeamformingFilter::beamformingSettings::DelayCalc::QuadApprox; + BFconfig.DelayCalculationMethod = mitk::BeamformingSettings::DelayCalc::QuadApprox; } else if ("Spherical Wave" == m_Controls.DelayCalculation->currentText()) { - BFconfig.DelayCalculationMethod = mitk::BeamformingFilter::beamformingSettings::DelayCalc::Spherical; + BFconfig.DelayCalculationMethod = mitk::BeamformingSettings::DelayCalc::Spherical; } if ("Von Hann" == m_Controls.Apodization->currentText()) { - BFconfig.Apod = mitk::BeamformingFilter::beamformingSettings::Apodization::Hann; + BFconfig.Apod = mitk::BeamformingSettings::Apodization::Hann; } else if ("Hamming" == m_Controls.Apodization->currentText()) { - BFconfig.Apod = mitk::BeamformingFilter::beamformingSettings::Apodization::Hamm; + BFconfig.Apod = mitk::BeamformingSettings::Apodization::Hamm; } else if ("Box" == m_Controls.Apodization->currentText()) { - BFconfig.Apod = mitk::BeamformingFilter::beamformingSettings::Apodization::Box; + BFconfig.Apod = mitk::BeamformingSettings::Apodization::Box; } BFconfig.Pitch = m_Controls.Pitch->value() / 1000; // [m] BFconfig.SpeedOfSound = m_Controls.SpeedOfSound->value(); // [m/s] BFconfig.SamplesPerLine = m_Controls.Samples->value(); BFconfig.ReconstructionLines = m_Controls.Lines->value(); BFconfig.TransducerElements = m_Controls.ElementCount->value(); BFconfig.Angle = m_Controls.Angle->value(); // [deg] BFconfig.UseBP = m_Controls.UseBP->isChecked(); - BFconfig.UseGPU = m_Controls.UseGPU->isChecked(); + BFconfig.UseGPU = m_Controls.UseGPUBf->isChecked(); + BFconfig.upperCutoff = m_Controls.Cutoff->value(); UpdateRecordTime(image); UpdateBounds(); } void PAImageProcessing::UpdateRecordTime(mitk::Image::Pointer image) { if (m_Controls.UseImageSpacing->isChecked()) { BFconfig.RecordTime = image->GetDimension(1)*image->GetGeometry()->GetSpacing()[1] / 1000000; // [s] BFconfig.TimeSpacing = image->GetGeometry()->GetSpacing()[1] / 1000000; - MITK_INFO << "Calculated Scan Depth of " << BFconfig.RecordTime * BFconfig.SpeedOfSound * 100 << "cm"; + MITK_INFO << "Calculated Scan Depth of " << BFconfig.RecordTime * BFconfig.SpeedOfSound * 100 / 2<< "cm"; } else { BFconfig.RecordTime = 2 * m_Controls.ScanDepth->value() / 1000 / BFconfig.SpeedOfSound; // [s] BFconfig.TimeSpacing = BFconfig.RecordTime / image->GetDimension(1); } if ("US Image" == m_Controls.ImageType->currentText()) { - BFconfig.Photoacoustic = false; + BFconfig.isPhotoacousticImage = false; } else if ("PA Image" == m_Controls.ImageType->currentText()) { - BFconfig.Photoacoustic = true; + BFconfig.isPhotoacousticImage = true; } } void PAImageProcessing::EnableControls() { m_Controls.DoResampling->setEnabled(true); UseResampling(); m_Controls.Logfilter->setEnabled(true); m_Controls.buttonApplyBModeFilter->setEnabled(true); m_Controls.CutoffAbove->setEnabled(true); m_Controls.CutoffBelow->setEnabled(true); m_Controls.Cutoff->setEnabled(true); m_Controls.buttonApplyCropFilter->setEnabled(true); m_Controls.buttonApplyBandpass->setEnabled(true); + m_Controls.Partial->setEnabled(true); + m_Controls.boundHigh->setEnabled(true); + m_Controls.boundLow->setEnabled(true); m_Controls.BFAlgorithm->setEnabled(true); m_Controls.DelayCalculation->setEnabled(true); m_Controls.ImageType->setEnabled(true); m_Controls.Apodization->setEnabled(true); m_Controls.UseBP->setEnabled(true); - m_Controls.UseGPU->setEnabled(true); + + #ifdef PHOTOACOUSTICS_USE_GPU + m_Controls.UseGPUBf->setEnabled(true); + m_Controls.UseGPUBmode->setEnabled(true); + #endif + m_Controls.BPhigh->setEnabled(true); m_Controls.BPlow->setEnabled(true); m_Controls.BPFalloff->setEnabled(true); m_Controls.UseImageSpacing->setEnabled(true); UseImageSpacing(); m_Controls.Pitch->setEnabled(true); m_Controls.ElementCount->setEnabled(true); m_Controls.SpeedOfSound->setEnabled(true); m_Controls.Samples->setEnabled(true); m_Controls.Lines->setEnabled(true); m_Controls.Angle->setEnabled(true); m_Controls.buttonApplyBeamforming->setEnabled(true); } void PAImageProcessing::DisableControls() { m_Controls.DoResampling->setEnabled(false); m_Controls.ResamplingValue->setEnabled(false); m_Controls.Logfilter->setEnabled(false); m_Controls.buttonApplyBModeFilter->setEnabled(false); m_Controls.CutoffAbove->setEnabled(false); m_Controls.CutoffBelow->setEnabled(false); m_Controls.Cutoff->setEnabled(false); m_Controls.buttonApplyCropFilter->setEnabled(false); m_Controls.buttonApplyBandpass->setEnabled(false); + m_Controls.Partial->setEnabled(false); + m_Controls.boundHigh->setEnabled(false); + m_Controls.boundLow->setEnabled(false); m_Controls.BFAlgorithm->setEnabled(false); m_Controls.DelayCalculation->setEnabled(false); m_Controls.ImageType->setEnabled(false); m_Controls.Apodization->setEnabled(false); m_Controls.UseBP->setEnabled(false); - m_Controls.UseGPU->setEnabled(false); + + #ifdef PHOTOACOUSTICS_USE_GPU + m_Controls.UseGPUBf->setEnabled(false); + m_Controls.UseGPUBmode->setEnabled(false); + #endif + m_Controls.BPhigh->setEnabled(false); m_Controls.BPlow->setEnabled(false); m_Controls.BPFalloff->setEnabled(false); m_Controls.UseImageSpacing->setEnabled(false); m_Controls.ScanDepth->setEnabled(false); m_Controls.Pitch->setEnabled(false); m_Controls.ElementCount->setEnabled(false); m_Controls.SpeedOfSound->setEnabled(false); m_Controls.Samples->setEnabled(false); m_Controls.Lines->setEnabled(false); m_Controls.Angle->setEnabled(false); m_Controls.buttonApplyBeamforming->setEnabled(false); } void PAImageProcessing::UseImageSpacing() { if (m_Controls.UseImageSpacing->isChecked()) { m_Controls.ScanDepth->setDisabled(true); } else { m_Controls.ScanDepth->setEnabled(true); } } void BeamformingThread::run() { mitk::Image::Pointer resultImage; - mitk::PhotoacousticImage::Pointer filterbank = mitk::PhotoacousticImage::New(); std::function progressHandle = [this](int progress, std::string progressInfo) { emit updateProgress(progress, progressInfo); }; - resultImage = filterbank->ApplyBeamforming(m_InputImage, m_BFconfig, m_Cutoff, progressHandle); + resultImage = m_FilterBank->ApplyBeamforming(m_InputImage, m_BFconfig, progressHandle); emit result(resultImage); } -void BeamformingThread::setConfig(mitk::BeamformingFilter::beamformingSettings BFconfig) +void BeamformingThread::setConfig(mitk::BeamformingSettings BFconfig) { m_BFconfig = BFconfig; } void BeamformingThread::setInputImage(mitk::Image::Pointer image) { m_InputImage = image; } -void BeamformingThread::setCutoff(int cutoff) -{ - m_Cutoff = cutoff; -} - void BmodeThread::run() { mitk::Image::Pointer resultImage; - mitk::PhotoacousticImage::Pointer filterbank = mitk::PhotoacousticImage::New(); - resultImage = filterbank->ApplyBmodeFilter(m_InputImage, m_Method, m_UseLogfilter, m_ResampleSpacing); + resultImage = m_FilterBank->ApplyBmodeFilter(m_InputImage, m_Method, m_UseGPU, m_UseLogfilter, m_ResampleSpacing); emit result(resultImage); } -void BmodeThread::setConfig(bool useLogfilter, double resampleSpacing, mitk::PhotoacousticImage::BModeMethod method) +void BmodeThread::setConfig(bool useLogfilter, double resampleSpacing, mitk::PhotoacousticImage::BModeMethod method, bool useGPU) { m_UseLogfilter = useLogfilter; m_ResampleSpacing = resampleSpacing; m_Method = method; + m_UseGPU = useGPU; } void BmodeThread::setInputImage(mitk::Image::Pointer image) { m_InputImage = image; } void CropThread::run() { mitk::Image::Pointer resultImage; - mitk::PhotoacousticImage::Pointer filterbank = mitk::PhotoacousticImage::New(); - resultImage = filterbank->ApplyCropping(m_InputImage, m_CutAbove, m_CutBelow, 0, 0, 0, m_InputImage->GetDimension(2) - 1); + resultImage = m_FilterBank->ApplyCropping(m_InputImage, m_CutAbove, m_CutBelow, 0, 0, 0, m_InputImage->GetDimension(2) - 1); emit result(resultImage); } void CropThread::setConfig(unsigned int CutAbove, unsigned int CutBelow) { m_CutAbove = CutAbove; m_CutBelow = CutBelow; } void CropThread::setInputImage(mitk::Image::Pointer image) { m_InputImage = image; } void BandpassThread::run() { mitk::Image::Pointer resultImage; - mitk::PhotoacousticImage::Pointer filterbank = mitk::PhotoacousticImage::New(); - resultImage = filterbank->BandpassFilter(m_InputImage, m_RecordTime, m_BPHighPass, m_BPLowPass, m_TukeyAlpha); + resultImage = m_FilterBank->BandpassFilter(m_InputImage, m_RecordTime, m_BPHighPass, m_BPLowPass, m_TukeyAlpha); emit result(resultImage); } void BandpassThread::setConfig(float BPHighPass, float BPLowPass, float TukeyAlpha, float recordTime) { m_BPHighPass = BPHighPass; m_BPLowPass = BPLowPass; m_TukeyAlpha = TukeyAlpha; m_RecordTime = recordTime; } void BandpassThread::setInputImage(mitk::Image::Pointer image) { m_InputImage = image; } diff --git a/Plugins/org.mitk.gui.qt.photoacoustics.imageprocessing/src/internal/PAImageProcessing.h b/Plugins/org.mitk.gui.qt.photoacoustics.imageprocessing/src/internal/PAImageProcessing.h index 42270dd4c7..59b649bde8 100644 --- a/Plugins/org.mitk.gui.qt.photoacoustics.imageprocessing/src/internal/PAImageProcessing.h +++ b/Plugins/org.mitk.gui.qt.photoacoustics.imageprocessing/src/internal/PAImageProcessing.h @@ -1,182 +1,214 @@ /*=================================================================== 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 PAImageProcessing_h #define PAImageProcessing_h #include + #include #include #include #include "ui_PAImageProcessingControls.h" #include "mitkPhotoacousticBeamformingFilter.h" +#include "mitkPhotoacousticBeamformingSettings.h" Q_DECLARE_METATYPE(mitk::Image::Pointer) Q_DECLARE_METATYPE(std::string) class PAImageProcessing : public QmitkAbstractView { // this is needed for all Qt objects that should have a Qt meta-object // (everything that derives from QObject and wants to have signal/slots) Q_OBJECT public: static const std::string VIEW_ID; PAImageProcessing(); protected slots: /// \brief Called when the user clicks the GUI button void UpdateBounds(); void UseResampling(); void UseLogfilter(); void SetResampling(); void UseImageSpacing(); void UpdateImageInfo(); void UpdateRecordTime(mitk::Image::Pointer image); void HandleBeamformingResults(mitk::Image::Pointer image); void StartBeamformingThread(); void HandleBmodeResults(mitk::Image::Pointer image); void StartBmodeThread(); void HandleCropResults(mitk::Image::Pointer image); void StartCropThread(); void HandleBandpassResults(mitk::Image::Pointer image); void StartBandpassThread(); void UpdateProgress(int progress, std::string progressInfo); void BatchProcessing(); void UpdateSaveBoxes(); + void ChangedSOSBandpass(); + void ChangedSOSBeamforming(); + protected: virtual void CreateQtPartControl(QWidget *parent) override; virtual void SetFocus() override; /// \brief called by QmitkFunctionality when DataManager's selection has changed virtual void OnSelectionChanged( berry::IWorkbenchPart::Pointer source, const QList& nodes ) override; Ui::PAImageProcessingControls m_Controls; float m_ResampleSpacing; bool m_UseLogfilter; std::string m_OldNodeName; - mitk::BeamformingFilter::beamformingSettings BFconfig; + mitk::BeamformingSettings BFconfig; void UpdateBFSettings(mitk::Image::Pointer image); void EnableControls(); void DisableControls(); + + mitk::PhotoacousticImage::Pointer m_FilterBank; }; class BeamformingThread : public QThread { Q_OBJECT void run() Q_DECL_OVERRIDE; signals: void result(mitk::Image::Pointer); void updateProgress(int, std::string); public: - void setConfig(mitk::BeamformingFilter::beamformingSettings BFconfig); + void setConfig(mitk::BeamformingSettings BFconfig); void setInputImage(mitk::Image::Pointer image); - void setCutoff(int cutoff); + void setFilterBank(mitk::PhotoacousticImage::Pointer filterBank) + { + m_FilterBank = filterBank; + } + protected: - mitk::BeamformingFilter::beamformingSettings m_BFconfig; + mitk::BeamformingSettings m_BFconfig; mitk::Image::Pointer m_InputImage; int m_Cutoff; + + mitk::PhotoacousticImage::Pointer m_FilterBank; }; class BmodeThread : public QThread { Q_OBJECT void run() Q_DECL_OVERRIDE; signals: void result(mitk::Image::Pointer); public: enum BModeMethod { ShapeDetection, Abs }; - void setConfig(bool useLogfilter, double resampleSpacing, mitk::PhotoacousticImage::BModeMethod method); + void setConfig(bool useLogfilter, double resampleSpacing, mitk::PhotoacousticImage::BModeMethod method, bool useGPU); void setInputImage(mitk::Image::Pointer image); + void setFilterBank(mitk::PhotoacousticImage::Pointer filterBank) + { + m_FilterBank = filterBank; + } + protected: mitk::Image::Pointer m_InputImage; mitk::PhotoacousticImage::BModeMethod m_Method; bool m_UseLogfilter; double m_ResampleSpacing; + bool m_UseGPU; + + mitk::PhotoacousticImage::Pointer m_FilterBank; }; class CropThread : public QThread { Q_OBJECT void run() Q_DECL_OVERRIDE; signals: void result(mitk::Image::Pointer); public: void setConfig(unsigned int CutAbove, unsigned int CutBelow); void setInputImage(mitk::Image::Pointer image); + void setFilterBank(mitk::PhotoacousticImage::Pointer filterBank) + { + m_FilterBank = filterBank; + } protected: mitk::Image::Pointer m_InputImage; unsigned int m_CutAbove; unsigned int m_CutBelow; + + mitk::PhotoacousticImage::Pointer m_FilterBank; }; class BandpassThread : public QThread { Q_OBJECT void run() Q_DECL_OVERRIDE; signals: void result(mitk::Image::Pointer); public: void setConfig(float BPHighPass, float BPLowPass, float TukeyAlpha, float recordTime); void setInputImage(mitk::Image::Pointer image); + void setFilterBank(mitk::PhotoacousticImage::Pointer filterBank) + { + m_FilterBank = filterBank; + } protected: mitk::Image::Pointer m_InputImage; float m_BPHighPass; float m_BPLowPass; float m_TukeyAlpha; float m_RecordTime; + + mitk::PhotoacousticImage::Pointer m_FilterBank; }; #endif // PAImageProcessing_h diff --git a/Plugins/org.mitk.gui.qt.photoacoustics.imageprocessing/src/internal/PAImageProcessingControls.ui b/Plugins/org.mitk.gui.qt.photoacoustics.imageprocessing/src/internal/PAImageProcessingControls.ui index bbc2c814b3..dfd5d7636c 100644 --- a/Plugins/org.mitk.gui.qt.photoacoustics.imageprocessing/src/internal/PAImageProcessingControls.ui +++ b/Plugins/org.mitk.gui.qt.photoacoustics.imageprocessing/src/internal/PAImageProcessingControls.ui @@ -1,960 +1,987 @@ PAImageProcessingControls 0 0 382 1217 0 0 QmitkTemplate <html><head/><body><p><span style=" font-weight:600;">Batch Processing</span></p></body></html> Start Batch Processing Bandpass true Crop true Save true Save true Save Beamform true BMode true Save true <html><head/><body><p><span style=" font-weight:600;">B-mode Filter Settings</span></p></body></html> Do Resampling Absolute Filter Envelope Detection 0 0 13 0 11 0.010000000000000 1.000000000000000 0.010000000000000 0.100000000000000 [mm] Depth Spacing Add Logfilter + + + + Use GPU + + + QLabel { color: rgb(255, 0, 0) } <html><head/><body><p align="center"><span style=" font-size:10pt; font-weight:600;">Please select an image!</span></p></body></html> 0 0 Do image processing Apply B-mode Filter Qt::Horizontal <html><head/><body><p><span style=" font-weight:600;">Bandpass Filter Settings</span></p></body></html> QLayout::SetDefaultConstraint 0 0 0 3 0.010000000000000 200.000000000000000 15.000000000000000 [MHz] f High Pass [MHz] f Low Pass 0 0 3 200.000000000000000 Thukey window alpha 1 200.000000000000000 3000.000000000000000 5.000000000000000 1480.000000000000000 [m/s] Speed of Sound 2 1.000000000000000 0.100000000000000 0.500000000000000 <html><head/><body><p align="center"><span style=" font-size:10pt; font-weight:600; color:#ff0000;">Please select an image!</span></p></body></html> Apply Bandpass Qt::Horizontal <html><head/><body><p><span style=" font-weight:600;">Crop Filter Settings</span></p></body></html> 99999 5 165 Cut Top 99999 10 Cut Bottom <html><head/><body><p align="center"><span style=" font-size:10pt; font-weight:600; color:#ff0000;">Please select an image!</span></p></body></html> Apply Crop Filer Qt::Horizontal <html><head/><body><p><span style=" font-weight:600;">Beamforming Filter Settings</span></p></body></html> 5 2 Delay Calculation Auto Get Depth true Apply Beamforming Beamforming Method [mm] Scan Depth 0 0 - 2 + 3 0.010000000000000 9.000000000000000 0.050000000000000 0.300000000000000 Transducer Elements 0 0 4 300.000000000000000 0.100000000000000 50.000000000000000 [mm] Transducer Pitch 0 0 64 1024 128 128 0 0 256 16384 256 2048 0 0 64 2048 128 256 Samples Reconstruction Lines true 0 0 100 0 0 0 900 10 0 - - - - Cutoff Upper Voxels - - - 0 0 DAS DMAS 0 0 Quad. Approx. Spherical Wave 0 0 PA Image US Image Image Type 0 0 Von Hann Hamming Box Apodization 0 0 1 200.000000000000000 3000.000000000000000 5.000000000000000 1480.000000000000000 [m/s] Speed of Sound false 99999 minimal beamformed slice min false 99999 10 Maximal beamformed slice max select slices - + Compute On GPU true true Auto Use Bandpass 0 0 1 1.000000000000000 180.000000000000000 27.000000000000000 [°] Element Angle + + + + Cutoff Upper Voxels + + + + + + + <html><head/><body><p align="center"><span style=" font-size:10pt; font-weight:600; color:#ff0000;">Please select an image!</span></p></body></html> + + + + + + + Qt::Vertical + + + + 20 + 40 + + + +