diff --git a/Modules/PhotoacousticsAlgorithms/Resources/DAS.cl b/Modules/PhotoacousticsAlgorithms/Resources/DAS.cl index 42d954137b..34c70224f0 100644 --- a/Modules/PhotoacousticsAlgorithms/Resources/DAS.cl +++ b/Modules/PhotoacousticsAlgorithms/Resources/DAS.cl @@ -1,116 +1,124 @@ /*=================================================================== 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 = 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]; float apod_mult = (float)apodArraySize / (float)curUsedLines; unsigned short Delay = 0; float output = 0; float mult = 0; - for (short l_s = 0; l_s < inputL; ++l_s) + for (short l_s = minLine; l_s < maxLine; ++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)*apod_mult)] * dSource[(int)(globalPosZ * inputL * inputS + Delay * inputL + l_s)]; + output += apodArray[(int)((l_s - minLine)*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 + char isPAImage, + __global unsigned short* usedLines // parameters ) { // 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 = 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]; apod_mult = (float)apodArraySize / curUsedLines; - for (int l_s = 0; l_s < inputL; ++l_s) + for (int l_s = minLine; l_s < maxLine; ++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)*apod_mult)]; + output += dSource[(int)(globalPosZ * inputL * inputS + l_s + AddSample*inputL)] * apodArray[(int)((l_s - minLine)*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 cb7ac05f4e..b71208e6e8 100644 --- a/Modules/PhotoacousticsAlgorithms/Resources/DelayCalculation.cl +++ b/Modules/PhotoacousticsAlgorithms/Resources/DelayCalculation.cl @@ -1,43 +1,44 @@ /*=================================================================== 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 new file mode 100644 index 0000000000..132c7d88fd --- /dev/null +++ b/Modules/PhotoacousticsAlgorithms/Resources/UsedLinesCalculation.cl @@ -0,0 +1,170 @@ +/*=================================================================== + +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 7b900861de..87f8116092 100644 --- a/Modules/PhotoacousticsAlgorithms/files.cmake +++ b/Modules/PhotoacousticsAlgorithms/files.cmake @@ -1,29 +1,31 @@ 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 00799f0ac6..bd9bf91913 100644 --- a/Modules/PhotoacousticsAlgorithms/include/OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.h +++ b/Modules/PhotoacousticsAlgorithms/include/OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.h @@ -1,145 +1,147 @@ /*=================================================================== 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 69960e3356..3e60b4ec5c 100644 --- a/Modules/PhotoacousticsAlgorithms/include/OpenCLFilter/mitkPhotoacousticOCLDelayCalculation.h +++ b/Modules/PhotoacousticsAlgorithms/include/OpenCLFilter/mitkPhotoacousticOCLDelayCalculation.h @@ -1,78 +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 _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/mitkPhotoacousticOCLDelayCalculation.h b/Modules/PhotoacousticsAlgorithms/include/OpenCLFilter/mitkPhotoacousticOCLUsedLinesCalculation.h similarity index 68% copy from Modules/PhotoacousticsAlgorithms/include/OpenCLFilter/mitkPhotoacousticOCLDelayCalculation.h copy to Modules/PhotoacousticsAlgorithms/include/OpenCLFilter/mitkPhotoacousticOCLUsedLinesCalculation.h index 69960e3356..52021ea7c5 100644 --- a/Modules/PhotoacousticsAlgorithms/include/OpenCLFilter/mitkPhotoacousticOCLDelayCalculation.h +++ b/Modules/PhotoacousticsAlgorithms/include/OpenCLFilter/mitkPhotoacousticOCLUsedLinesCalculation.h @@ -1,78 +1,84 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ -#ifndef _MITKPHOTOACOUSTICSDELAYCALC_H_ -#define _MITKPHOTOACOUSTICSDELAYCALC_H_ +#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 the delays used for beamforming. + * \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 OCLDelayCalculation : public OclDataSetToDataSetFilter, public itk::Object + class OCLUsedLinesCalculation : public OclDataSetToDataSetFilter, public itk::Object { public: - mitkClassMacroItkParent(OCLDelayCalculation, itk::Object); + mitkClassMacroItkParent(OCLUsedLinesCalculation, itk::Object); mitkNewMacro1Param(Self, mitk::BeamformingSettings::Pointer); void Update(); + void SetElementHeightsBuffer(cl_mem elementHeightsBuffer); + void SetElementPositionsBuffer(cl_mem elementPositionsBuffer); + protected: - OCLDelayCalculation(mitk::BeamformingSettings::Pointer settings); - virtual ~OCLDelayCalculation(); + /** 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; - unsigned int m_BufferSize; - float m_DelayMultiplicatorRaw; - char m_IsPAImage; + 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 8724808d36..6ee498e470 100644 --- a/Modules/PhotoacousticsAlgorithms/include/mitkBeamformingFilter.h +++ b/Modules/PhotoacousticsAlgorithms/include/mitkBeamformingFilter.h @@ -1,89 +1,90 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #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 891dbf0d88..8684d9632e 100644 --- a/Modules/PhotoacousticsAlgorithms/include/mitkBeamformingSettings.h +++ b/Modules/PhotoacousticsAlgorithms/include/mitkBeamformingSettings.h @@ -1,242 +1,248 @@ /*=================================================================== 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 6fe1eecc09..fff7126538 100644 --- a/Modules/PhotoacousticsAlgorithms/include/mitkBeamformingUtils.h +++ b/Modules/PhotoacousticsAlgorithms/include/mitkBeamformingUtils.h @@ -1,68 +1,72 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #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 ba58a646ef..dc3ad02f83 100644 --- a/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.cpp +++ b/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.cpp @@ -1,314 +1,325 @@ /*=================================================================== 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_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)); + 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)); } 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 ec8eb77fb4..b6f4e2f0f7 100644 --- a/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLDelayCalculation.cpp +++ b/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLDelayCalculation.cpp @@ -1,127 +1,128 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #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_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)); + 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)); 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 new file mode 100644 index 0000000000..8da3c4c197 --- /dev/null +++ b/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLUsedLinesCalculation.cpp @@ -0,0 +1,167 @@ +/*=================================================================== + +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 223778796c..894802fe76 100644 --- a/Modules/PhotoacousticsAlgorithms/source/filters/mitkBeamformingSettings.cpp +++ b/Modules/PhotoacousticsAlgorithms/source/filters/mitkBeamformingSettings.cpp @@ -1,131 +1,141 @@ /*=================================================================== 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_ProbeRadius(probeRadius), + m_MinMaxLines(nullptr) { 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 854473ea8f..b367bb57f3 100644 --- a/Modules/PhotoacousticsAlgorithms/source/utils/mitkBeamformingUtils.cpp +++ b/Modules/PhotoacousticsAlgorithms/source/utils/mitkBeamformingUtils.cpp @@ -1,277 +1,420 @@ /*=================================================================== 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; - usedLines = inputL; + + minLine = config->GetMinMaxLines()[2*sample*(short)outputL + 2*line]; + maxLine = config->GetMinMaxLines()[2*sample*(short)outputL + 2*line + 1]; + usedLines = (maxLine - minLine); + apod_mult = (float)apodArraySize / (float)usedLines; - for (short l_s = 0; l_s < inputL; ++l_s) + for (short l_s = minLine; l_s < maxLine; ++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)]; + apodisation[(short)((l_s - minLine)*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 = inputL; + 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; - usedLines = inputL; + + minLine = config->GetMinMaxLines()[2 * sample*(short)outputL + 2 * line]; + maxLine = config->GetMinMaxLines()[2 * sample*(short)outputL + 2 * line + 1]; + usedLines = (maxLine - minLine); + 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) + short* AddSample = new short[maxLine - minLine]; + for (short l_s = 0; l_s < maxLine - minLine; ++l_s) { AddSample[l_s] = (int)sqrt( - pow(s_i - elementHeights[l_s] / (config->GetSpeedOfSound()*config->GetTimeSpacing()), 2) + pow(s_i - elementHeights[l_s + minLine] / (config->GetSpeedOfSound()*config->GetTimeSpacing()), 2) + - pow((1 / (config->GetTimeSpacing()*config->GetSpeedOfSound())) * (l_p - elementPositions[l_s]), 2) + pow((1 / (config->GetTimeSpacing()*config->GetSpeedOfSound())) * (l_p - elementPositions[l_s + minLine]), 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) + for (short l_s1 = minLine; l_s1 < maxLine - 1; ++l_s1) { - if (AddSample[l_s1] < inputS && AddSample[l_s1] >= 0) + if (AddSample[l_s1 - minLine] < inputS && AddSample[l_s1 - minLine] >= 0) { - for (short l_s2 = l_s1 + 1; l_s2 < inputL; ++l_s2) + for (short l_s2 = l_s1 + 1; l_s2 < maxLine; ++l_s2) { - if (AddSample[l_s2] < inputS && AddSample[l_s2] >= 0) + if (AddSample[l_s2 - minLine] < inputS && AddSample[l_s2 - minLine] >= 0) { - s_2 = input[l_s2 + AddSample[l_s2] * (short)inputL]; - s_1 = input[l_s1 + AddSample[l_s1] * (short)inputL]; + s_2 = input[l_s2 + AddSample[l_s2 - minLine] * (short)inputL]; + s_1 = input[l_s1 + AddSample[l_s1 - minLine] * (short)inputL]; - mult = s_2 * apodisation[(int)((l_s2)*apod_mult)] * s_1 * apodisation[(int)((l_s1)*apod_mult)]; + mult = s_2 * apodisation[(int)((l_s2 - minLine)*apod_mult)] * s_1 * apodisation[(int)((l_s1 - minLine)*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 = inputL; + 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; - usedLines = inputL; + + minLine = config->GetMinMaxLines()[2 * sample*(short)outputL + 2 * line]; + maxLine = config->GetMinMaxLines()[2 * sample*(short)outputL + 2 * line + 1]; + usedLines = (maxLine - minLine); + 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) + short* AddSample = new short[maxLine - minLine]; + for (short l_s = 0; l_s < maxLine - minLine; ++l_s) { AddSample[l_s] = (int)sqrt( - pow(s_i - elementHeights[l_s] / (config->GetSpeedOfSound()*config->GetTimeSpacing()), 2) + pow(s_i - elementHeights[l_s + minLine] / (config->GetSpeedOfSound()*config->GetTimeSpacing()), 2) + - pow((1 / (config->GetTimeSpacing()*config->GetSpeedOfSound())) * (l_p - elementPositions[l_s]), 2) + pow((1 / (config->GetTimeSpacing()*config->GetSpeedOfSound())) * (l_p - elementPositions[l_s + minLine]), 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) + for (short l_s1 = minLine; l_s1 < maxLine - 1; ++l_s1) { - if (AddSample[l_s1] < inputS && AddSample[l_s1] >= 0) + if (AddSample[l_s1 - minLine] < inputS && AddSample[l_s1 - minLine] >= 0) { - s_1 = input[l_s1 + AddSample[l_s1] * (short)inputL]; + s_1 = input[l_s1 + AddSample[l_s1 - minLine] * (short)inputL]; sign += s_1; - for (short l_s2 = l_s1 + 1; l_s2 < inputL; ++l_s2) + for (short l_s2 = l_s1 + 1; l_s2 < maxLine; ++l_s2) { - if (AddSample[l_s2] < inputS && AddSample[l_s2] >= 0) + if (AddSample[l_s2 - minLine] < inputS && AddSample[l_s2 - minLine] >= 0) { - s_2 = input[l_s2 + AddSample[l_s2] * (short)inputL]; + s_2 = input[l_s2 + AddSample[l_s2 - minLine] * (short)inputL]; - mult = s_2 * apodisation[(int)((l_s2)*apod_mult)] * s_1 * apodisation[(int)((l_s1)*apod_mult)]; + mult = s_2 * apodisation[(int)((l_s2 - minLine)*apod_mult)] * s_1 * apodisation[(int)((l_s1 - minLine)*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 deleted file mode 100644 index 7194bd8458..0000000000 --- a/Modules/PhotoacousticsAlgorithms/source/utils/mitkBeamformingUtils.cpp.autosave +++ /dev/null @@ -1,276 +0,0 @@ -/*=================================================================== -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 d0d2edb679..28a016d489 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 - 1 + 0 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 - 306 + 301 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.000000000000000 + 0.010000000000000 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