diff --git a/Modules/PhotoacousticsAlgorithms/Resources/DAS.cl b/Modules/PhotoacousticsAlgorithms/Resources/DAS.cl index 34c70224f0..42d954137b 100644 --- a/Modules/PhotoacousticsAlgorithms/Resources/DAS.cl +++ b/Modules/PhotoacousticsAlgorithms/Resources/DAS.cl @@ -1,124 +1,116 @@ /*=================================================================== 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 short* delays, __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 ) { float l_i = (float)globalPosX / (float)outputL * (float)inputL; - 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]; + unsigned short curUsedLines = inputL; float apod_mult = (float)apodArraySize / (float)curUsedLines; unsigned short Delay = 0; float output = 0; float mult = 0; - for (short l_s = minLine; l_s < maxLine; ++l_s) + for (short l_s = 0; l_s < inputL; ++l_s) { Delay = delays[globalPosY * (outputL / 2) + (int)(fabs(l_s - l_i)/(float)inputL * (float)outputL)]; if (Delay < inputS && Delay >= 0) { - output += apodArray[(int)((l_s - minLine)*apod_mult)] * dSource[(int)(globalPosZ * inputL * inputS + Delay * inputL + l_s)]; + output += apodArray[(int)((l_s)*apod_mult)] * dSource[(int)(globalPosZ * inputL * inputS + Delay * inputL + l_s)]; } else --curUsedLines; } dDest[ globalPosZ * outputL * outputS + globalPosY * outputL + globalPosX ] = output / (float)curUsedLines; } } __kernel void ckDAS_g( __global float* dSource, // input image __global float* dDest, // output buffer __global float* elementHeights, __global float* elementPositions, __constant float* apodArray, unsigned short apodArraySize, unsigned int inputL, unsigned int inputS, int Slices, int outputL, int outputS, float totalSamples_i, float horizontalExtent, float mult, - char isPAImage, - __global unsigned short* usedLines // parameters + char isPAImage ) { // get thread identifier int globalPosX = get_global_id(0); int globalPosY = get_global_id(1); int globalPosZ = get_global_id(2); // terminate non-valid threads if ( globalPosX < outputL && globalPosY < outputS && globalPosZ < Slices ) { int AddSample = 0; float l_p = 0; float s_i = 0; - float apod_mult = 1; - float output = 0; l_p = (float)globalPosX / outputL * horizontalExtent; s_i = (float)globalPosY / outputS * totalSamples_i; - 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]; + unsigned short curUsedLines = inputL; apod_mult = (float)apodArraySize / curUsedLines; - for (int l_s = minLine; l_s < maxLine; ++l_s) + for (int l_s = 0; l_s < inputL; ++l_s) { AddSample = (int)sqrt( pow(s_i-elementHeights[l_s]*mult, 2) + pow(mult * (l_p - elementPositions[l_s]), 2) ) + (1 - isPAImage)*s_i; if (AddSample < inputS && AddSample >= 0) - output += dSource[(int)(globalPosZ * inputL * inputS + l_s + AddSample*inputL)] * apodArray[(int)((l_s - minLine)*apod_mult)]; + output += dSource[(int)(globalPosZ * inputL * inputS + l_s + AddSample*inputL)] * apodArray[(int)((l_s)*apod_mult)]; else --curUsedLines; } dDest[ globalPosZ * outputL * outputS + globalPosY * outputL + globalPosX ] = output / curUsedLines; } } diff --git a/Modules/PhotoacousticsAlgorithms/Resources/DelayCalculation.cl b/Modules/PhotoacousticsAlgorithms/Resources/DelayCalculation.cl index b71208e6e8..cb7ac05f4e 100644 --- a/Modules/PhotoacousticsAlgorithms/Resources/DelayCalculation.cl +++ b/Modules/PhotoacousticsAlgorithms/Resources/DelayCalculation.cl @@ -1,44 +1,43 @@ /*=================================================================== 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 ckDelayCalculationSphe( __global unsigned short *gDest, - __global unsigned short *usedLines, unsigned int inputL, unsigned int inputS, unsigned int outputL, unsigned int outputS, char isPAImage, float delayMultiplicatorRaw, float totalSamples_i // parameters ) { uint globalPosX = get_global_id(0); uint globalPosY = get_global_id(1); if (globalPosX * 2 < outputL && globalPosY < outputS) { float l_i = 0; // we calculate the delays relative to line zero float s_i = (float)globalPosY / (float)outputS * totalSamples_i; float l_s = (float)globalPosX / (float)outputL * (float)inputL; // the currently calculated line gDest[globalPosY * (outputL / 2) + globalPosX] = 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/UsedLinesCalculation.cl b/Modules/PhotoacousticsAlgorithms/Resources/UsedLinesCalculation.cl deleted file mode 100644 index 132c7d88fd..0000000000 --- a/Modules/PhotoacousticsAlgorithms/Resources/UsedLinesCalculation.cl +++ /dev/null @@ -1,170 +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. - -===================================================================*/ - -__kernel void ckUsedLines( - __global unsigned short* dDest, // output buffer - float partMult, - unsigned int inputL, - unsigned int inputS, - unsigned int outputL, - unsigned int outputS, - float totalSamples_i // parameters -) -{ - // get thread identifier - unsigned int globalPosX = get_global_id(0); - unsigned int globalPosY = get_global_id(1); - - // terminate non-valid threads - if ( globalPosX < outputL && globalPosY < outputS) - { - float l_i = (float)globalPosX / outputL * inputL; - float s_i = (float)globalPosY / (float)outputS * totalSamples_i; - - 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 - } -} - -__kernel void ckUsedLines_g( - __global unsigned short* dDest, // output buffer - __global float* elementHeights, - __global float* elementPositions, - float cos_deg, - float probeRadius, - unsigned int inputL, - unsigned int inputS, - unsigned int outputL, - unsigned int outputS, - float horizontalExtent, - float verticalExtent -) -{ - // get thread identifier - unsigned int globalPosX = get_global_id(0); - unsigned int globalPosY = get_global_id(1); - - // terminate non-valid threads - if ( globalPosX < outputL && globalPosY < outputS) - { - float cos = 0; - float a = 0; - float d = 0; - - float l_p = (float)globalPosX / outputL * horizontalExtent; - float s_p = (float)globalPosY / (float)outputS * verticalExtent; - - int maxLine = inputL; - int minLine = 0; - - for(int l_s = 0; l_s < inputL; l_s+=32) - { - a = sqrt((probeRadius - s_p)*(probeRadius - s_p) + (l_p - horizontalExtent/2)*(l_p - horizontalExtent/2)); - d = sqrt((s_p - elementHeights[l_s])*(s_p - elementHeights[l_s]) + (l_p - elementPositions[l_s])*(l_p - elementPositions[l_s])); - cos = (d*d + probeRadius*probeRadius - a*a)/(2*probeRadius*d); - - if(cos > cos_deg) - { - minLine = l_s-32; - if(minLine < 0) - minLine = 0; - break; - } - } - for(int l_s = minLine; l_s < inputL; l_s+=8) - { - a = sqrt((probeRadius - s_p)*(probeRadius - s_p) + (l_p - horizontalExtent/2)*(l_p - horizontalExtent/2)); - d = sqrt((s_p - elementHeights[l_s])*(s_p - elementHeights[l_s]) + (l_p - elementPositions[l_s])*(l_p - elementPositions[l_s])); - cos = (d*d + probeRadius*probeRadius - a*a)/(2*probeRadius*d); - - if(cos > cos_deg) - { - minLine = l_s-8; - if(minLine < 0) - minLine = 0; - break; - } - } - for(int l_s = minLine; l_s < inputL; l_s+=1) - { - a = sqrt((probeRadius - s_p)*(probeRadius - s_p) + (l_p - horizontalExtent/2)*(l_p - horizontalExtent/2)); - d = sqrt((s_p - elementHeights[l_s])*(s_p - elementHeights[l_s]) + (l_p - elementPositions[l_s])*(l_p - elementPositions[l_s])); - cos = (d*d + probeRadius*probeRadius - a*a)/(2*probeRadius*d); - - if(cos > cos_deg) - { - minLine = l_s; - break; - } - } - - for(int l_s = inputL; l_s >= 0 ; l_s-=32) - { - a = sqrt((probeRadius - s_p)*(probeRadius - s_p) + (l_p - horizontalExtent/2)*(l_p - horizontalExtent/2)); - d = sqrt((s_p - elementHeights[l_s])*(s_p - elementHeights[l_s]) + (l_p - elementPositions[l_s])*(l_p - elementPositions[l_s])); - cos = (d*d + probeRadius*probeRadius - a*a)/(2*probeRadius*d); - cos = 0; - - if(cos > cos_deg) - { - maxLine = l_s+32; - if(maxLine > inputL) - minLine = inputL; - break; - } - } - for(int l_s = maxLine; l_s >= 0 ; l_s-=8) - { - a = sqrt((probeRadius - s_p)*(probeRadius - s_p) + (l_p - horizontalExtent/2)*(l_p - horizontalExtent/2)); - d = sqrt((s_p - elementHeights[l_s])*(s_p - elementHeights[l_s]) + (l_p - elementPositions[l_s])*(l_p - elementPositions[l_s])); - cos = (d*d + probeRadius*probeRadius - a*a)/(2*probeRadius*d); - cos = 0; - - if(cos > cos_deg) - { - maxLine = l_s+8; - if(maxLine > inputL) - minLine = inputL; - break; - } - } - for(int l_s = maxLine; l_s >= 0 ; l_s-=1) - { - a = sqrt((probeRadius - s_p)*(probeRadius - s_p) + (l_p - horizontalExtent/2)*(l_p - horizontalExtent/2)); - d = sqrt((s_p - elementHeights[l_s])*(s_p - elementHeights[l_s]) + (l_p - elementPositions[l_s])*(l_p - elementPositions[l_s])); - cos = (d*d + probeRadius*probeRadius - a*a)/(2*probeRadius*d); - cos = 0; - - if(cos > cos_deg) - { - maxLine = l_s; - break; - } - } - - 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 87f8116092..7b900861de 100644 --- a/Modules/PhotoacousticsAlgorithms/files.cmake +++ b/Modules/PhotoacousticsAlgorithms/files.cmake @@ -1,31 +1,29 @@ file(GLOB_RECURSE H_FILES RELATIVE "${CMAKE_CURRENT_SOURCE_DIR}" "${CMAKE_CURRENT_SOURCE_DIR}/include/*") set(CPP_FILES source/filters/mitkBeamformingFilter.cpp source/filters/mitkBeamformingSettings.cpp source/filters/mitkImageSliceSelectionFilter.cpp source/filters/mitkCastToFloatImageFilter.cpp source/filters/mitkCropImageFilter.cpp source/filters/mitkBandpassFilter.cpp source/OpenCLFilter/mitkPhotoacousticBModeFilter.cpp source/utils/mitkPhotoacousticFilterService.cpp source/utils/mitkBeamformingUtils.cpp source/mitkPhotoacousticMotionCorrectionFilter.cpp ) IF(MITK_USE_OpenCL) list(APPEND CPP_FILES source/OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.cpp - source/OpenCLFilter/mitkPhotoacousticOCLUsedLinesCalculation.cpp source/OpenCLFilter/mitkPhotoacousticOCLDelayCalculation.cpp ) ENDIF(MITK_USE_OpenCL) set(RESOURCE_FILES - UsedLinesCalculation.cl DelayCalculation.cl DMAS.cl DAS.cl sDMAS.cl ) diff --git a/Modules/PhotoacousticsAlgorithms/include/OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.h b/Modules/PhotoacousticsAlgorithms/include/OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.h index bd9bf91913..00799f0ac6 100644 --- a/Modules/PhotoacousticsAlgorithms/include/OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.h +++ b/Modules/PhotoacousticsAlgorithms/include/OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.h @@ -1,147 +1,145 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef _MITKPHOTOACOUSTICSOCLBEAMFORMER_H_ #define _MITKPHOTOACOUSTICSOCLBEAMFORMER_H_ #include #if defined(PHOTOACOUSTICS_USE_GPU) || DOXYGEN #include "mitkOclDataSetToDataSetFilter.h" #include "mitkBeamformingSettings.h" #include "mitkPhotoacousticOCLDelayCalculation.h" -#include "mitkPhotoacousticOCLUsedLinesCalculation.h" #include namespace mitk { /*! * \brief Class implementing a mitk::OclDataSetToDataSetFilter for beamforming on GPU * * The class must be given a configuration class instance of mitk::BeamformingSettings for beamforming parameters through mitk::PhotoacousticOCLBeamformingFilter::SetConfig(BeamformingSettings settings) * Additional configuration of the apodisation function is needed. */ class PhotoacousticOCLBeamformingFilter : public OclDataSetToDataSetFilter, public itk::Object { public: mitkClassMacroItkParent(PhotoacousticOCLBeamformingFilter, itk::Object); mitkNewMacro1Param(Self, BeamformingSettings::Pointer); /** * @brief SetInput Set the input data through an image. Arbitrary images are supported */ void SetInput(Image::Pointer image); /** * brief SetInput Manually set the input data while providing 3 dimensions and memory size of the input data (Bytes per element). */ void SetInput(void* data, unsigned int* dimensions, unsigned int BpE); /** * @brief GetOutput Get a pointer to the processed data. The standard datatype is float. */ void* GetOutput(); /** * @brief GetOutputAsImage Returns an mitk::Image constructed from the processed data */ mitk::Image::Pointer GetOutputAsImage(); /** \brief Update the filter */ void Update(); /** \brief Set the Apodisation function to apply when beamforming */ void SetApodisation(const float* apodisation, unsigned short apodArraySize) { m_ApodArraySize = apodArraySize; m_Apodisation = apodisation; } protected: PhotoacousticOCLBeamformingFilter(BeamformingSettings::Pointer settings); virtual ~PhotoacousticOCLBeamformingFilter(); /** \brief Initialize the filter */ bool Initialize(); /** \brief Updated the used data for beamforming depending on whether the configuration has significantly changed */ void UpdateDataBuffers(); /** \brief Execute the filter */ void Execute(); mitk::PixelType GetOutputType() { return mitk::MakeScalarPixelType(); } int GetBytesPerElem() { return sizeof(float); } virtual us::Module* GetModule(); private: /** The OpenCL kernel for the filter */ cl_kernel m_PixelCalculation; unsigned int m_OutputDim[3]; const float* m_Apodisation; unsigned short m_ApodArraySize; unsigned int m_inputSlices; unsigned short m_PAImage; BeamformingSettings::Pointer m_Conf; mitk::Image::Pointer m_InputImage; size_t m_ChunkSize[3]; - mitk::OCLUsedLinesCalculation::Pointer m_UsedLinesCalculation; mitk::OCLDelayCalculation::Pointer m_DelayCalculation; cl_mem m_ApodizationBuffer; cl_mem m_DelaysBuffer; - cl_mem m_UsedLinesBuffer; + cl_mem m_ElementHeightsBuffer; cl_mem m_ElementPositionsBuffer; }; } #else namespace mitk { class PhotoacousticOCLBeamformingFilter : public itk::Object { public: mitkClassMacroItkParent(mitk::PhotoacousticOCLBeamformingFilter, itk::Object); itkNewMacro(Self); protected: /** Constructor */ PhotoacousticOCLBeamformingFilter() {} /** Destructor */ ~PhotoacousticOCLBeamformingFilter() override {} }; } #endif #endif diff --git a/Modules/PhotoacousticsAlgorithms/include/OpenCLFilter/mitkPhotoacousticOCLDelayCalculation.h b/Modules/PhotoacousticsAlgorithms/include/OpenCLFilter/mitkPhotoacousticOCLDelayCalculation.h index 3e60b4ec5c..69960e3356 100644 --- a/Modules/PhotoacousticsAlgorithms/include/OpenCLFilter/mitkPhotoacousticOCLDelayCalculation.h +++ b/Modules/PhotoacousticsAlgorithms/include/OpenCLFilter/mitkPhotoacousticOCLDelayCalculation.h @@ -1,89 +1,78 @@ /*=================================================================== 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 _MITKPHOTOACOUSTICSDELAYCALC_H_ #define _MITKPHOTOACOUSTICSDELAYCALC_H_ #if defined(PHOTOACOUSTICS_USE_GPU) || DOXYGEN #include "mitkOclDataSetToDataSetFilter.h" #include #include "mitkBeamformingSettings.h" namespace mitk { /*! * \brief Class implementing a mitk::OclDataSetToDataSetFilter to calculate the delays used for beamforming. * * The class must be given a configuration class instance of mitk::BeamformingSettings for beamforming parameters through mitk::OCLDelayCalculation::SetConfig(BeamformingSettings conf) - * Additionally the output of an instance of mitk::OCLUsedLinesCalculation is needed to calculate the delays. */ class OCLDelayCalculation : public OclDataSetToDataSetFilter, public itk::Object { public: mitkClassMacroItkParent(OCLDelayCalculation, itk::Object); mitkNewMacro1Param(Self, mitk::BeamformingSettings::Pointer); void Update(); - /** \brief Sets the usedLines buffer object to use for the calculation of the delays. - * - * @param usedLines An buffer generated as the output of an instance of mitk::OCLUsedLinesCalculation. - */ - void SetInputs(cl_mem usedLines) - { - m_UsedLines = usedLines; - } - protected: OCLDelayCalculation(mitk::BeamformingSettings::Pointer settings); virtual ~OCLDelayCalculation(); /** Initialize the filter */ bool Initialize(); void Execute(); mitk::PixelType GetOutputType() { return mitk::MakeScalarPixelType(); } int GetBytesPerElem() { return sizeof(unsigned short); } virtual us::Module* GetModule(); int m_sizeThis; private: /** The OpenCL kernel for the filter */ cl_kernel m_PixelCalculation; BeamformingSettings::Pointer m_Conf; - cl_mem m_UsedLines; unsigned int m_BufferSize; float m_DelayMultiplicatorRaw; char m_IsPAImage; size_t m_ChunkSize[3]; }; } #endif #endif diff --git a/Modules/PhotoacousticsAlgorithms/include/OpenCLFilter/mitkPhotoacousticOCLUsedLinesCalculation.h b/Modules/PhotoacousticsAlgorithms/include/OpenCLFilter/mitkPhotoacousticOCLUsedLinesCalculation.h deleted file mode 100644 index 52021ea7c5..0000000000 --- a/Modules/PhotoacousticsAlgorithms/include/OpenCLFilter/mitkPhotoacousticOCLUsedLinesCalculation.h +++ /dev/null @@ -1,84 +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. - -===================================================================*/ - -#ifndef _MITKPHOTOACOUSTICSOCLUSEDLINESCALCULATION_H_ -#define _MITKPHOTOACOUSTICSOCLUSEDLINESCALCULATION_H_ - -#if defined(PHOTOACOUSTICS_USE_GPU) || DOXYGEN - -#include "mitkOclDataSetToDataSetFilter.h" -#include -#include "mitkBeamformingSettings.h" - -namespace mitk -{ - /*! - * \brief Class implementing a mitk::OclDataSetToDataSetFilter to calculate which lines each sample should use when beamforming. - * - * The class must be given a configuration class instance of mitk::BeamformingSettings for beamforming parameters through mitk::OCLDelayCalculation::SetConfig(BeamformingSettings conf) - */ - - class OCLUsedLinesCalculation : public OclDataSetToDataSetFilter, public itk::Object - { - public: - mitkClassMacroItkParent(OCLUsedLinesCalculation, itk::Object); - mitkNewMacro1Param(Self, mitk::BeamformingSettings::Pointer); - - void Update(); - - void SetElementHeightsBuffer(cl_mem elementHeightsBuffer); - void SetElementPositionsBuffer(cl_mem elementPositionsBuffer); - - protected: - - /** Constructor */ - OCLUsedLinesCalculation(mitk::BeamformingSettings::Pointer settings); - - /** Destructor */ - virtual ~OCLUsedLinesCalculation(); - - /** Initialize the filter */ - bool Initialize(); - - void Execute(); - - mitk::PixelType GetOutputType() - { - return mitk::MakeScalarPixelType(); - } - - int GetBytesPerElem() - { - return sizeof(unsigned short); - } - - virtual us::Module* GetModule(); - - int m_sizeThis; - - private: - /** The OpenCL kernel for the filter */ - cl_kernel m_PixelCalculation; - - BeamformingSettings::Pointer m_Conf; - float m_part; - size_t m_ChunkSize[3]; - cl_mem m_ElementHeightsBuffer; - cl_mem m_ElementPositionsBuffer; - }; -} -#endif -#endif diff --git a/Modules/PhotoacousticsAlgorithms/include/mitkBeamformingFilter.h b/Modules/PhotoacousticsAlgorithms/include/mitkBeamformingFilter.h index 6ee498e470..8724808d36 100644 --- a/Modules/PhotoacousticsAlgorithms/include/mitkBeamformingFilter.h +++ b/Modules/PhotoacousticsAlgorithms/include/mitkBeamformingFilter.h @@ -1,90 +1,89 @@ /*=================================================================== 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 "mitkBeamformingSettings.h" #include "mitkBeamformingUtils.h" #include "MitkPhotoacousticsAlgorithmsExports.h" namespace mitk { /*! * \brief Class implementing an mitk::ImageToImageFilter for beamforming on both CPU and GPU * * The class must be given a configuration class instance of mitk::BeamformingSettings for beamforming parameters through mitk::BeamformingFilter::Configure(BeamformingSettings settings) * Whether the GPU is used can be set in the configuration. * For significant problems or important messages a string is written, which can be accessed via GetMessageString(). */ class MITKPHOTOACOUSTICSALGORITHMS_EXPORT BeamformingFilter : public ImageToImageFilter { public: mitkClassMacro(BeamformingFilter, ImageToImageFilter); mitkNewMacro1Param(Self, mitk::BeamformingSettings::Pointer); itkCloneMacro(Self); /** \brief Sets a callback for progress checking * * An std::function can be set, through which progress of the currently updating filter is reported. * The integer argument is a number between 0 an 100 to indicate how far completion has been achieved, the std::string argument indicates what the filter is currently doing. */ void SetProgressHandle(std::function progressHandle); protected: BeamformingFilter(mitk::BeamformingSettings::Pointer settings); ~BeamformingFilter() override; void GenerateInputRequestedRegion() override; void GenerateOutputInformation() override; void GenerateData() override; //##Description //## @brief Time when Header was last initialized itk::TimeStamp m_TimeOfHeaderInitialization; /** \brief The std::function, through which progress of the currently updating filter is reported. */ std::function m_ProgressHandle; float* m_OutputData; float* m_InputData; float* m_InputDataPuffer; - unsigned short* m_MinMaxLines; /** \brief Current configuration set */ BeamformingSettings::Pointer m_Conf; /** * The size of the apodization array when it last changed. */ int m_LastApodizationArraySize; /** \brief Pointer to the GPU beamforming filter class; for performance reasons the filter is initialized within the constructor and kept for all later computations. */ mitk::PhotoacousticOCLBeamformingFilter::Pointer m_BeamformingOclFilter; }; } // namespace mitk #endif //MITK_PHOTOACOUSTICS_BEAMFORMING_FILTER diff --git a/Modules/PhotoacousticsAlgorithms/include/mitkBeamformingSettings.h b/Modules/PhotoacousticsAlgorithms/include/mitkBeamformingSettings.h index 8684d9632e..891dbf0d88 100644 --- a/Modules/PhotoacousticsAlgorithms/include/mitkBeamformingSettings.h +++ b/Modules/PhotoacousticsAlgorithms/include/mitkBeamformingSettings.h @@ -1,248 +1,242 @@ /*=================================================================== 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_BEAMFORMING_SETTINGS #define MITK_BEAMFORMING_SETTINGS #include #include #include #include namespace mitk { /*! * \brief Class holding the configuration data for the beamforming filters mitk::BeamformingFilter and mitk::PhotoacousticOCLBeamformingFilter * * A detailed description can be seen below. All parameters should be set manually for successfull beamforming. */ class MITKPHOTOACOUSTICSALGORITHMS_EXPORT BeamformingSettings : public itk::Object { public: mitkClassMacroItkParent(BeamformingSettings, itk::Object); itkCloneMacro(Self); /** \brief Available apodization functions: * - Hamming function. * - Von-Hann function. * - Box function. */ enum Apodization { Hamm, Hann, Box }; /** \brief Available beamforming algorithms: * - DAS (Delay and sum). * - DMAS (Delay multiply and sum). */ enum BeamformingAlgorithm { DMAS, DAS, sDMAS }; /** \brief Available geometries for Probes: * - Linear * - Concave */ enum ProbeGeometry { Linear, Concave}; itkGetConstMacro(PitchInMeters, float); itkGetConstMacro(SpeedOfSound, float); itkGetConstMacro(TimeSpacing, float); itkGetConstMacro(Angle, float); itkGetConstMacro(IsPhotoacousticImage, bool); itkGetConstMacro(TransducerElements, unsigned int); itkGetConstMacro(SamplesPerLine, unsigned int); itkGetConstMacro(ReconstructionLines, unsigned int); itkGetConstMacro(InputDim, const unsigned int*); itkGetConstMacro(UseGPU, bool); itkGetConstMacro(GPUBatchSize, unsigned int); itkGetConstMacro(ApodizationFunction, const float*); itkGetConstMacro(Apod, Apodization); itkGetConstMacro(ApodizationArraySize, int); itkGetConstMacro(Algorithm, BeamformingAlgorithm); itkGetConstMacro(ReconstructionDepth, float); itkGetConstMacro(Geometry, ProbeGeometry); itkGetConstMacro(ProbeRadius, float); itkGetConstMacro(ElementHeights, float*); itkGetConstMacro(ElementPositions, float*); itkGetConstMacro(HorizontalExtent, float); /** \brief function for mitk::PhotoacousticOCLBeamformingFilter to check whether buffers need to be updated * this method only checks parameters relevant for the openCL implementation */ static bool SettingsChangedOpenCL(const BeamformingSettings::Pointer lhs, const BeamformingSettings::Pointer rhs) { return !((std::abs(lhs->GetAngle() - rhs->GetAngle()) < 0.01f) && // 0.01 degree error margin (lhs->GetApod() == rhs->GetApod()) && (lhs->GetGeometry() == rhs->GetGeometry()) && (std::abs(lhs->GetProbeRadius() - rhs->GetProbeRadius()) < 0.001f) && (lhs->GetIsPhotoacousticImage() == rhs->GetIsPhotoacousticImage()) && (std::abs(lhs->GetPitchInMeters() - rhs->GetPitchInMeters()) < 0.000001f) && // 0.0001 mm error margin (lhs->GetReconstructionLines() == rhs->GetReconstructionLines()) && (lhs->GetSamplesPerLine() == rhs->GetSamplesPerLine()) && (lhs->GetReconstructionDepth() == rhs->GetReconstructionDepth()) && (std::abs(lhs->GetSpeedOfSound() - rhs->GetSpeedOfSound()) < 0.01f) && (std::abs(lhs->GetTimeSpacing() - rhs->GetTimeSpacing()) < 0.00000000001f) && //0.01 ns error margin (lhs->GetTransducerElements() == rhs->GetTransducerElements())); } static Pointer New(float pitchInMeters, float speedOfSound, float timeSpacing, float angle, bool isPhotoacousticImage, unsigned int samplesPerLine, unsigned int reconstructionLines, unsigned int* inputDim, float reconstructionDepth, bool useGPU, unsigned int GPUBatchSize, Apodization apod, unsigned int apodizationArraySize, BeamformingAlgorithm algorithm, ProbeGeometry geometry, float probeRadius) { Pointer smartPtr = new BeamformingSettings(pitchInMeters, speedOfSound, timeSpacing, angle, isPhotoacousticImage, samplesPerLine, reconstructionLines, inputDim, reconstructionDepth, useGPU, GPUBatchSize, apod, apodizationArraySize, algorithm, geometry, probeRadius); smartPtr->UnRegister(); return smartPtr; } - unsigned short* GetMinMaxLines(); - protected: /** */ BeamformingSettings(float pitchInMeters, float speedOfSound, float timeSpacing, float angle, bool isPhotoacousticImage, unsigned int samplesPerLine, unsigned int reconstructionLines, unsigned int* inputDim, float reconstructionDepth, bool useGPU, unsigned int GPUBatchSize, Apodization apod, unsigned int apodizationArraySize, BeamformingAlgorithm algorithm, ProbeGeometry geometry, float probeRadius ); ~BeamformingSettings() override; /** \brief Pitch of the used transducer in [m]. */ float m_PitchInMeters; /** \brief Speed of sound in the used medium in [m/s]. */ float m_SpeedOfSound; /** \brief The time spacing of the input image */ float m_TimeSpacing; // [s] /** \brief The angle of the transducer elements */ float m_Angle; /** \brief Flag whether processed image is a photoacoustic image or an ultrasound image */ bool m_IsPhotoacousticImage; /** \brief How many transducer elements the used transducer had. */ unsigned int m_TransducerElements; /** \brief How many vertical samples should be used in the final image. */ unsigned int m_SamplesPerLine; /** \brief How many lines should be reconstructed in the final image. */ unsigned int m_ReconstructionLines; /** \brief Sets the dimensions of the inputImage. */ const unsigned int* m_InputDim; /** \brief The Depth up to which the filter should reconstruct the image [m] */ float m_ReconstructionDepth; /** \brief Decides whether GPU computing should be used */ bool m_UseGPU; unsigned int m_GPUBatchSize; /** \brief Sets the amount of image slices in batches when GPU is used */ const float* m_ApodizationFunction; /** \brief Sets the used apodization function. */ Apodization m_Apod; /** \brief Sets the resolution of the apodization array (must be greater than 0). */ int m_ApodizationArraySize; /** \brief Sets the used beamforming algorithm. */ BeamformingAlgorithm m_Algorithm; /** \brief Sets the used probe geometry */ ProbeGeometry m_Geometry; /** \brief Sets the radius of the curved probe [m] */ float m_ProbeRadius; /** */ float *m_ElementHeights; /** */ float *m_ElementPositions; /** */ float m_HorizontalExtent; - - /** - */ - unsigned short* m_MinMaxLines; }; } #endif //MITK_BEAMFORMING_SETTINGS diff --git a/Modules/PhotoacousticsAlgorithms/include/mitkBeamformingUtils.h b/Modules/PhotoacousticsAlgorithms/include/mitkBeamformingUtils.h index fff7126538..6fe1eecc09 100644 --- a/Modules/PhotoacousticsAlgorithms/include/mitkBeamformingUtils.h +++ b/Modules/PhotoacousticsAlgorithms/include/mitkBeamformingUtils.h @@ -1,72 +1,68 @@ /*=================================================================== 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_BEAMFORMING_FILTER_UTILS #define MITK_BEAMFORMING_FILTER_UTILS #include "mitkImageToImageFilter.h" #include #include "./OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.h" #include "mitkBeamformingSettings.h" namespace mitk { /*! * \brief Class implementing util functionality for beamforming on CPU * */ class BeamformingUtils final { public: /** \brief Function to perform beamforming on CPU for a single line, using DAS and spherical delay */ static void DASSphericalLine(float* input, float* output, float inputDim[2], float outputDim[2], const short& line, const mitk::BeamformingSettings::Pointer config); /** \brief Function to perform beamforming on CPU for a single line, using DMAS and spherical delay */ static void DMASSphericalLine(float* input, float* output, float inputDim[2], float outputDim[2], const short& line, const mitk::BeamformingSettings::Pointer config); /** \brief Function to perform beamforming on CPU for a single line, using signed DMAS and spherical delay */ static void sDMASSphericalLine(float* input, float* output, float inputDim[2], float outputDim[2], const short& line, const mitk::BeamformingSettings::Pointer config); /** \brief Pointer holding the Von-Hann apodization window for beamforming * @param samples the resolution at which the window is created */ static float* VonHannFunction(int samples); /** \brief Function to create a Hamming apodization window * @param samples the resolution at which the window is created */ static float* HammFunction(int samples); /** \brief Function to create a Box apodization window * @param samples the resolution at which the window is created */ static float* BoxFunction(int samples); - /** \brief - */ - static unsigned short* MinMaxLines(const mitk::BeamformingSettings::Pointer config); - protected: BeamformingUtils(); ~BeamformingUtils(); }; } // namespace mitk #endif //MITK_BEAMFORMING_FILTER_UTILS diff --git a/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.cpp b/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.cpp index dc3ad02f83..ba58a646ef 100644 --- a/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.cpp +++ b/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.cpp @@ -1,325 +1,314 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #if defined(PHOTOACOUSTICS_USE_GPU) || DOXYGEN #include "./OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.h" #include "usServiceReference.h" mitk::PhotoacousticOCLBeamformingFilter::PhotoacousticOCLBeamformingFilter(BeamformingSettings::Pointer settings) : m_PixelCalculation(NULL), m_inputSlices(1), m_Conf(settings), m_InputImage(mitk::Image::New()), m_ApodizationBuffer(nullptr), m_DelaysBuffer(nullptr), - m_UsedLinesBuffer(nullptr), m_ElementHeightsBuffer(nullptr), m_ElementPositionsBuffer(nullptr) { MITK_INFO << "Instantiating OCL beamforming Filter..."; this->AddSourceFile("DAS.cl"); this->AddSourceFile("DMAS.cl"); this->AddSourceFile("sDMAS.cl"); this->m_FilterID = "OpenCLBeamformingFilter"; this->Initialize(); unsigned int dim[] = { 128, 2048, 2 }; m_InputImage->Initialize(mitk::MakeScalarPixelType(), 3, dim); m_ChunkSize[0] = 128; m_ChunkSize[1] = 128; m_ChunkSize[2] = 8; - m_UsedLinesCalculation = mitk::OCLUsedLinesCalculation::New(m_Conf); m_DelayCalculation = mitk::OCLDelayCalculation::New(m_Conf); MITK_INFO << "Instantiating OCL beamforming Filter...[Done]"; } mitk::PhotoacousticOCLBeamformingFilter::~PhotoacousticOCLBeamformingFilter() { if (this->m_PixelCalculation) { clReleaseKernel(m_PixelCalculation); } if (m_ApodizationBuffer) clReleaseMemObject(m_ApodizationBuffer); if (m_ElementHeightsBuffer) clReleaseMemObject(m_ElementHeightsBuffer); if (m_ElementPositionsBuffer) clReleaseMemObject(m_ElementPositionsBuffer); } void mitk::PhotoacousticOCLBeamformingFilter::Update() { //Check if context & program available if (!this->Initialize()) { us::ServiceReference ref = GetModuleContext()->GetServiceReference(); OclResourceService* resources = GetModuleContext()->GetService(ref); // clean-up also the resources resources->InvalidateStorage(); mitkThrow() << "Filter is not initialized. Cannot update."; } else { // Execute this->Execute(); } } void mitk::PhotoacousticOCLBeamformingFilter::UpdateDataBuffers() { us::ServiceReference ref = GetModuleContext()->GetServiceReference(); OclResourceService* resources = GetModuleContext()->GetService(ref); cl_context gpuContext = resources->GetContext(); //Initialize the Output try { size_t outputSize = (size_t)m_Conf->GetReconstructionLines() * (size_t)m_Conf->GetSamplesPerLine() * (size_t)m_inputSlices; m_OutputDim[0] = m_Conf->GetReconstructionLines(); m_OutputDim[1] = m_Conf->GetSamplesPerLine(); m_OutputDim[2] = m_inputSlices; this->InitExec(this->m_PixelCalculation, m_OutputDim, outputSize, sizeof(float)); } catch (const mitk::Exception& e) { MITK_ERROR << "Caught exception while initializing filter: " << e.what(); return; } cl_int clErr = 0; MITK_DEBUG << "Updating GPU Buffers for new configuration"; // create the apodisation buffer if (m_Apodisation == nullptr) { MITK_INFO << "No apodisation function set; Beamforming will be done without any apodisation."; m_Apodisation = new float[1]{ 1 }; m_ApodArraySize = 1; } if (m_ApodizationBuffer) clReleaseMemObject(m_ApodizationBuffer); this->m_ApodizationBuffer = clCreateBuffer(gpuContext, CL_MEM_READ_ONLY | CL_MEM_USE_HOST_PTR, sizeof(float) * m_ApodArraySize, const_cast(m_Apodisation), &clErr); CHECK_OCL_ERR(clErr); if (m_ElementHeightsBuffer) clReleaseMemObject(m_ElementHeightsBuffer); this->m_ElementHeightsBuffer = clCreateBuffer(gpuContext, CL_MEM_READ_ONLY | CL_MEM_USE_HOST_PTR, sizeof(float) * m_Conf->GetTransducerElements(), const_cast(m_Conf->GetElementHeights()), &clErr); CHECK_OCL_ERR(clErr); if (m_ElementPositionsBuffer) clReleaseMemObject(m_ElementPositionsBuffer); this->m_ElementPositionsBuffer = clCreateBuffer(gpuContext, CL_MEM_READ_ONLY | CL_MEM_USE_HOST_PTR, sizeof(float) * m_Conf->GetTransducerElements(), const_cast(m_Conf->GetElementPositions()), &clErr); CHECK_OCL_ERR(clErr); - // calculate used lines - - m_UsedLinesCalculation->SetElementPositionsBuffer(m_ElementPositionsBuffer); - m_UsedLinesCalculation->SetElementHeightsBuffer(m_ElementHeightsBuffer); - m_UsedLinesCalculation->Update(); - m_UsedLinesBuffer = m_UsedLinesCalculation->GetGPUOutput()->GetGPUBuffer(); // calculate the Delays - m_DelayCalculation->SetInputs(m_UsedLinesBuffer); m_DelayCalculation->Update(); m_DelaysBuffer = m_DelayCalculation->GetGPUOutput()->GetGPUBuffer(); } void mitk::PhotoacousticOCLBeamformingFilter::Execute() { cl_int clErr = 0; UpdateDataBuffers(); if (m_Conf->GetGeometry() == mitk::BeamformingSettings::ProbeGeometry::Linear) { unsigned int reconstructionLines = this->m_Conf->GetReconstructionLines(); unsigned int samplesPerLine = this->m_Conf->GetSamplesPerLine(); - clErr = clSetKernelArg(this->m_PixelCalculation, 2, sizeof(cl_mem), &(this->m_UsedLinesBuffer)); - clErr |= clSetKernelArg(this->m_PixelCalculation, 3, sizeof(cl_mem), &(this->m_DelaysBuffer)); - clErr |= clSetKernelArg(this->m_PixelCalculation, 4, sizeof(cl_mem), &(this->m_ApodizationBuffer)); - clErr |= clSetKernelArg(this->m_PixelCalculation, 5, sizeof(cl_ushort), &(this->m_ApodArraySize)); - clErr |= clSetKernelArg(this->m_PixelCalculation, 6, sizeof(cl_uint), &(this->m_Conf->GetInputDim()[0])); - clErr |= clSetKernelArg(this->m_PixelCalculation, 7, sizeof(cl_uint), &(this->m_Conf->GetInputDim()[1])); - clErr |= clSetKernelArg(this->m_PixelCalculation, 8, sizeof(cl_uint), &(m_inputSlices)); - clErr |= clSetKernelArg(this->m_PixelCalculation, 9, sizeof(cl_uint), &(reconstructionLines)); - clErr |= clSetKernelArg(this->m_PixelCalculation, 10, sizeof(cl_uint), &(samplesPerLine)); + clErr = clSetKernelArg(this->m_PixelCalculation, 2, sizeof(cl_mem), &(this->m_DelaysBuffer)); + clErr |= clSetKernelArg(this->m_PixelCalculation, 3, sizeof(cl_mem), &(this->m_ApodizationBuffer)); + clErr |= clSetKernelArg(this->m_PixelCalculation, 4, sizeof(cl_ushort), &(this->m_ApodArraySize)); + clErr |= clSetKernelArg(this->m_PixelCalculation, 5, sizeof(cl_uint), &(this->m_Conf->GetInputDim()[0])); + clErr |= clSetKernelArg(this->m_PixelCalculation, 6, sizeof(cl_uint), &(this->m_Conf->GetInputDim()[1])); + clErr |= clSetKernelArg(this->m_PixelCalculation, 7, sizeof(cl_uint), &(m_inputSlices)); + clErr |= clSetKernelArg(this->m_PixelCalculation, 8, sizeof(cl_uint), &(reconstructionLines)); + clErr |= clSetKernelArg(this->m_PixelCalculation, 9, sizeof(cl_uint), &(samplesPerLine)); } else { int reconstructionLines = this->m_Conf->GetReconstructionLines(); int samplesPerLine = this->m_Conf->GetSamplesPerLine(); float totalSamples_i = (float)(m_Conf->GetReconstructionDepth()) / (float)(m_Conf->GetSpeedOfSound() * m_Conf->GetTimeSpacing()); totalSamples_i = totalSamples_i <= m_Conf->GetInputDim()[1] ? totalSamples_i : m_Conf->GetInputDim()[1]; float horizontalExtent = m_Conf->GetHorizontalExtent(); float mult = 1 / (this->m_Conf->GetTimeSpacing() * this->m_Conf->GetSpeedOfSound()); char isPAImage = (char)m_Conf->GetIsPhotoacousticImage(); clErr = clSetKernelArg(this->m_PixelCalculation, 2, sizeof(cl_mem), &(this->m_ElementHeightsBuffer)); clErr |= clSetKernelArg(this->m_PixelCalculation, 3, sizeof(cl_mem), &(this->m_ElementPositionsBuffer)); clErr |= clSetKernelArg(this->m_PixelCalculation, 4, sizeof(cl_mem), &(this->m_ApodizationBuffer)); clErr |= clSetKernelArg(this->m_PixelCalculation, 5, sizeof(cl_ushort), &(this->m_ApodArraySize)); clErr |= clSetKernelArg(this->m_PixelCalculation, 6, sizeof(cl_uint), &(this->m_Conf->GetInputDim()[0])); clErr |= clSetKernelArg(this->m_PixelCalculation, 7, sizeof(cl_uint), &(this->m_Conf->GetInputDim()[1])); clErr |= clSetKernelArg(this->m_PixelCalculation, 8, sizeof(cl_int), &(m_inputSlices)); clErr |= clSetKernelArg(this->m_PixelCalculation, 9, sizeof(cl_int), &(reconstructionLines)); clErr |= clSetKernelArg(this->m_PixelCalculation, 10, sizeof(cl_int), &(samplesPerLine)); clErr |= clSetKernelArg(this->m_PixelCalculation, 11, sizeof(cl_float), &(totalSamples_i)); clErr |= clSetKernelArg(this->m_PixelCalculation, 12, sizeof(cl_float), &(horizontalExtent)); clErr |= clSetKernelArg(this->m_PixelCalculation, 13, sizeof(cl_float), &(mult)); clErr |= clSetKernelArg(this->m_PixelCalculation, 14, sizeof(cl_char), &(isPAImage)); - clErr |= clSetKernelArg(this->m_PixelCalculation, 15, sizeof(cl_mem), &(this->m_UsedLinesBuffer)); } // execute the filter on a 2D/3D NDRange if (m_OutputDim[2] == 1 || m_ChunkSize[2] == 1) { if (!this->ExecuteKernelChunksInBatches(m_PixelCalculation, 2, m_ChunkSize, m_inputSlices, 50)) mitkThrow() << "openCL Error when executing Kernel"; } else { if (!this->ExecuteKernelChunksInBatches(m_PixelCalculation, 3, m_ChunkSize, m_inputSlices, 50)) mitkThrow() << "openCL Error when executing Kernel"; } // signalize the GPU-side data changed m_Output->Modified(GPU_DATA); } us::Module *mitk::PhotoacousticOCLBeamformingFilter::GetModule() { return us::GetModuleContext()->GetModule(); } bool mitk::PhotoacousticOCLBeamformingFilter::Initialize() { bool buildErr = true; cl_int clErr = 0; if (OclFilter::Initialize()) { if (m_Conf->GetGeometry() == mitk::BeamformingSettings::ProbeGeometry::Linear) { switch (m_Conf->GetAlgorithm()) { case BeamformingSettings::BeamformingAlgorithm::DAS: { MITK_INFO << "DAS bf"; this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "ckDAS", &clErr); break; } case BeamformingSettings::BeamformingAlgorithm::DMAS: { MITK_INFO << "DMAS bf"; this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "ckDMAS", &clErr); break; } case BeamformingSettings::BeamformingAlgorithm::sDMAS: { MITK_INFO << "sDMAS bf"; this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "cksDMAS", &clErr); break; } default: { MITK_INFO << "No beamforming algorithm specified, setting to DAS"; this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "ckDAS", &clErr); break; } } } else { switch (m_Conf->GetAlgorithm()) { case BeamformingSettings::BeamformingAlgorithm::DAS: { MITK_INFO << "DAS bf"; this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "ckDAS_g", &clErr); break; } case BeamformingSettings::BeamformingAlgorithm::DMAS: { MITK_INFO << "DMAS bf"; this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "ckDMAS_g", &clErr); break; } case BeamformingSettings::BeamformingAlgorithm::sDMAS: { MITK_INFO << "sDMAS bf"; this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "cksDMAS_g", &clErr); break; } default: { MITK_INFO << "No beamforming algorithm specified, setting to DAS"; this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "ckDAS_g", &clErr); break; } } } buildErr |= CHECK_OCL_ERR(clErr); } CHECK_OCL_ERR(clErr); return (OclFilter::IsInitialized() && buildErr); } void mitk::PhotoacousticOCLBeamformingFilter::SetInput(mitk::Image::Pointer image) { OclDataSetToDataSetFilter::SetInput(image); m_InputImage = image; m_inputSlices = image->GetDimension(2); } void mitk::PhotoacousticOCLBeamformingFilter::SetInput(void* data, unsigned int* dimensions, unsigned int BpE) { OclDataSetToDataSetFilter::SetInput(data, dimensions[0] * dimensions[1] * dimensions[2], BpE); } mitk::Image::Pointer mitk::PhotoacousticOCLBeamformingFilter::GetOutputAsImage() { mitk::Image::Pointer outputImage = mitk::Image::New(); if (m_Output->IsModified(GPU_DATA)) { void* pData = m_Output->TransferDataToCPU(m_CommandQue); const unsigned int dimension = 3; unsigned int dimensions[3] = { (unsigned int)m_OutputDim[0], (unsigned int)m_OutputDim[1], (unsigned int)m_OutputDim[2] }; const mitk::SlicedGeometry3D::Pointer p_slg = m_InputImage->GetSlicedGeometry(); MITK_DEBUG << "Creating new MITK Image."; outputImage->Initialize(this->GetOutputType(), dimension, dimensions); outputImage->SetSpacing(p_slg->GetSpacing()); outputImage->SetImportVolume(pData, 0, 0, mitk::Image::ImportMemoryManagementType::CopyMemory); free(pData); } 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 index b6f4e2f0f7..ec8eb77fb4 100644 --- a/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLDelayCalculation.cpp +++ b/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLDelayCalculation.cpp @@ -1,128 +1,127 @@ /*=================================================================== 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(mitk::BeamformingSettings::Pointer settings) : m_PixelCalculation(NULL), m_Conf(settings) { this->AddSourceFile("DelayCalculation.cl"); this->m_FilterID = "DelayCalculation"; m_ChunkSize[0] = 128; m_ChunkSize[1] = 128; m_ChunkSize[2] = 8; 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->GetReconstructionLines() / 2, m_Conf->GetSamplesPerLine(), 1 }; m_BufferSize = gridDim[0] * gridDim[1] * 1; try { this->InitExecNoInput(this->m_PixelCalculation, gridDim, m_BufferSize, sizeof(unsigned short)); } catch (const mitk::Exception& e) { MITK_ERROR << "Caught exception while initializing Delay Calculation filter: " << e.what(); return; } // This calculation is the same for all kernels, so for performance reasons simply perform it here instead of within the kernels m_DelayMultiplicatorRaw = 1 / (m_Conf->GetTimeSpacing()*m_Conf->GetSpeedOfSound()) * (m_Conf->GetPitchInMeters()*(float)m_Conf->GetTransducerElements()); // as openCL does not support bool as a kernel argument, we need to buffer this value in a char... m_IsPAImage = m_Conf->GetIsPhotoacousticImage(); unsigned int reconstructionLines = this->m_Conf->GetReconstructionLines(); unsigned int samplesperLine = this->m_Conf->GetSamplesPerLine(); float totalSamples_i = (float)(m_Conf->GetReconstructionDepth()) / (float)(m_Conf->GetSpeedOfSound() * m_Conf->GetTimeSpacing()); totalSamples_i = totalSamples_i <= m_Conf->GetInputDim()[1] ? totalSamples_i : m_Conf->GetInputDim()[1]; - clErr = clSetKernelArg(this->m_PixelCalculation, 1, sizeof(cl_mem), &(this->m_UsedLines)); - clErr |= clSetKernelArg(this->m_PixelCalculation, 2, sizeof(cl_uint), &(this->m_Conf->GetInputDim()[0])); - clErr |= clSetKernelArg(this->m_PixelCalculation, 3, sizeof(cl_uint), &(this->m_Conf->GetInputDim()[1])); - clErr |= clSetKernelArg(this->m_PixelCalculation, 4, sizeof(cl_uint), &(reconstructionLines)); - clErr |= clSetKernelArg(this->m_PixelCalculation, 5, sizeof(cl_uint), &(samplesperLine)); - clErr |= clSetKernelArg(this->m_PixelCalculation, 6, sizeof(cl_char), &(this->m_IsPAImage)); - clErr |= clSetKernelArg(this->m_PixelCalculation, 7, sizeof(cl_float), &(this->m_DelayMultiplicatorRaw)); - clErr |= clSetKernelArg(this->m_PixelCalculation, 8, sizeof(cl_float), &(totalSamples_i)); + clErr = clSetKernelArg(this->m_PixelCalculation, 1, sizeof(cl_uint), &(this->m_Conf->GetInputDim()[0])); + clErr |= clSetKernelArg(this->m_PixelCalculation, 2, sizeof(cl_uint), &(this->m_Conf->GetInputDim()[1])); + clErr |= clSetKernelArg(this->m_PixelCalculation, 3, sizeof(cl_uint), &(reconstructionLines)); + clErr |= clSetKernelArg(this->m_PixelCalculation, 4, sizeof(cl_uint), &(samplesperLine)); + clErr |= clSetKernelArg(this->m_PixelCalculation, 5, sizeof(cl_char), &(this->m_IsPAImage)); + clErr |= clSetKernelArg(this->m_PixelCalculation, 6, sizeof(cl_float), &(this->m_DelayMultiplicatorRaw)); + clErr |= clSetKernelArg(this->m_PixelCalculation, 7, sizeof(cl_float), &(totalSamples_i)); CHECK_OCL_ERR(clErr); // execute the filter on a 3D NDRange if (!this->ExecuteKernelChunksInBatches(m_PixelCalculation, 2, m_ChunkSize, 16, 50)) mitkThrow() << "openCL Error when executing Kernel"; // 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()) { this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "ckDelayCalculationSphe", &clErr); buildErr |= CHECK_OCL_ERR(clErr); } return (OclFilter::IsInitialized() && buildErr); } diff --git a/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLUsedLinesCalculation.cpp b/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLUsedLinesCalculation.cpp deleted file mode 100644 index 8da3c4c197..0000000000 --- a/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLUsedLinesCalculation.cpp +++ /dev/null @@ -1,167 +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. - -===================================================================*/ -#if defined(PHOTOACOUSTICS_USE_GPU) || DOXYGEN - -#include "./OpenCLFilter/mitkPhotoacousticOCLUsedLinesCalculation.h" -#include "usServiceReference.h" -#include "mitkImageReadAccessor.h" - -mitk::OCLUsedLinesCalculation::OCLUsedLinesCalculation(mitk::BeamformingSettings::Pointer settings) - : m_PixelCalculation(NULL), - m_Conf(settings) -{ - this->AddSourceFile("UsedLinesCalculation.cl"); - this->m_FilterID = "UsedLinesCalculation"; - - m_ChunkSize[0] = 128; - m_ChunkSize[1] = 128; - m_ChunkSize[2] = 8; - - this->Initialize(); -} - -mitk::OCLUsedLinesCalculation::~OCLUsedLinesCalculation() -{ - if (this->m_PixelCalculation) - { - clReleaseKernel(m_PixelCalculation); - } -} - -void mitk::OCLUsedLinesCalculation::SetElementHeightsBuffer(cl_mem elementHeightsBuffer) -{ - m_ElementHeightsBuffer = elementHeightsBuffer; -} - -void mitk::OCLUsedLinesCalculation::SetElementPositionsBuffer(cl_mem elementPositionsBuffer) -{ - m_ElementPositionsBuffer = elementPositionsBuffer; -} - -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->GetReconstructionLines(), m_Conf->GetSamplesPerLine(), 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 UsedLines filter: " << e.what(); - return; - } - - if (m_Conf->GetGeometry() == mitk::BeamformingSettings::ProbeGeometry::Linear) - { - // This calculation is the same for all kernels, so for performance reasons simply perform it here instead of within the kernels - m_part = (tan(m_Conf->GetAngle() / 360 * 2 * itk::Math::pi) * - ((m_Conf->GetSpeedOfSound() * m_Conf->GetTimeSpacing())) / - (m_Conf->GetPitchInMeters() * m_Conf->GetTransducerElements())) * m_Conf->GetInputDim()[0]; - - unsigned int reconLines = this->m_Conf->GetReconstructionLines(); - unsigned int samplesPerLine = this->m_Conf->GetSamplesPerLine(); - - float totalSamples_i = (float)(m_Conf->GetReconstructionDepth()) / (float)(m_Conf->GetSpeedOfSound() * m_Conf->GetTimeSpacing()); - totalSamples_i = totalSamples_i <= m_Conf->GetInputDim()[1] ? totalSamples_i : m_Conf->GetInputDim()[1]; - - clErr = clSetKernelArg(this->m_PixelCalculation, 1, sizeof(cl_float), &(this->m_part)); - clErr |= clSetKernelArg(this->m_PixelCalculation, 2, sizeof(cl_uint), &(this->m_Conf->GetInputDim()[0])); - clErr |= clSetKernelArg(this->m_PixelCalculation, 3, sizeof(cl_uint), &(this->m_Conf->GetInputDim()[1])); - clErr |= clSetKernelArg(this->m_PixelCalculation, 4, sizeof(cl_uint), &(reconLines)); - clErr |= clSetKernelArg(this->m_PixelCalculation, 5, sizeof(cl_uint), &(samplesPerLine)); - clErr |= clSetKernelArg(this->m_PixelCalculation, 6, sizeof(cl_float), &(totalSamples_i)); - } - else - { - unsigned int reconLines = this->m_Conf->GetReconstructionLines(); - unsigned int samplesPerLine = this->m_Conf->GetSamplesPerLine(); - - float probeRadius = m_Conf->GetProbeRadius(); - float cos_deg = std::cos(m_Conf->GetAngle()/2.f / 360 * 2 * itk::Math::pi); - - float horizontalExtent = m_Conf->GetHorizontalExtent(); - float verticalExtent = m_Conf->GetReconstructionDepth(); - - clErr = clSetKernelArg(this->m_PixelCalculation, 1, sizeof(cl_mem), &(this->m_ElementHeightsBuffer)); - clErr |= clSetKernelArg(this->m_PixelCalculation, 2, sizeof(cl_mem), &(this->m_ElementPositionsBuffer)); - clErr |= clSetKernelArg(this->m_PixelCalculation, 3, sizeof(cl_float), &(cos_deg)); - clErr |= clSetKernelArg(this->m_PixelCalculation, 4, sizeof(cl_float), &(probeRadius)); - clErr |= clSetKernelArg(this->m_PixelCalculation, 5, sizeof(cl_uint), &(this->m_Conf->GetInputDim()[0])); - clErr |= clSetKernelArg(this->m_PixelCalculation, 6, sizeof(cl_uint), &(this->m_Conf->GetInputDim()[1])); - clErr |= clSetKernelArg(this->m_PixelCalculation, 7, sizeof(cl_uint), &(reconLines)); - clErr |= clSetKernelArg(this->m_PixelCalculation, 8, sizeof(cl_uint), &(samplesPerLine)); - clErr |= clSetKernelArg(this->m_PixelCalculation, 9, sizeof(cl_float), &(horizontalExtent)); - clErr |= clSetKernelArg(this->m_PixelCalculation, 10, sizeof(cl_float), &(verticalExtent)); - } - - CHECK_OCL_ERR(clErr); - - // execute the filter on a 2D NDRange - if (!this->ExecuteKernelChunksInBatches(m_PixelCalculation, 2, m_ChunkSize, 16, 50)) - mitkThrow() << "openCL Error when executing Kernel"; - - // 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()) - { - if (m_Conf->GetGeometry() == mitk::BeamformingSettings::ProbeGeometry::Linear) - { - this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "ckUsedLines", &clErr); - buildErr |= CHECK_OCL_ERR(clErr); - } - else - { - this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "ckUsedLines_g", &clErr); - buildErr |= CHECK_OCL_ERR(clErr); - } - } - return (OclFilter::IsInitialized() && buildErr); -} -#endif diff --git a/Modules/PhotoacousticsAlgorithms/source/filters/mitkBeamformingSettings.cpp b/Modules/PhotoacousticsAlgorithms/source/filters/mitkBeamformingSettings.cpp index 894802fe76..223778796c 100644 --- a/Modules/PhotoacousticsAlgorithms/source/filters/mitkBeamformingSettings.cpp +++ b/Modules/PhotoacousticsAlgorithms/source/filters/mitkBeamformingSettings.cpp @@ -1,141 +1,131 @@ /*=================================================================== mitkBeamformingSettings The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkBeamformingSettings.h" #include "mitkBeamformingUtils.h" #include "itkMutexLock.h" mitk::BeamformingSettings::BeamformingSettings(float pitchInMeters, float speedOfSound, float timeSpacing, float angle, bool isPhotoacousticImage, unsigned int samplesPerLine, unsigned int reconstructionLines, unsigned int* inputDim, float reconstructionDepth, bool useGPU, unsigned int GPUBatchSize, Apodization apod, unsigned int apodizationArraySize, BeamformingAlgorithm algorithm, ProbeGeometry geometry, float probeRadius ) : m_PitchInMeters(pitchInMeters), m_SpeedOfSound(speedOfSound), m_TimeSpacing(timeSpacing), m_Angle(angle), m_IsPhotoacousticImage(isPhotoacousticImage), m_SamplesPerLine(samplesPerLine), m_ReconstructionLines(reconstructionLines), m_ReconstructionDepth(reconstructionDepth), m_UseGPU(useGPU), m_GPUBatchSize(GPUBatchSize), m_Apod(apod), m_ApodizationArraySize(apodizationArraySize), m_Algorithm(algorithm), m_Geometry(geometry), - m_ProbeRadius(probeRadius), - m_MinMaxLines(nullptr) + m_ProbeRadius(probeRadius) { if (inputDim == nullptr) { MITK_ERROR << "No input dimension given."; mitkThrow() << "No input dimension given."; } switch (GetApod()) { case BeamformingSettings::Apodization::Hann: m_ApodizationFunction = mitk::BeamformingUtils::VonHannFunction(GetApodizationArraySize()); break; case BeamformingSettings::Apodization::Hamm: m_ApodizationFunction = mitk::BeamformingUtils::HammFunction(GetApodizationArraySize()); break; case BeamformingSettings::Apodization::Box: default: m_ApodizationFunction = mitk::BeamformingUtils::BoxFunction(GetApodizationArraySize()); break; } m_InputDim = new unsigned int[3]{ inputDim[0], inputDim[1], inputDim[2] }; m_TransducerElements = m_InputDim[0]; m_ElementHeights = new float[m_TransducerElements]; m_ElementPositions = new float[m_TransducerElements]; if (m_Geometry == ProbeGeometry::Concave) { float openingAngle = (m_TransducerElements * m_PitchInMeters) / (probeRadius * 2 * itk::Math::pi) * 2 * itk::Math::pi; m_HorizontalExtent = std::sin(openingAngle / 2.f) * probeRadius * 2.f; float elementAngle = 0; for (unsigned int i = 0; i < m_TransducerElements; ++i) { elementAngle = ((i- m_TransducerElements /2.f) * m_PitchInMeters) / (probeRadius * 2 * itk::Math::pi) * 2 * itk::Math::pi; m_ElementHeights[i] = probeRadius - std::cos(elementAngle) * probeRadius; m_ElementPositions[i] = m_HorizontalExtent/2.f + std::sin(elementAngle) * probeRadius; } } else { m_HorizontalExtent = m_PitchInMeters * m_TransducerElements; for (unsigned int i = 0; i < m_TransducerElements; ++i) { m_ElementHeights[i] = 0; m_ElementPositions[i] = i * m_PitchInMeters; } } } mitk::BeamformingSettings::~BeamformingSettings() { MITK_INFO << "Destructing beamforming settings..."; //Free memory if (m_ApodizationFunction != nullptr) { MITK_INFO << "Deleting apodization function..."; delete[] m_ApodizationFunction; MITK_INFO << "Deleting apodization function...[Done]"; } if (m_InputDim != nullptr) { MITK_INFO << "Deleting input dim..."; delete[] m_InputDim; MITK_INFO << "Deleting input dim...[Done]"; } if (m_ElementHeights != nullptr || m_ElementPositions != nullptr) { MITK_INFO << "Deleting element geometry..."; if (m_ElementHeights != nullptr) delete[] m_ElementHeights; if (m_ElementPositions != nullptr) delete[] m_ElementPositions; MITK_INFO << "Destructing beamforming settings...[Done]"; } - if (m_MinMaxLines) - delete[] m_MinMaxLines; } - -unsigned short* mitk::BeamformingSettings::GetMinMaxLines() -{ - if (!m_MinMaxLines) - m_MinMaxLines = mitk::BeamformingUtils::MinMaxLines(this); - return m_MinMaxLines; -} \ No newline at end of file diff --git a/Modules/PhotoacousticsAlgorithms/source/utils/mitkBeamformingUtils.cpp b/Modules/PhotoacousticsAlgorithms/source/utils/mitkBeamformingUtils.cpp index b367bb57f3..854473ea8f 100644 --- a/Modules/PhotoacousticsAlgorithms/source/utils/mitkBeamformingUtils.cpp +++ b/Modules/PhotoacousticsAlgorithms/source/utils/mitkBeamformingUtils.cpp @@ -1,420 +1,277 @@ /*=================================================================== 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. ===================================================================*/ #include "mitkProperties.h" #include "mitkImageReadAccessor.h" #include #include #include #include #include #include "mitkImageCast.h" #include "mitkBeamformingUtils.h" mitk::BeamformingUtils::BeamformingUtils() { } mitk::BeamformingUtils::~BeamformingUtils() { } float* mitk::BeamformingUtils::VonHannFunction(int samples) { float* ApodWindow = new float[samples]; for (int n = 0; n < samples; ++n) { ApodWindow[n] = (1 - cos(2 * itk::Math::pi * n / (samples - 1))) / 2; } return ApodWindow; } float* mitk::BeamformingUtils::HammFunction(int samples) { float* ApodWindow = new float[samples]; for (int n = 0; n < samples; ++n) { ApodWindow[n] = 0.54 - 0.46*cos(2 * itk::Math::pi*n / (samples - 1)); } return ApodWindow; } float* mitk::BeamformingUtils::BoxFunction(int samples) { float* ApodWindow = new float[samples]; for (int n = 0; n < samples; ++n) { ApodWindow[n] = 1; } return ApodWindow; } -unsigned short* mitk::BeamformingUtils::MinMaxLines(const mitk::BeamformingSettings::Pointer config) -{ - int outputL = (int)config->GetReconstructionLines(); - int outputS = (int)config->GetSamplesPerLine(); - - unsigned short* dDest = new unsigned short[outputL * outputS * 2]; - - int inputL = (int)config->GetInputDim()[0]; - - float horizontalExtent = config->GetHorizontalExtent(); - float verticalExtent = config->GetReconstructionDepth(); - - float probeRadius = config->GetProbeRadius(); - float* elementHeights = config->GetElementHeights(); - float* elementPositions = config->GetElementPositions(); - - float cos_deg = std::cos(config->GetAngle() / 2.f / 360 * 2 * itk::Math::pi); - - float cos = 0; - float a = 0; - float d = 0; - - for (int x = 0; x < outputL; ++x) - { - for (int y = 0; y < outputS; ++y) - { - float l_p = (float)x / outputL * horizontalExtent; - float s_p = (float)y / (float)outputS * verticalExtent; - - int maxLine = inputL; - int minLine = 0; - - for (int l_s = 0; l_s < inputL; l_s += 32) - { - a = sqrt((probeRadius - s_p)*(probeRadius - s_p) + (l_p - horizontalExtent / 2)*(l_p - horizontalExtent / 2)); - d = sqrt((s_p - elementHeights[l_s])*(s_p - elementHeights[l_s]) + (l_p - elementPositions[l_s])*(l_p - elementPositions[l_s])); - cos = (d*d + probeRadius * probeRadius - a * a) / (2 * probeRadius*d); - - if (cos > cos_deg) - { - minLine = l_s - 32; - if (minLine < 0) - minLine = 0; - break; - } - } - for (int l_s = minLine; l_s < inputL; l_s += 8) - { - a = sqrt((probeRadius - s_p)*(probeRadius - s_p) + (l_p - horizontalExtent / 2)*(l_p - horizontalExtent / 2)); - d = sqrt((s_p - elementHeights[l_s])*(s_p - elementHeights[l_s]) + (l_p - elementPositions[l_s])*(l_p - elementPositions[l_s])); - cos = (d*d + probeRadius * probeRadius - a * a) / (2 * probeRadius*d); - - if (cos > cos_deg) - { - minLine = l_s - 8; - if (minLine < 0) - minLine = 0; - break; - } - } - for (int l_s = minLine; l_s < inputL; l_s += 1) - { - a = sqrt((probeRadius - s_p)*(probeRadius - s_p) + (l_p - horizontalExtent / 2)*(l_p - horizontalExtent / 2)); - d = sqrt((s_p - elementHeights[l_s])*(s_p - elementHeights[l_s]) + (l_p - elementPositions[l_s])*(l_p - elementPositions[l_s])); - cos = (d*d + probeRadius * probeRadius - a * a) / (2 * probeRadius*d); - - if (cos > cos_deg) - { - minLine = l_s; - break; - } - } - - for (int l_s = inputL; l_s >= 0; l_s -= 32) - { - a = sqrt((probeRadius - s_p)*(probeRadius - s_p) + (l_p - horizontalExtent / 2)*(l_p - horizontalExtent / 2)); - d = sqrt((s_p - elementHeights[l_s])*(s_p - elementHeights[l_s]) + (l_p - elementPositions[l_s])*(l_p - elementPositions[l_s])); - cos = (d*d + probeRadius * probeRadius - a * a) / (2 * probeRadius*d); - cos = 0; - - if (cos > cos_deg) - { - maxLine = l_s + 32; - if (maxLine > inputL) - minLine = inputL; - break; - } - } - for (int l_s = maxLine; l_s >= 0; l_s -= 8) - { - a = sqrt((probeRadius - s_p)*(probeRadius - s_p) + (l_p - horizontalExtent / 2)*(l_p - horizontalExtent / 2)); - d = sqrt((s_p - elementHeights[l_s])*(s_p - elementHeights[l_s]) + (l_p - elementPositions[l_s])*(l_p - elementPositions[l_s])); - cos = (d*d + probeRadius * probeRadius - a * a) / (2 * probeRadius*d); - cos = 0; - - if (cos > cos_deg) - { - maxLine = l_s + 8; - if (maxLine > inputL) - minLine = inputL; - break; - } - } - for (int l_s = maxLine; l_s >= 0; l_s -= 1) - { - a = sqrt((probeRadius - s_p)*(probeRadius - s_p) + (l_p - horizontalExtent / 2)*(l_p - horizontalExtent / 2)); - d = sqrt((s_p - elementHeights[l_s])*(s_p - elementHeights[l_s]) + (l_p - elementPositions[l_s])*(l_p - elementPositions[l_s])); - cos = (d*d + probeRadius * probeRadius - a * a) / (2 * probeRadius*d); - cos = 0; - - if (cos > cos_deg) - { - maxLine = l_s; - break; - } - } - - dDest[y * 2 * outputL + 2 * x] = (unsigned short)minLine; //minLine - dDest[y * 2 * outputL + 2 * x + 1] = (unsigned short)maxLine; //maxLine - } - } - return dDest; -} - void mitk::BeamformingUtils::DASSphericalLine( float* input, float* output, float inputDim[2], float outputDim[2], const short& line, const mitk::BeamformingSettings::Pointer config) { const float* apodisation = config->GetApodizationFunction(); const short apodArraySize = config->GetApodizationArraySize(); - + MITK_INFO << "=============== apodisation = " << apodisation; + MITK_INFO << "=============== apodArraySize = " << apodArraySize; const float* elementHeights = config->GetElementHeights(); const float* elementPositions = config->GetElementPositions(); 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_p = 0; float s_i = 0; + short usedLines = inputL; float apod_mult = 1; - short usedLines = (maxLine - minLine); - float totalSamples_i = (float)(config->GetReconstructionDepth()) / (float)(config->GetSpeedOfSound() * config->GetTimeSpacing()); totalSamples_i = totalSamples_i <= inputS ? totalSamples_i : inputS; l_p = (float)line / outputL * config->GetHorizontalExtent(); for (short sample = 0; sample < outputS; ++sample) { s_i = (float)sample / outputS * totalSamples_i; - - minLine = config->GetMinMaxLines()[2*sample*(short)outputL + 2*line]; - maxLine = config->GetMinMaxLines()[2*sample*(short)outputL + 2*line + 1]; - usedLines = (maxLine - minLine); - + usedLines = inputL; apod_mult = (float)apodArraySize / (float)usedLines; - for (short l_s = minLine; l_s < maxLine; ++l_s) + for (short l_s = 0; l_s < inputL; ++l_s) { AddSample = (int)sqrt( pow(s_i-elementHeights[l_s]/(config->GetSpeedOfSound()*config->GetTimeSpacing()), 2) + pow((1 / (config->GetTimeSpacing()*config->GetSpeedOfSound())) * (l_p - elementPositions[l_s]), 2) ) + (1 - config->GetIsPhotoacousticImage())*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)]; + apodisation[(short)((l_s)*apod_mult)]; else --usedLines; } output[sample*(short)outputL + line] = output[sample*(short)outputL + line] / usedLines; } } void mitk::BeamformingUtils::DMASSphericalLine( float* input, float* output, float inputDim[2], float outputDim[2], const short& line, const mitk::BeamformingSettings::Pointer config) { const float* apodisation = config->GetApodizationFunction(); const short apodArraySize = config->GetApodizationArraySize(); const float* elementHeights = config->GetElementHeights(); const float* elementPositions = config->GetElementPositions(); float& inputS = inputDim[1]; float& inputL = inputDim[0]; float& outputS = outputDim[1]; float& outputL = outputDim[0]; - short maxLine = 0; - short minLine = 0; float l_p = 0; float s_i = 0; float apod_mult = 1; - float mult = 0; - short usedLines = (maxLine - minLine); + short usedLines = inputL; float totalSamples_i = (float)(config->GetReconstructionDepth()) / (float)(config->GetSpeedOfSound() * config->GetTimeSpacing()); totalSamples_i = totalSamples_i <= inputS ? totalSamples_i : inputS; l_p = (float)line / outputL * config->GetHorizontalExtent(); for (short sample = 0; sample < outputS; ++sample) { s_i = (float)sample / outputS * totalSamples_i; - - minLine = config->GetMinMaxLines()[2 * sample*(short)outputL + 2 * line]; - maxLine = config->GetMinMaxLines()[2 * sample*(short)outputL + 2 * line + 1]; - usedLines = (maxLine - minLine); - + usedLines = inputL; 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) + short* AddSample = new short[usedLines]; + for (short l_s = 0; l_s < usedLines; ++l_s) { AddSample[l_s] = (int)sqrt( - pow(s_i - elementHeights[l_s + minLine] / (config->GetSpeedOfSound()*config->GetTimeSpacing()), 2) + pow(s_i - elementHeights[l_s] / (config->GetSpeedOfSound()*config->GetTimeSpacing()), 2) + - pow((1 / (config->GetTimeSpacing()*config->GetSpeedOfSound())) * (l_p - elementPositions[l_s + minLine]), 2) + pow((1 / (config->GetTimeSpacing()*config->GetSpeedOfSound())) * (l_p - elementPositions[l_s]), 2) ) + (1 - config->GetIsPhotoacousticImage())*s_i; } float s_1 = 0; float s_2 = 0; - for (short l_s1 = minLine; l_s1 < maxLine - 1; ++l_s1) + for (short l_s1 = 0; l_s1 < usedLines - 1; ++l_s1) { - if (AddSample[l_s1 - minLine] < inputS && AddSample[l_s1 - minLine] >= 0) + if (AddSample[l_s1] < inputS && AddSample[l_s1] >= 0) { - for (short l_s2 = l_s1 + 1; l_s2 < maxLine; ++l_s2) + for (short l_s2 = l_s1 + 1; l_s2 < inputL; ++l_s2) { - if (AddSample[l_s2 - minLine] < inputS && AddSample[l_s2 - minLine] >= 0) + if (AddSample[l_s2] < inputS && AddSample[l_s2] >= 0) { - s_2 = input[l_s2 + AddSample[l_s2 - minLine] * (short)inputL]; - s_1 = input[l_s1 + AddSample[l_s1 - minLine] * (short)inputL]; + s_2 = input[l_s2 + AddSample[l_s2] * (short)inputL]; + s_1 = input[l_s1 + AddSample[l_s1] * (short)inputL]; - mult = s_2 * apodisation[(int)((l_s2 - minLine)*apod_mult)] * s_1 * apodisation[(int)((l_s1 - minLine)*apod_mult)]; + mult = s_2 * apodisation[(int)((l_s2)*apod_mult)] * s_1 * apodisation[(int)((l_s1)*apod_mult)]; output[sample*(short)outputL + line] += sqrt(fabs(mult)) * ((mult > 0) - (mult < 0)); } } } else --usedLines; } output[sample*(short)outputL + line] = output[sample*(short)outputL + line] / (float)(pow(usedLines, 2) - (usedLines - 1)); delete[] AddSample; } } void mitk::BeamformingUtils::sDMASSphericalLine( float* input, float* output, float inputDim[2], float outputDim[2], const short& line, const mitk::BeamformingSettings::Pointer config) { const float* apodisation = config->GetApodizationFunction(); const short apodArraySize = config->GetApodizationArraySize(); const float* elementHeights = config->GetElementHeights(); const float* elementPositions = config->GetElementPositions(); float& inputS = inputDim[1]; float& inputL = inputDim[0]; float& outputS = outputDim[1]; float& outputL = outputDim[0]; - short maxLine = 0; - short minLine = 0; float l_p = 0; float s_i = 0; float apod_mult = 1; float mult = 0; - short usedLines = (maxLine - minLine); + short usedLines = inputL; float totalSamples_i = (float)(config->GetReconstructionDepth()) / (float)(config->GetSpeedOfSound() * config->GetTimeSpacing()); totalSamples_i = totalSamples_i <= inputS ? totalSamples_i : inputS; l_p = (float)line / outputL * config->GetHorizontalExtent(); for (short sample = 0; sample < outputS; ++sample) { s_i = (float)sample / outputS * totalSamples_i; - - minLine = config->GetMinMaxLines()[2 * sample*(short)outputL + 2 * line]; - maxLine = config->GetMinMaxLines()[2 * sample*(short)outputL + 2 * line + 1]; - usedLines = (maxLine - minLine); - + usedLines = inputL; 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) + short* AddSample = new short[usedLines]; + for (short l_s = 0; l_s < usedLines; ++l_s) { AddSample[l_s] = (int)sqrt( - pow(s_i - elementHeights[l_s + minLine] / (config->GetSpeedOfSound()*config->GetTimeSpacing()), 2) + pow(s_i - elementHeights[l_s] / (config->GetSpeedOfSound()*config->GetTimeSpacing()), 2) + - pow((1 / (config->GetTimeSpacing()*config->GetSpeedOfSound())) * (l_p - elementPositions[l_s + minLine]), 2) + pow((1 / (config->GetTimeSpacing()*config->GetSpeedOfSound())) * (l_p - elementPositions[l_s]), 2) ) + (1 - config->GetIsPhotoacousticImage())*s_i; } float s_1 = 0; float s_2 = 0; float sign = 0; - for (short l_s1 = minLine; l_s1 < maxLine - 1; ++l_s1) + for (short l_s1 = 0; l_s1 < inputL - 1; ++l_s1) { - if (AddSample[l_s1 - minLine] < inputS && AddSample[l_s1 - minLine] >= 0) + if (AddSample[l_s1] < inputS && AddSample[l_s1] >= 0) { - s_1 = input[l_s1 + AddSample[l_s1 - minLine] * (short)inputL]; + s_1 = input[l_s1 + AddSample[l_s1] * (short)inputL]; sign += s_1; - for (short l_s2 = l_s1 + 1; l_s2 < maxLine; ++l_s2) + for (short l_s2 = l_s1 + 1; l_s2 < inputL; ++l_s2) { - if (AddSample[l_s2 - minLine] < inputS && AddSample[l_s2 - minLine] >= 0) + if (AddSample[l_s2] < inputS && AddSample[l_s2] >= 0) { - s_2 = input[l_s2 + AddSample[l_s2 - minLine] * (short)inputL]; + s_2 = input[l_s2 + AddSample[l_s2] * (short)inputL]; - mult = s_2 * apodisation[(int)((l_s2 - minLine)*apod_mult)] * s_1 * apodisation[(int)((l_s1 - minLine)*apod_mult)]; + mult = s_2 * apodisation[(int)((l_s2)*apod_mult)] * s_1 * apodisation[(int)((l_s1)*apod_mult)]; output[sample*(short)outputL + line] += sqrt(fabs(mult)) * ((mult > 0) - (mult < 0)); } } } else --usedLines; } output[sample*(short)outputL + line] = output[sample*(short)outputL + line] / (float)(pow(usedLines, 2) - (usedLines - 1)) * ((sign > 0) - (sign < 0)); delete[] AddSample; } } diff --git a/Modules/PhotoacousticsAlgorithms/source/utils/mitkBeamformingUtils.cpp.autosave b/Modules/PhotoacousticsAlgorithms/source/utils/mitkBeamformingUtils.cpp.autosave new file mode 100644 index 0000000000..7194bd8458 --- /dev/null +++ b/Modules/PhotoacousticsAlgorithms/source/utils/mitkBeamformingUtils.cpp.autosave @@ -0,0 +1,276 @@ +/*=================================================================== +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. + +===================================================================*/ + +#include "mitkProperties.h" +#include "mitkImageReadAccessor.h" +#include +#include +#include +#include +#include +#include "mitkImageCast.h" +#include "mitkBeamformingUtils.h" + +mitk::BeamformingUtils::BeamformingUtils() +{ +} + +mitk::BeamformingUtils::~BeamformingUtils() +{ +} + +float* mitk::BeamformingUtils::VonHannFunction(int samples) +{ + float* ApodWindow = new float[samples]; + + for (int n = 0; n < samples; ++n) + { + ApodWindow[n] = (1 - cos(2 * itk::Math::pi * n / (samples - 1))) / 2; + } + + return ApodWindow; +} + +float* mitk::BeamformingUtils::HammFunction(int samples) +{ + float* ApodWindow = new float[samples]; + + for (int n = 0; n < samples; ++n) + { + ApodWindow[n] = 0.54 - 0.46*cos(2 * itk::Math::pi*n / (samples - 1)); + } + + return ApodWindow; +} + +float* mitk::BeamformingUtils::BoxFunction(int samples) +{ + float* ApodWindow = new float[samples]; + + for (int n = 0; n < samples; ++n) + { + ApodWindow[n] = 1; + } + + return ApodWindow; +} + +void mitk::BeamformingUtils::DASSphericalLine( + float* input, float* output, float inputDim[2], float outputDim[2], + const short& line, const mitk::BeamformingSettings::Pointer config) +{ + const float* apodisation = config->GetApodizationFunction(); + const short apodArraySize = config->GetApodizationArraySize(); + MITK_INFO << "=============== apodisation = " << apodisation; + MITK_INFO << "=============== apodArraySize = " << apodArraySize; + const float* elementHeights = config->GetElementHeights(); + const float* elementPositions = config->GetElementPositions(); + + float& inputS = inputDim[1]; + float& inputL = inputDim[0]; + + float& outputS = outputDim[1]; + float& outputL = outputDim[0]; + + short AddSample = 0; + float l_p = 0; + float s_i = 0; + + float apod_mult = 1; + + float totalSamples_i = (float)(config->GetReconstructionDepth()) / (float)(config->GetSpeedOfSound() * config->GetTimeSpacing()); + totalSamples_i = totalSamples_i <= inputS ? totalSamples_i : inputS; + + l_p = (float)line / outputL * config->GetHorizontalExtent(); + + for (short sample = 0; sample < outputS; ++sample) + { + s_i = (float)sample / outputS * totalSamples_i; + usedLines = inputL; + apod_mult = (float)apodArraySize / (float)usedLines; + + for (short l_s = 0; l_s < inputL; ++l_s) + { + AddSample = (int)sqrt( + pow(s_i-elementHeights[l_s]/(config->GetSpeedOfSound()*config->GetTimeSpacing()), 2) + + + pow((1 / (config->GetTimeSpacing()*config->GetSpeedOfSound())) * (l_p - elementPositions[l_s]), 2) + ) + (1 - config->GetIsPhotoacousticImage())*s_i; + if (AddSample < inputS && AddSample >= 0) + output[sample*(short)outputL + line] += input[l_s + AddSample*(short)inputL] * + apodisation[(short)((l_s)*apod_mult)]; + else + --usedLines; + } + output[sample*(short)outputL + line] = output[sample*(short)outputL + line] / usedLines; + } +} + +void mitk::BeamformingUtils::DMASSphericalLine( + float* input, float* output, float inputDim[2], float outputDim[2], + const short& line, const mitk::BeamformingSettings::Pointer config) +{ + const float* apodisation = config->GetApodizationFunction(); + const short apodArraySize = config->GetApodizationArraySize(); + + const float* elementHeights = config->GetElementHeights(); + const float* elementPositions = config->GetElementPositions(); + + float& inputS = inputDim[1]; + float& inputL = inputDim[0]; + + float& outputS = outputDim[1]; + float& outputL = outputDim[0]; + + float l_p = 0; + float s_i = 0; + + float apod_mult = 1; + float mult = 0; + + short usedLines = inputL; + + float totalSamples_i = (float)(config->GetReconstructionDepth()) / + (float)(config->GetSpeedOfSound() * config->GetTimeSpacing()); + totalSamples_i = totalSamples_i <= inputS ? totalSamples_i : inputS; + + l_p = (float)line / outputL * config->GetHorizontalExtent(); + + for (short sample = 0; sample < outputS; ++sample) + { + s_i = (float)sample / outputS * totalSamples_i; + usedLines = inputL; + apod_mult = (float)apodArraySize / (float)usedLines; + + //calculate the AddSamples beforehand to save some time + short* AddSample = new short[usedLines]; + for (short l_s = 0; l_s < usedLines; ++l_s) + { + AddSample[l_s] = (int)sqrt( + pow(s_i - elementHeights[l_s] / (config->GetSpeedOfSound()*config->GetTimeSpacing()), 2) + + + pow((1 / (config->GetTimeSpacing()*config->GetSpeedOfSound())) * (l_p - elementPositions[l_s]), 2) + ) + (1 - config->GetIsPhotoacousticImage())*s_i; + } + + float s_1 = 0; + float s_2 = 0; + + for (short l_s1 = 0; l_s1 < usedLines - 1; ++l_s1) + { + if (AddSample[l_s1] < inputS && AddSample[l_s1] >= 0) + { + for (short l_s2 = l_s1 + 1; l_s2 < inputL; ++l_s2) + { + if (AddSample[l_s2] < inputS && AddSample[l_s2] >= 0) + { + s_2 = input[l_s2 + AddSample[l_s2] * (short)inputL]; + s_1 = input[l_s1 + AddSample[l_s1] * (short)inputL]; + + mult = s_2 * apodisation[(int)((l_s2)*apod_mult)] * s_1 * apodisation[(int)((l_s1)*apod_mult)]; + output[sample*(short)outputL + line] += sqrt(fabs(mult)) * ((mult > 0) - (mult < 0)); + } + } + } + else + --usedLines; + } + + output[sample*(short)outputL + line] = output[sample*(short)outputL + line] / (float)(pow(usedLines, 2) - (usedLines - 1)); + + delete[] AddSample; + } +} + +void mitk::BeamformingUtils::sDMASSphericalLine( + float* input, float* output, float inputDim[2], float outputDim[2], + const short& line, const mitk::BeamformingSettings::Pointer config) +{ + const float* apodisation = config->GetApodizationFunction(); + const short apodArraySize = config->GetApodizationArraySize(); + + const float* elementHeights = config->GetElementHeights(); + const float* elementPositions = config->GetElementPositions(); + + float& inputS = inputDim[1]; + float& inputL = inputDim[0]; + + float& outputS = outputDim[1]; + float& outputL = outputDim[0]; + + float l_p = 0; + float s_i = 0; + + float apod_mult = 1; + + float mult = 0; + + short usedLines = inputL; + + float totalSamples_i = (float)(config->GetReconstructionDepth()) / + (float)(config->GetSpeedOfSound() * config->GetTimeSpacing()); + totalSamples_i = totalSamples_i <= inputS ? totalSamples_i : inputS; + + l_p = (float)line / outputL * config->GetHorizontalExtent(); + + for (short sample = 0; sample < outputS; ++sample) + { + s_i = (float)sample / outputS * totalSamples_i; + usedLines = inputL; + apod_mult = (float)apodArraySize / (float)usedLines; + + //calculate the AddSamples beforehand to save some time + short* AddSample = new short[usedLines]; + for (short l_s = 0; l_s < usedLines; ++l_s) + { + AddSample[l_s] = (int)sqrt( + pow(s_i - elementHeights[l_s] / (config->GetSpeedOfSound()*config->GetTimeSpacing()), 2) + + + pow((1 / (config->GetTimeSpacing()*config->GetSpeedOfSound())) * (l_p - elementPositions[l_s]), 2) + ) + (1 - config->GetIsPhotoacousticImage())*s_i; + } + + float s_1 = 0; + float s_2 = 0; + float sign = 0; + + for (short l_s1 = 0; l_s1 < inputL - 1; ++l_s1) + { + if (AddSample[l_s1] < inputS && AddSample[l_s1] >= 0) + { + s_1 = input[l_s1 + AddSample[l_s1] * (short)inputL]; + sign += s_1; + + for (short l_s2 = l_s1 + 1; l_s2 < inputL; ++l_s2) + { + if (AddSample[l_s2] < inputS && AddSample[l_s2] >= 0) + { + s_2 = input[l_s2 + AddSample[l_s2] * (short)inputL]; + + mult = s_2 * apodisation[(int)((l_s2)*apod_mult)] * s_1 * apodisation[(int)((l_s1)*apod_mult)]; + output[sample*(short)outputL + line] += sqrt(fabs(mult)) * ((mult > 0) - (mult < 0)); + } + } + } + else + --usedLines; + } + + output[sample*(short)outputL + line] = output[sample*(short)outputL + line] / (float)(pow(usedLines, 2) - (usedLines - 1)) * ((sign > 0) - (sign < 0)); + + delete[] AddSample; + } +} 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 28a016d489..d0d2edb679 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,1317 +1,1317 @@ PAImageProcessingControls 0 0 601 890 0 0 QmitkTemplate 0 0 - 0 + 1 0 0 583 732 Bandpass 10 0 301 529 QLayout::SetDefaultConstraint 0 0 0 3 0.010000000000000 200.000000000000000 0.100000000000000 8.000000000000000 [MHz] f High Pass [MHz] f Low Pass 0 0 3 200.000000000000000 0.100000000000000 0.100000000000000 1 200.000000000000000 3000.000000000000000 5.000000000000000 1500.000000000000000 [m/s] Speed of Sound 1.000000000000000 0.100000000000000 0.500000000000000 <html><head/><body><p><span style=" font-family:'sans-serif'; font-size:14px; color:#222222; background-color:#ffffff;">At </span><span style=" font-family:'sans-serif'; font-size:14px; font-style:italic; color:#222222; background-color:#ffffff;">α</span><span style=" font-family:'sans-serif'; font-size:14px; color:#222222; background-color:#ffffff;"> = 0 it's rectangular, at </span><span style=" font-family:'sans-serif'; font-size:14px; font-style:italic; color:#222222; background-color:#ffffff;">α</span><span style=" font-family:'sans-serif'; font-size:14px; color:#222222; background-color:#ffffff;"> = 1 it's a Hann window. Both flanks can be set independently.</span></p></body></html> Tukey Window α High Pass 2 1.000000000000000 0.100000000000000 0.500000000000000 <html><head/><body><p><span style=" font-family:'sans-serif'; font-size:14px; color:#222222; background-color:#ffffff;">At </span><span style=" font-family:'sans-serif'; font-size:14px; font-style:italic; color:#222222; background-color:#ffffff;">α</span><span style=" font-family:'sans-serif'; font-size:14px; color:#222222; background-color:#ffffff;"> = 0 it's rectangular, at </span><span style=" font-family:'sans-serif'; font-size:14px; font-style:italic; color:#222222; background-color:#ffffff;">α</span><span style=" font-family:'sans-serif'; font-size:14px; color:#222222; background-color:#ffffff;"> = 1 it's a Hann window. Both flanks can be set independently.</span></p></body></html> Tukey Window α Low Pass 0 25 16777215 25 <html><head/><body><p><span style=" font-size:10pt; font-weight:600; color:#ff0000;">Select input data in datamanager!</span></p></body></html> Apply Bandpass <html><head/><body><p>If not checked, treat input as raw US/PA data with Y-axis as a time coordinate in microseconds</p></body></html> Assume Spatial Coordinates Qt::Vertical 20 40 0 0 583 732 Beamforming 10 0 - 301 + 306 661 5 0 0 PA Image US Image <html><head/><body><p>which type of image do you plan to reconstruct?</p></body></html> Raw Data Type 0 0 DAS DMAS sDMAS 0 0 Von Hann Hamming Box [mm] Transducer Pitch Apodization Function Beamforming Algorithm 0 0 3 0.010000000000000 9.000000000000000 0.050000000000000 0.340000000000000 Probe Geomentry 0 0 1 1.000000000000000 180.000000000000000 45.000000000000000 <html><head/><body><p>... of the transducer elements.</p></body></html> [°] Sensitive Angle [mm] Reconstruction Depth 4 300.000000000000000 0.100000000000000 60.000000000000000 Linear Linear Concave 2 - 0.010000000000000 + 0.000000000000000 399.990000000000009 40.000000000000000 [mm] Concave Probe Radius <html><head/><body><p>Some setups' hardware produces signal delays that need to be cropped out of the image before performing beamforming. To do this, select this box.</p></body></html> Consider Hardware Delay [µs] false 0.100000000000000 1.000000000000000 <html><head/><body><p>... from y-spacing in the selected raw data. If this is switched of &quot;Manual Scan Depth&quot; is used.</p></body></html> Automatic Get Depth true [mm] Manual Scan Depth <html><head/><body><p>Using GPU is recommended - It is so much faster.</p></body></html> Compute On GPU true 0 0 64 2048 128 512 0 0 2 2048 64 256 Reconstructed Lines 75 true Advanced Options 0 0 4 300.000000000000000 0.100000000000000 60.000000000000000 0 0 1 200.000000000000000 3000.000000000000000 5.000000000000000 1500.000000000000000 <html><head/><body><p>... Good default - change only if you know what you are doing.</p></body></html> Samples [m/s] Speed of Sound 0 0 256 16384 256 2048 <html><head/><body><p>Set automatically from selected raw data x-Geometry.</p></body></html> Transducer Elements 75 true 0 25 16777215 25 <html><head/><body><p><span style=" font-size:10pt; font-weight:600; color:#ff0000;">Select input data in datamanager!</span></p></body></html> true 0 0 100 0 Qt::Vertical 20 40 Apply Beamforming 0 0 583 732 Cropping 10 0 330 301 999999999 Left Bottom minimal beamformed slice First Slice 999999999 Right Select Slices Apply Crop Filer 0 25 16777215 25 <html><head/><body><p><span style=" font-size:10pt; font-weight:600; color:#ff0000;">Select input data in datamanager!</span></p></body></html> Top false 99999 Cut N Pixels from ... 999999999 1 0 999999999 0 false 99999 10 Maximal beamformed slice Last Slice Qt::Vertical 20 40 0 0 583 732 B-mode Generation 10 0 301 261 0 25 16777215 25 <html><head/><body><p><span style=" font-size:10pt; font-weight:600; color:#ff0000;">Select input data in datamanager!</span></p></body></html> Envelope Detection Envelope Detection Absolute Filter <html><head/><body><p>The image will still be oversampled after B-Mode generation. This default performs a sensible downsampling.</p></body></html> Resample Image after B-mode true 0 0 13 0 11 3 0.010000000000000 1.000000000000000 0.010000000000000 0.050000000000000 [mm] Resampled y-Spacing <html><head/><body><p>Adds a log() filter after the B-mode filter. Often used in US B-mode images.</p></body></html> Logarithmic Compression Use GPU 0 0 Do image processing Apply B-mode Filter Qt::Vertical 20 40 0 0 583 732 Legacy Batch Processing 10 0 309 231 300 100 <html><head/><body><p><span style=" font-weight:600;">Note:</span> This batch processing tool is depricated <br/>and is no longer tested. It is recommended to <br/>use the PA command line tool to process large <br/>amounts of files with consistent settings.</p></body></html> Bandpass true Crop true Save false Save false Save false Beamform true B-Mode true Save true Start Batch Processing Qt::Vertical 20 40