diff --git a/Modules/PhotoacousticsAlgorithms/include/OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.h b/Modules/PhotoacousticsAlgorithms/include/OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.h index ff7d62b1be..508c4d05ca 100644 --- a/Modules/PhotoacousticsAlgorithms/include/OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.h +++ b/Modules/PhotoacousticsAlgorithms/include/OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.h @@ -1,145 +1,146 @@ /*=================================================================== 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 "mitkPhotoacousticOCLDelayCalculation.h" #include "mitkPhotoacousticOCLUsedLinesCalculation.h" #include "mitkBeamformingSettings.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); + void SetInput(Image::Pointer image, unsigned int inputSlices); /** * 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_MemoryLocationsBuffer; cl_mem m_DelaysBuffer; cl_mem m_UsedLinesBuffer; }; } #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/mitkBeamformingSettings.h b/Modules/PhotoacousticsAlgorithms/include/mitkBeamformingSettings.h index 718652f9bb..dd837b4f10 100644 --- a/Modules/PhotoacousticsAlgorithms/include/mitkBeamformingSettings.h +++ b/Modules/PhotoacousticsAlgorithms/include/mitkBeamformingSettings.h @@ -1,254 +1,214 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef 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); mitkNewMacro1Param(BeamformingSettings, std::string); itkCloneMacro(Self); /** \brief Available delay calculation methods: * - Spherical delay for best results. * - A quadratic Taylor approximation for slightly faster results with hardly any quality loss. */ enum DelayCalc { QuadApprox, Spherical }; /** \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 }; 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(UpperCutoff, unsigned int); - itkGetConstMacro(Partial, bool); - itkGetConstMacro(CropBounds, const unsigned int*); itkGetConstMacro(InputDim, const unsigned int*); itkGetConstMacro(UseGPU, bool); + itkGetConstMacro(GPUBatchSize, unsigned int); itkGetConstMacro(DelayCalculationMethod, DelayCalc); itkGetConstMacro(ApodizationFunction, const float*); itkGetConstMacro(Apod, Apodization); itkGetConstMacro(ApodizationArraySize, int); itkGetConstMacro(Algorithm, BeamformingAlgorithm); - itkGetConstMacro(UseBP, bool); - itkGetConstMacro(BPHighPass, float); - itkGetConstMacro(BPLowPass, 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 !((abs(lhs->GetAngle() - rhs->GetAngle()) < 0.01f) && // 0.01 degree error margin (lhs->GetApod() == rhs->GetApod()) && (lhs->GetDelayCalculationMethod() == rhs->GetDelayCalculationMethod()) && (lhs->GetIsPhotoacousticImage() == rhs->GetIsPhotoacousticImage()) && (abs(lhs->GetPitchInMeters() - rhs->GetPitchInMeters()) < 0.000001f) && // 0.0001 mm error margin (lhs->GetReconstructionLines() == rhs->GetReconstructionLines()) && (lhs->GetSamplesPerLine() == rhs->GetSamplesPerLine()) && (abs(lhs->GetSpeedOfSound() - rhs->GetSpeedOfSound()) < 0.01f) && (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 upperCutoff, - bool partial, - unsigned int* cropBounds, unsigned int* inputDim, bool useGPU, + unsigned int GPUBatchSize, DelayCalc delayCalculationMethod, Apodization apod, unsigned int apodizationArraySize, - BeamformingAlgorithm algorithm, - bool useBP, - float BPHighPass, - float BPLowPass) + BeamformingAlgorithm algorithm) { Pointer smartPtr = new BeamformingSettings(pitchInMeters, speedOfSound, timeSpacing, angle, isPhotoacousticImage, samplesPerLine, reconstructionLines, - upperCutoff, - partial, - cropBounds, inputDim, useGPU, + GPUBatchSize, delayCalculationMethod, apod, apodizationArraySize, - algorithm, - useBP, - BPHighPass, - BPLowPass); + algorithm); smartPtr->UnRegister(); return smartPtr; } protected: /** */ BeamformingSettings(std::string xmlFile); /** */ BeamformingSettings(float pitchInMeters, float speedOfSound, float timeSpacing, float angle, bool isPhotoacousticImage, unsigned int samplesPerLine, unsigned int reconstructionLines, - unsigned int upperCutoff, - bool partial, - unsigned int* cropBounds, unsigned int* inputDim, bool useGPU, + unsigned int GPUBatchSize, DelayCalc delayCalculationMethod, Apodization apod, unsigned int apodizationArraySize, - BeamformingAlgorithm algorithm, - bool useBP, - float BPHighPass, - float BPLowPass + BeamformingAlgorithm algorithm ); ~BeamformingSettings(); /** \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 how many voxels should be cut off from the top of the image before beamforming, to potentially avoid artifacts. - */ - unsigned int m_UpperCutoff; - - /** \brief Sets whether only the slices selected by mitk::BeamformingSettings::CropBounds should be beamformed. - */ - bool m_Partial; - - /** \brief Sets the first and last slice to be beamformed. - */ - const unsigned int* m_CropBounds; - /** \brief Sets the dimensions of the inputImage. */ const unsigned int* m_InputDim; /** \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 + */ + /** \brief Sets how the delays for beamforming should be calculated. */ DelayCalc m_DelayCalculationMethod; 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 whether after beamforming a bandpass should be automatically applied - */ - bool m_UseBP; - - /** \brief Sets the position at which lower frequencies are completely cut off in Hz. - */ - float m_BPHighPass; - - /** \brief Sets the position at which higher frequencies are completely cut off in Hz. - */ - float m_BPLowPass; }; } #endif //MITK_BEAMFORMING_SETTINGS diff --git a/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.cpp b/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.cpp index 16960e2e47..4d02104ace 100644 --- a/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.cpp +++ b/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.cpp @@ -1,260 +1,259 @@ /*=================================================================== 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_InputImage(mitk::Image::New()), m_ApodizationBuffer(nullptr), m_MemoryLocationsBuffer(nullptr), m_DelaysBuffer(nullptr), m_UsedLinesBuffer(nullptr), + m_inputSlices(1), m_Conf(settings) { 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); } 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_ulong globalMemSize = oclGetGlobalMemSize(resources->GetCurrentDevice());*/ //Initialize the Output try { size_t outputSize = (size_t)m_Conf->GetReconstructionLines() * (size_t)m_Conf->GetSamplesPerLine() * - (size_t)m_Conf->GetInputDim()[2]; + (size_t)m_inputSlices; m_OutputDim[0] = m_Conf->GetReconstructionLines(); m_OutputDim[1] = m_Conf->GetSamplesPerLine(); - m_OutputDim[2] = m_Conf->GetInputDim()[2]; + 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; } //TODO FIXME 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; } us::ServiceReference ref = GetModuleContext()->GetServiceReference(); OclResourceService* resources = GetModuleContext()->GetService(ref); cl_context gpuContext = resources->GetContext(); if (m_ApodizationBuffer) clReleaseMemObject(m_ApodizationBuffer); this->m_ApodizationBuffer = clCreateBuffer(gpuContext, CL_MEM_READ_ONLY | CL_MEM_USE_HOST_PTR, sizeof(float) * m_ApodArraySize, const_cast(m_Apodisation), &clErr); CHECK_OCL_ERR(clErr); // calculate used lines 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(); //m_ConfOld = m_Conf; } void mitk::PhotoacousticOCLBeamformingFilter::Execute() { cl_int clErr = 0; UpdateDataBuffers(); unsigned int reconstructionLines = this->m_Conf->GetReconstructionLines(); unsigned int samplesPerLine = this->m_Conf->GetSamplesPerLine(); clErr = clSetKernelArg(this->m_PixelCalculation, 2, sizeof(cl_mem), &(this->m_UsedLinesBuffer)); clErr |= clSetKernelArg(this->m_PixelCalculation, 3, sizeof(cl_mem), &(this->m_DelaysBuffer)); clErr |= clSetKernelArg(this->m_PixelCalculation, 4, sizeof(cl_mem), &(this->m_ApodizationBuffer)); clErr |= clSetKernelArg(this->m_PixelCalculation, 5, sizeof(cl_ushort), &(this->m_ApodArraySize)); clErr |= clSetKernelArg(this->m_PixelCalculation, 6, sizeof(cl_uint), &(this->m_Conf->GetInputDim()[0])); clErr |= clSetKernelArg(this->m_PixelCalculation, 7, sizeof(cl_uint), &(this->m_Conf->GetInputDim()[1])); - clErr |= clSetKernelArg(this->m_PixelCalculation, 8, sizeof(cl_uint), &(this->m_Conf->GetInputDim()[2])); + clErr |= clSetKernelArg(this->m_PixelCalculation, 8, sizeof(cl_uint), &(m_OutputDim[2])); clErr |= clSetKernelArg(this->m_PixelCalculation, 9, sizeof(cl_uint), &(reconstructionLines)); clErr |= clSetKernelArg(this->m_PixelCalculation, 10, sizeof(cl_uint), &(samplesPerLine)); - // execute the filter on a 3D NDRange + // 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, 16, 50)) mitkThrow() << "openCL Error when executing Kernel"; } else { if (!this->ExecuteKernelChunksInBatches(m_PixelCalculation, 3, m_ChunkSize, 16, 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() { - MITK_INFO << "Initializing OCL Beamforming Filter..."; bool buildErr = true; cl_int clErr = 0; if (OclFilter::Initialize()) { switch (m_Conf->GetAlgorithm()) { case BeamformingSettings::BeamformingAlgorithm::DAS: { this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "ckDAS", &clErr); break; } case BeamformingSettings::BeamformingAlgorithm::DMAS: { this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "ckDMAS", &clErr); break; } case BeamformingSettings::BeamformingAlgorithm::sDMAS: { 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; } } buildErr |= CHECK_OCL_ERR(clErr); } CHECK_OCL_ERR(clErr); - MITK_INFO << "Initializing OCL Beamforming Filter...[Done]"; - return (OclFilter::IsInitialized() && buildErr); } -void mitk::PhotoacousticOCLBeamformingFilter::SetInput(mitk::Image::Pointer image) +void mitk::PhotoacousticOCLBeamformingFilter::SetInput(mitk::Image::Pointer image, unsigned int inputSlices) { OclDataSetToDataSetFilter::SetInput(image); m_InputImage = image; + m_inputSlices = inputSlices; } 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); delete[] pData; } MITK_DEBUG << "Image Initialized."; return outputImage; } void* mitk::PhotoacousticOCLBeamformingFilter::GetOutput() { return OclDataSetToDataSetFilter::GetOutput(); } #endif diff --git a/Modules/PhotoacousticsAlgorithms/source/filters/mitkBeamformingFilter.cpp b/Modules/PhotoacousticsAlgorithms/source/filters/mitkBeamformingFilter.cpp index ad4ee3b1a8..ae0b924a76 100644 --- a/Modules/PhotoacousticsAlgorithms/source/filters/mitkBeamformingFilter.cpp +++ b/Modules/PhotoacousticsAlgorithms/source/filters/mitkBeamformingFilter.cpp @@ -1,299 +1,294 @@ /*=================================================================== mitkBeamformingFilter 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 "mitkBeamformingFilter.h" #include "mitkBeamformingUtils.h" mitk::BeamformingFilter::BeamformingFilter(mitk::BeamformingSettings::Pointer settings) : m_OutputData(nullptr), m_InputData(nullptr), m_Conf(settings) { MITK_INFO << "Instantiating BeamformingFilter..."; this->SetNumberOfIndexedInputs(1); this->SetNumberOfRequiredInputs(1); m_ProgressHandle = [](int, std::string) {}; m_BeamformingOclFilter = mitk::PhotoacousticOCLBeamformingFilter::New(m_Conf); MITK_INFO << "Instantiating BeamformingFilter...[Done]"; } void mitk::BeamformingFilter::SetProgressHandle(std::function progressHandle) { m_ProgressHandle = progressHandle; } mitk::BeamformingFilter::~BeamformingFilter() { MITK_INFO << "Destructed BeamformingFilter"; } void mitk::BeamformingFilter::GenerateInputRequestedRegion() { Superclass::GenerateInputRequestedRegion(); mitk::Image* output = this->GetOutput(); mitk::Image* input = const_cast (this->GetInput()); if (!output->IsInitialized()) { return; } input->SetRequestedRegionToLargestPossibleRegion(); } void mitk::BeamformingFilter::GenerateOutputInformation() { mitk::Image::ConstPointer input = this->GetInput(); mitk::Image::Pointer output = this->GetOutput(); if ((output->IsInitialized()) && (this->GetMTime() <= m_TimeOfHeaderInitialization.GetMTime())) return; unsigned int dim[] = { m_Conf->GetReconstructionLines(), m_Conf->GetSamplesPerLine(), input->GetDimension(2) }; output->Initialize(mitk::MakeScalarPixelType(), 3, dim); mitk::Vector3D spacing; spacing[0] = m_Conf->GetPitchInMeters() * m_Conf->GetTransducerElements() * 1000 / m_Conf->GetReconstructionLines(); spacing[1] = (m_Conf->GetTimeSpacing() * m_Conf->GetInputDim()[1]) / (2-(int)m_Conf->GetIsPhotoacousticImage()) * m_Conf->GetSpeedOfSound() * 1000 / m_Conf->GetSamplesPerLine(); spacing[2] = 1; output->GetGeometry()->SetSpacing(spacing); output->GetGeometry()->Modified(); output->SetPropertyList(input->GetPropertyList()->Clone()); m_TimeOfHeaderInitialization.Modified(); } void mitk::BeamformingFilter::GenerateData() { mitk::Image::Pointer input = this->GetInput(); if (!(input->GetPixelType().GetTypeAsString() == "scalar (float)" || input->GetPixelType().GetTypeAsString() == " (float)")) { MITK_ERROR << "Pixel type of input needs to be float for this filter to work."; mitkThrow() << "Pixel type of input needs to be float for this filter to work."; } GenerateOutputInformation(); mitk::Image::Pointer output = this->GetOutput(); if (!output->IsInitialized()) return; auto begin = std::chrono::high_resolution_clock::now(); // debbuging the performance... if (!m_Conf->GetUseGPU()) { int progInterval = output->GetDimension(2) / 20 > 1 ? output->GetDimension(2) / 20 : 1; // the interval at which we update the gui progress bar float inputDim[2] = { (float)input->GetDimension(0), (float)input->GetDimension(1) }; float outputDim[2] = { (float)output->GetDimension(0), (float)output->GetDimension(1) }; for (unsigned int i = 0; i < output->GetDimension(2); ++i) // seperate Slices should get Beamforming seperately applied { mitk::ImageReadAccessor inputReadAccessor(input, input->GetSliceData(i)); m_InputData = (float*)inputReadAccessor.GetData(); m_OutputData = new float[m_Conf->GetReconstructionLines()*m_Conf->GetSamplesPerLine()]; // fill the image with zeros for (int l = 0; l < outputDim[0]; ++l) { for (int s = 0; s < outputDim[1]; ++s) { m_OutputData[l*(short)outputDim[1] + s] = 0; } } std::thread *threads = new std::thread[(short)outputDim[0]]; // every line will be beamformed in a seperate thread if (m_Conf->GetAlgorithm() == BeamformingSettings::BeamformingAlgorithm::DAS) { if (m_Conf->GetDelayCalculationMethod() == BeamformingSettings::DelayCalc::QuadApprox) { for (short line = 0; line < outputDim[0]; ++line) { threads[line] = std::thread(&BeamformingUtils::DASQuadraticLine, m_InputData, m_OutputData, inputDim, outputDim, line, m_Conf); } } else if (m_Conf->GetDelayCalculationMethod() == BeamformingSettings::DelayCalc::Spherical) { for (short line = 0; line < outputDim[0]; ++line) { threads[line] = std::thread(&BeamformingUtils::DASSphericalLine, m_InputData, m_OutputData, inputDim, outputDim, line, m_Conf); } } } else if (m_Conf->GetAlgorithm() == BeamformingSettings::BeamformingAlgorithm::DMAS) { if (m_Conf->GetDelayCalculationMethod() == BeamformingSettings::DelayCalc::QuadApprox) { for (short line = 0; line < outputDim[0]; ++line) { threads[line] = std::thread(&BeamformingUtils::DMASQuadraticLine, m_InputData, m_OutputData, inputDim, outputDim, line, m_Conf); } } else if (m_Conf->GetDelayCalculationMethod() == BeamformingSettings::DelayCalc::Spherical) { for (short line = 0; line < outputDim[0]; ++line) { threads[line] = std::thread(&BeamformingUtils::DMASSphericalLine, m_InputData, m_OutputData, inputDim, outputDim, line, m_Conf); } } } else if (m_Conf->GetAlgorithm() == BeamformingSettings::BeamformingAlgorithm::sDMAS) { if (m_Conf->GetDelayCalculationMethod() == BeamformingSettings::DelayCalc::QuadApprox) { for (short line = 0; line < outputDim[0]; ++line) { threads[line] = std::thread(&BeamformingUtils::sDMASQuadraticLine, m_InputData, m_OutputData, inputDim, outputDim, line, m_Conf); } } else if (m_Conf->GetDelayCalculationMethod() == BeamformingSettings::DelayCalc::Spherical) { for (short line = 0; line < outputDim[0]; ++line) { threads[line] = std::thread(&BeamformingUtils::sDMASSphericalLine, m_InputData, m_OutputData, inputDim, outputDim, line, m_Conf); } } } // wait for all lines to finish for (short line = 0; line < outputDim[0]; ++line) { threads[line].join(); } output->SetSlice(m_OutputData, i); if (i % progInterval == 0) m_ProgressHandle((int)((i + 1) / (float)output->GetDimension(2) * 100), "performing reconstruction"); delete[] m_OutputData; m_OutputData = nullptr; m_InputData = nullptr; } } #if defined(PHOTOACOUSTICS_USE_GPU) || DOXYGEN else { try { // first, we check whether the data is float, other formats are unsupported if (!(input->GetPixelType().GetTypeAsString() == "scalar (float)" || input->GetPixelType().GetTypeAsString() == " (float)")) { MITK_ERROR << "Pixel type is not float, abort"; mitkThrow() << "Pixel type is not float, abort"; } unsigned long availableMemory = m_BeamformingOclFilter->GetDeviceMemory(); - unsigned int batchSize = 16; + unsigned int batchSize = m_Conf->GetGPUBatchSize(); unsigned int batches = (unsigned int)((float)input->GetDimension(2) / batchSize) + (input->GetDimension(2) % batchSize > 0); unsigned int batchDim[] = { input->GetDimension(0), input->GetDimension(1), batchSize }; unsigned int batchDimLast[] = { input->GetDimension(0), input->GetDimension(1), input->GetDimension(2) % batchSize }; // the following safeguard is probably only needed for absurdly small GPU memory - for (batchSize = 16; - (unsigned long)batchSize * + if((unsigned long)batchSize * ((unsigned long)(batchDim[0] * batchDim[1]) * 4 + // single input image (float) (unsigned long)(m_Conf->GetReconstructionLines() * m_Conf->GetSamplesPerLine()) * 4) // single output image (float) > availableMemory - (unsigned long)(m_Conf->GetReconstructionLines() / 2 * m_Conf->GetSamplesPerLine()) * 2 - // Delays buffer (unsigned short) (unsigned long)(m_Conf->GetReconstructionLines() * m_Conf->GetSamplesPerLine()) * 3 * 2 - // UsedLines buffer (unsigned short) - 50 * 1024 * 1024; // 50 MB buffer for local data, system purposes etc - --batchSize) + 50 * 1024 * 1024)// 50 MB buffer for local data, system purposes etc { - } - if (batchSize < 1) - { - MITK_ERROR << "device memory too small for GPU beamforming"; + MITK_ERROR << "device memory too small for GPU beamforming; try decreasing the batch size"; return; } mitk::ImageReadAccessor copy(input); for (unsigned int i = 0; i < batches; ++i) { - m_ProgressHandle(input->GetDimension(2) / batches * i, "performing reconstruction"); + m_ProgressHandle(100.f * (float)i / (float)batches, "performing reconstruction"); mitk::Image::Pointer inputBatch = mitk::Image::New(); unsigned int num_Slices = 1; if (i == batches - 1 && (input->GetDimension(2) % batchSize > 0)) { inputBatch->Initialize(mitk::MakeScalarPixelType(), 3, batchDimLast); num_Slices = batchDimLast[2]; } else { inputBatch->Initialize(mitk::MakeScalarPixelType(), 3, batchDim); num_Slices = batchDim[2]; } inputBatch->SetSpacing(input->GetGeometry()->GetSpacing()); inputBatch->SetImportVolume(&(((float*)copy.GetData())[input->GetDimension(0) * input->GetDimension(1) * batchSize * i])); - + m_BeamformingOclFilter->SetApodisation(m_Conf->GetApodizationFunction(), m_Conf->GetApodizationArraySize()); - m_BeamformingOclFilter->SetInput(inputBatch); + m_BeamformingOclFilter->SetInput(inputBatch, num_Slices); m_BeamformingOclFilter->Update(); void* out = m_BeamformingOclFilter->GetOutput(); for (unsigned int slice = 0; slice < num_Slices; ++slice) { output->SetImportSlice( &(((float*)out)[m_Conf->GetReconstructionLines() * m_Conf->GetSamplesPerLine() * slice]), batchSize * i + slice, 0, 0, mitk::Image::ImportMemoryManagementType::CopyMemory); } } } catch (mitk::Exception &e) { std::string errorMessage = "Caught unexpected exception "; errorMessage.append(e.what()); MITK_ERROR << errorMessage; float* dummyData = new float[m_Conf->GetReconstructionLines() * m_Conf->GetSamplesPerLine() * m_Conf->GetInputDim()[2]]; output->SetImportVolume(dummyData, 0, 0, mitk::Image::ImportMemoryManagementType::CopyMemory); } } #endif m_TimeOfHeaderInitialization.Modified(); auto end = std::chrono::high_resolution_clock::now(); MITK_INFO << "Beamforming of " << output->GetDimension(2) << " Images completed in " << ((float)std::chrono::duration_cast(end - begin).count()) / 1000000 << "ms" << std::endl; } diff --git a/Modules/PhotoacousticsAlgorithms/source/filters/mitkBeamformingSettings.cpp b/Modules/PhotoacousticsAlgorithms/source/filters/mitkBeamformingSettings.cpp index 6aa5f38981..d956ecfdcd 100644 --- a/Modules/PhotoacousticsAlgorithms/source/filters/mitkBeamformingSettings.cpp +++ b/Modules/PhotoacousticsAlgorithms/source/filters/mitkBeamformingSettings.cpp @@ -1,120 +1,100 @@ /*=================================================================== 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(std::string xmlFile) { MITK_ERROR << "Not implemented yet."; mitkThrow() << "Not implemented yet."; } mitk::BeamformingSettings::BeamformingSettings(float pitchInMeters, float speedOfSound, float timeSpacing, float angle, bool isPhotoacousticImage, unsigned int samplesPerLine, unsigned int reconstructionLines, - unsigned int upperCutoff, - bool partial, - unsigned int* cropBounds, unsigned int* inputDim, bool useGPU, + unsigned int GPUBatchSize, DelayCalc delayCalculationMethod, Apodization apod, unsigned int apodizationArraySize, - BeamformingAlgorithm algorithm, - bool useBP, - float BPHighPass, - float BPLowPass + BeamformingAlgorithm algorithm ) : m_PitchInMeters(pitchInMeters), m_SpeedOfSound(speedOfSound), m_TimeSpacing(timeSpacing), m_Angle(angle), m_IsPhotoacousticImage(isPhotoacousticImage), m_SamplesPerLine(samplesPerLine), m_ReconstructionLines(reconstructionLines), - m_UpperCutoff(upperCutoff), - m_Partial(partial), - m_CropBounds(cropBounds), m_UseGPU(useGPU), + m_GPUBatchSize(GPUBatchSize), m_DelayCalculationMethod(delayCalculationMethod), m_Apod(apod), m_ApodizationArraySize(apodizationArraySize), - m_Algorithm(algorithm), - m_UseBP(useBP), - m_BPHighPass(BPHighPass), - m_BPLowPass(BPLowPass) + m_Algorithm(algorithm) { 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; } - if (m_CropBounds == nullptr) - m_CropBounds = new unsigned int[2]{ 0, 0 }; - m_InputDim = new unsigned int[3]{ inputDim[0], inputDim[1], inputDim[2] }; m_TransducerElements = m_InputDim[0]; } 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_CropBounds != nullptr) - { - MITK_INFO << "Deleting crop bounds..."; - delete[] m_CropBounds; - MITK_INFO << "Deleting crop bounds...[Done]"; - } - if (m_InputDim != nullptr) { MITK_INFO << "Deleting input dim..."; delete[] m_InputDim; MITK_INFO << "Deleting input dim...[Done]"; } MITK_INFO << "Destructing beamforming settings...[Done]"; } diff --git a/Modules/PhotoacousticsAlgorithms/test/mitkBeamformingFilterTest.cpp b/Modules/PhotoacousticsAlgorithms/test/mitkBeamformingFilterTest.cpp index dd54563f5c..54e4504acf 100644 --- a/Modules/PhotoacousticsAlgorithms/test/mitkBeamformingFilterTest.cpp +++ b/Modules/PhotoacousticsAlgorithms/test/mitkBeamformingFilterTest.cpp @@ -1,283 +1,279 @@ /*=================================================================== 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 #include #include #include #include #include #include #include class SyntheticPAImageData { public: SyntheticPAImageData(float spacing_x, float spacing_y, unsigned int samples, unsigned int num_transducers, float speedOfSound) { m_Spacing_x = spacing_x; m_Spacing_y = spacing_y; m_TransducerElements = num_transducers; m_Samples = samples; m_SpeedOfSound = speedOfSound; m_Data = new float[m_Samples*m_TransducerElements]; for (size_t i = 0; i < m_Samples * m_TransducerElements; ++i) { m_Data[i] = 0; } } ~SyntheticPAImageData() { delete[] m_Data; } void AddWave(float origin_depth, float origin_x, float base_value= 10000) { AddLine(origin_depth, origin_x, base_value); AddLine(origin_depth + 0.0001f, origin_x, -base_value); } const float* GetData() { return (const float*)m_Data; } private: void AddLine(float origin_depth, unsigned int origin_x, float base_value = 10000) { for (unsigned int x = 0; x < m_TransducerElements; ++x) { float distance = std::abs((int)x - (int)origin_x); float delay_in_seconds = std::sqrt(std::pow(origin_depth, 2) + std::pow(distance*m_Spacing_x, 2)) / m_SpeedOfSound; int pixels = std::round(delay_in_seconds / m_Spacing_y); for (int index = -4; index < 9; ++index) { if ((int)pixels + index < (int)m_Samples && (int)pixels + index > 0) { m_Data[(size_t)(x + (pixels + index)*m_TransducerElements)] += base_value / std::sqrt(distance + 1); } } } } float* m_Data; float m_Spacing_x; float m_Spacing_y; unsigned int m_TransducerElements; unsigned int m_Samples; float m_SpeedOfSound; }; class mitkBeamformingFilterTestSuite : public mitk::TestFixture { CPPUNIT_TEST_SUITE(mitkBeamformingFilterTestSuite); MITK_TEST(testBeamformingCPU_DAS); MITK_TEST(testBeamformingGPU_DAS); MITK_TEST(testBeamformingCPU_DMAS); MITK_TEST(testBeamformingGPU_DMAS); MITK_TEST(testBeamformingCPU_sDMAS); MITK_TEST(testBeamformingGPU_sDMAS); CPPUNIT_TEST_SUITE_END(); private: mitk::BeamformingFilter::Pointer m_BeamformingFilter; const unsigned int NUM_ITERATIONS = 15; const unsigned int SAMPLES = 5000 * 2; const unsigned int RECONSTRUCTED_SAMPLES = 2048; const unsigned int ELEMENTS = 128; const unsigned int RECONSTRUCTED_LINES = 128; const float SPEED_OF_SOUND = 1540; // m/s const float SPACING_X = 0.3; // mm const float SPACING_Y = 0.00625 / 2; // us + const unsigned int GPU_BATCH_SIZE = 16; public: void setUp() override { } void test() { std::random_device r; std::default_random_engine randGen(r()); float maxDepthInMeters = SAMPLES * (SPACING_Y / 1000000) * SPEED_OF_SOUND / 2; std::uniform_real_distribution randDistrDepth(maxDepthInMeters*0.1, maxDepthInMeters*0.8); for (unsigned int iteration = 0; iteration < NUM_ITERATIONS; ++iteration) { // create some synthetic input data float depth_in_meters1 = randDistrDepth(randGen); float depth_in_meters2 = randDistrDepth(randGen); float depth_in_meters3 = randDistrDepth(randGen); unsigned int element1 = 29; unsigned int element2 = 63; unsigned int element3 = 98; SyntheticPAImageData image(SPACING_X / 1000.f, SPACING_Y / 1000000.f, SAMPLES, ELEMENTS, SPEED_OF_SOUND); image.AddWave(depth_in_meters1, 29); image.AddWave(depth_in_meters2, 63); image.AddWave(depth_in_meters3, 98); mitk::Image::Pointer inputImage = mitk::Image::New(); unsigned int dimension[3]{ ELEMENTS, SAMPLES, 1 }; inputImage->Initialize(mitk::MakeScalarPixelType(), 3, dimension); mitk::Vector3D spacing; spacing[0] = SPACING_X; spacing[1] = SPACING_Y; spacing[2] = 1; inputImage->SetSpacing(spacing); inputImage->SetImportVolume((const void*)image.GetData(), mitk::Image::RtlCopyMemory); // setup the beamforming filter m_BeamformingFilter->SetInput(inputImage); m_BeamformingFilter->Update(); mitk::Image::Pointer outputImage = m_BeamformingFilter->GetOutput(); mitk::ImageReadAccessor readAccess(outputImage); const float* outputData = (const float*)readAccess.GetData(); unsigned int pos1[3] = { element1, (unsigned int)std::round(depth_in_meters1 * 1000.f / outputImage->GetGeometry()->GetSpacing()[1]), 0 }; unsigned int pos2[3] = { element2, (unsigned int)std::round(depth_in_meters2 * 1000.f / outputImage->GetGeometry()->GetSpacing()[1]), 0 }; unsigned int pos3[3] = { element3, (unsigned int)std::round(depth_in_meters3 * 1000.f / outputImage->GetGeometry()->GetSpacing()[1]), 0 }; double average = 0; for (unsigned int i = 0; i < RECONSTRUCTED_LINES*RECONSTRUCTED_SAMPLES; ++i) { average += outputData[i] / (RECONSTRUCTED_LINES*RECONSTRUCTED_SAMPLES); } CPPUNIT_ASSERT_MESSAGE(std::string("Iteration " + std::to_string(iteration) + ": first point source incorrectly reconstructed; should be > average*100, is " + std::to_string(abs(outputData[pos1[0] + pos1[1] * RECONSTRUCTED_LINES]))) + " < " + std::to_string(average) + "*100" , abs(outputData[pos1[0] + pos1[1]*RECONSTRUCTED_LINES] / average) > 100); CPPUNIT_ASSERT_MESSAGE(std::string("Iteration " + std::to_string(iteration) + ": second point source incorrectly reconstructed; should be > average*100, is " + std::to_string(abs(outputData[pos2[0] + pos2[1] * RECONSTRUCTED_LINES]))) + " < " + std::to_string(average) + "*100" , abs(outputData[pos2[0] + pos2[1] * RECONSTRUCTED_LINES] / average) > 100); CPPUNIT_ASSERT_MESSAGE(std::string("Iteration " + std::to_string(iteration) + ": third point source incorrectly reconstructed; should be > average*100, is " + std::to_string(abs(outputData[pos3[0] + pos3[1] * RECONSTRUCTED_LINES]))) + " < " + std::to_string(average) + "*100" , abs(outputData[pos3[0] + pos3[1] * RECONSTRUCTED_LINES] / average) > 100); } } mitk::BeamformingSettings::Pointer createConfig(bool UseGPU, unsigned int* inputDim, mitk::BeamformingSettings::BeamformingAlgorithm alg) { return mitk::BeamformingSettings::New(SPACING_X / 1000, SPEED_OF_SOUND, SPACING_Y / 1000000, 27.f, true, RECONSTRUCTED_SAMPLES, RECONSTRUCTED_LINES, - 0, - false, - nullptr, inputDim, UseGPU, + GPU_BATCH_SIZE, mitk::BeamformingSettings::DelayCalc::Spherical, mitk::BeamformingSettings::Apodization::Box, ELEMENTS * 2, - alg, - false, - 0, - 0); + alg); } void testBeamformingCPU_DAS() { MITK_INFO << "Started DAS test on CPU"; unsigned int* inputDim = new unsigned int[3]; inputDim[0] = ELEMENTS; inputDim[1] = SAMPLES; inputDim[2] = 1; m_BeamformingFilter = mitk::BeamformingFilter::New(createConfig(false, inputDim, mitk::BeamformingSettings::BeamformingAlgorithm::DAS)); test(); delete[] inputDim; } void testBeamformingGPU_DAS() { MITK_INFO << "Started DAS test on GPU"; unsigned int* inputDim = new unsigned int[3]; inputDim[0] = ELEMENTS; inputDim[1] = SAMPLES; inputDim[2] = 1; m_BeamformingFilter = mitk::BeamformingFilter::New(createConfig(true, inputDim, mitk::BeamformingSettings::BeamformingAlgorithm::DAS)); test(); delete[] inputDim; } void testBeamformingCPU_sDMAS() { MITK_INFO << "Started sDMAS test on CPU"; unsigned int* inputDim = new unsigned int[3]; inputDim[0] = ELEMENTS; inputDim[1] = SAMPLES; inputDim[2] = 1; m_BeamformingFilter = mitk::BeamformingFilter::New(createConfig(false, inputDim, mitk::BeamformingSettings::BeamformingAlgorithm::sDMAS)); test(); delete[] inputDim; } void testBeamformingGPU_sDMAS() { MITK_INFO << "Started sDMAS test on GPU"; unsigned int* inputDim = new unsigned int[3]; inputDim[0] = ELEMENTS; inputDim[1] = SAMPLES; inputDim[2] = 1; m_BeamformingFilter = mitk::BeamformingFilter::New(createConfig(true, inputDim, mitk::BeamformingSettings::BeamformingAlgorithm::sDMAS)); test(); delete[] inputDim; } void testBeamformingCPU_DMAS() { MITK_INFO << "Started DMAS test on CPU"; unsigned int* inputDim = new unsigned int[3]; inputDim[0] = ELEMENTS; inputDim[1] = SAMPLES; inputDim[2] = 1; m_BeamformingFilter = mitk::BeamformingFilter::New(createConfig(false, inputDim, mitk::BeamformingSettings::BeamformingAlgorithm::DMAS)); test(); delete[] inputDim; } void testBeamformingGPU_DMAS() { MITK_INFO << "Started DMAS test on GPU"; unsigned int* inputDim = new unsigned int[3]; inputDim[0] = ELEMENTS; inputDim[1] = SAMPLES; inputDim[2] = 1; m_BeamformingFilter = mitk::BeamformingFilter::New(createConfig(true, inputDim, mitk::BeamformingSettings::BeamformingAlgorithm::DMAS)); test(); delete[] inputDim; } }; MITK_TEST_SUITE_REGISTRATION(mitkBeamformingFilter) diff --git a/Modules/PhotoacousticsAlgorithms/test/mitkPAFilterServiceTest.cpp b/Modules/PhotoacousticsAlgorithms/test/mitkPAFilterServiceTest.cpp index 9ce85b81ce..17820cf95f 100644 --- a/Modules/PhotoacousticsAlgorithms/test/mitkPAFilterServiceTest.cpp +++ b/Modules/PhotoacousticsAlgorithms/test/mitkPAFilterServiceTest.cpp @@ -1,127 +1,122 @@ /*=================================================================== 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 #include #include #include #include #include // us #include #include #include #include #include #include #include class mitkPAFilterServiceTestSuite : public mitk::TestFixture { CPPUNIT_TEST_SUITE(mitkPAFilterServiceTestSuite); MITK_TEST(testRunning); CPPUNIT_TEST_SUITE_END(); private: mitk::PhotoacousticFilterService::Pointer m_PhotoacousticFilterService; mitk::BeamformingSettings::Pointer m_BeamformingSettings; unsigned int* inputDimensions; const int xDim = 16; const int yDim = 128; const int length = yDim * xDim; public: void setUp() override { m_PhotoacousticFilterService = mitk::PhotoacousticFilterService::New(); m_BeamformingSettings = CreateBeamformingSettings(); } mitk::BeamformingSettings::Pointer CreateBeamformingSettings() { inputDimensions = new unsigned int[length]; inputDimensions[0] = yDim; inputDimensions[1] = xDim; mitk::BeamformingSettings::Pointer outputSettings = mitk::BeamformingSettings::New( (float)(0.3 / 1000), (float)(1500), (float)(0.0125 / 1000000), 27, true, 3000, 128, - 0, - false, - (unsigned int*) nullptr, inputDimensions, false, + 16, mitk::BeamformingSettings::DelayCalc::Spherical, mitk::BeamformingSettings::Apodization::Box, 128, - mitk::BeamformingSettings::BeamformingAlgorithm::DAS, - false, - 0.0f, - 0.0f); + mitk::BeamformingSettings::BeamformingAlgorithm::DAS); return outputSettings; } void testRunning() { float* testArray = new float[length]; for (int i = 0; i < length; ++i) { testArray[i] = 0; } mitk::PixelType pixelType = mitk::MakeScalarPixelType(); mitk::Image::Pointer testImage = mitk::Image::New(); testImage->Initialize(pixelType, 2, inputDimensions); testImage->SetImportSlice(testArray, 0, 0, 0, mitk::Image::ImportMemoryManagementType::CopyMemory); delete[] testArray; mitk::ImageReadAccessor readAccessInput(testImage); const float* inputArray = (const float*)readAccessInput.GetData(); for (int i = 0; i < length; ++i) { CPPUNIT_ASSERT_MESSAGE(std::string("Input array already not correct: " + std::to_string(inputArray[i])), abs(inputArray[i]) < 1e-5f); } auto output = m_PhotoacousticFilterService->ApplyBeamforming(testImage, m_BeamformingSettings); mitk::ImageReadAccessor readAccess(output); const float* outputArray = (const float*)readAccess.GetData(); for (int i = 0; i < length; ++i) { CPPUNIT_ASSERT_MESSAGE(std::string("Output array not correct: " + std::to_string(abs(outputArray[i]))), abs(outputArray[i]) < 1e-5f); } } void tearDown() override { m_PhotoacousticFilterService = nullptr; delete[] inputDimensions; } }; MITK_TEST_SUITE_REGISTRATION(mitkPAFilterService) diff --git a/Plugins/org.mitk.gui.qt.photoacoustics.imageprocessing/src/internal/PAImageProcessing.cpp b/Plugins/org.mitk.gui.qt.photoacoustics.imageprocessing/src/internal/PAImageProcessing.cpp index 6fb103678e..d35d5508b0 100644 --- a/Plugins/org.mitk.gui.qt.photoacoustics.imageprocessing/src/internal/PAImageProcessing.cpp +++ b/Plugins/org.mitk.gui.qt.photoacoustics.imageprocessing/src/internal/PAImageProcessing.cpp @@ -1,1021 +1,1019 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ // Blueberry #include #include // Qmitk #include "PAImageProcessing.h" // Qt #include #include #include #include //mitk image #include #include "mitkPhotoacousticFilterService.h" #include "mitkCastToFloatImageFilter.h" #include "mitkBeamformingFilter.h" //other #include #include #include const std::string PAImageProcessing::VIEW_ID = "org.mitk.views.paimageprocessing"; PAImageProcessing::PAImageProcessing() : m_ResampleSpacing(0), m_UseLogfilter(false), m_FilterBank(mitk::PhotoacousticFilterService::New()) { qRegisterMetaType(); qRegisterMetaType(); } void PAImageProcessing::SetFocus() { m_Controls.buttonApplyBModeFilter->setFocus(); } void PAImageProcessing::CreateQtPartControl(QWidget *parent) { // create GUI widgets from the Qt Designer's .ui file m_Controls.setupUi(parent); connect(m_Controls.buttonApplyBModeFilter, SIGNAL(clicked()), this, SLOT(StartBmodeThread())); connect(m_Controls.DoResampling, SIGNAL(clicked()), this, SLOT(UseResampling())); connect(m_Controls.Logfilter, SIGNAL(clicked()), this, SLOT(UseLogfilter())); connect(m_Controls.ResamplingValue, SIGNAL(valueChanged(double)), this, SLOT(SetResampling())); connect(m_Controls.buttonApplyBeamforming, SIGNAL(clicked()), this, SLOT(StartBeamformingThread())); connect(m_Controls.buttonApplyCropFilter, SIGNAL(clicked()), this, SLOT(StartCropThread())); connect(m_Controls.buttonApplyBandpass, SIGNAL(clicked()), this, SLOT(StartBandpassThread())); connect(m_Controls.UseImageSpacing, SIGNAL(clicked()), this, SLOT(UseImageSpacing())); connect(m_Controls.ScanDepth, SIGNAL(valueChanged(double)), this, SLOT(UpdateImageInfo())); connect(m_Controls.SpeedOfSound, SIGNAL(valueChanged(double)), this, SLOT(UpdateImageInfo())); connect(m_Controls.SpeedOfSound, SIGNAL(valueChanged(double)), this, SLOT(ChangedSOSBeamforming())); connect(m_Controls.BPSpeedOfSound, SIGNAL(valueChanged(double)), this, SLOT(ChangedSOSBandpass())); connect(m_Controls.Samples, SIGNAL(valueChanged(int)), this, SLOT(UpdateImageInfo())); connect(m_Controls.UseImageSpacing, SIGNAL(clicked()), this, SLOT(UpdateImageInfo())); connect(m_Controls.boundLow, SIGNAL(valueChanged(int)), this, SLOT(LowerSliceBoundChanged())); connect(m_Controls.boundHigh, SIGNAL(valueChanged(int)), this, SLOT(UpperSliceBoundChanged())); connect(m_Controls.Partial, SIGNAL(clicked()), this, SLOT(SliceBoundsEnabled())); connect(m_Controls.BatchProcessing, SIGNAL(clicked()), this, SLOT(BatchProcessing())); connect(m_Controls.StepBeamforming, SIGNAL(clicked()), this, SLOT(UpdateSaveBoxes())); connect(m_Controls.StepCropping, SIGNAL(clicked()), this, SLOT(UpdateSaveBoxes())); connect(m_Controls.StepBandpass, SIGNAL(clicked()), this, SLOT(UpdateSaveBoxes())); connect(m_Controls.StepBMode, SIGNAL(clicked()), this, SLOT(UpdateSaveBoxes())); UpdateSaveBoxes(); m_Controls.DoResampling->setChecked(false); m_Controls.ResamplingValue->setEnabled(false); m_Controls.progressBar->setMinimum(0); m_Controls.progressBar->setMaximum(100); m_Controls.progressBar->setVisible(false); m_Controls.UseImageSpacing->setToolTip("Image spacing of y-Axis must be in us, x-Axis in mm."); m_Controls.UseImageSpacing->setToolTipDuration(5000); m_Controls.ProgressInfo->setVisible(false); m_Controls.UseBP->hide(); m_Controls.UseGPUBmode->hide(); #ifndef PHOTOACOUSTICS_USE_GPU m_Controls.UseGPUBf->setEnabled(false); m_Controls.UseGPUBf->setChecked(false); m_Controls.UseGPUBmode->setEnabled(false); m_Controls.UseGPUBmode->setChecked(false); #endif UseImageSpacing(); } void PAImageProcessing::ChangedSOSBandpass() { m_Controls.SpeedOfSound->setValue(m_Controls.BPSpeedOfSound->value()); } void PAImageProcessing::ChangedSOSBeamforming() { m_Controls.BPSpeedOfSound->setValue(m_Controls.SpeedOfSound->value()); } std::vector splitpath( const std::string& str , const std::set delimiters) { std::vector result; char const* pch = str.c_str(); char const* start = pch; for (; *pch; ++pch) { if (delimiters.find(*pch) != delimiters.end()) { if (start != pch) { std::string str(start, pch); result.push_back(str); } else { result.push_back(""); } start = pch + 1; } } result.push_back(start); return result; } void PAImageProcessing::UpdateSaveBoxes() { if (m_Controls.StepBeamforming->isChecked()) m_Controls.SaveBeamforming->setEnabled(true); else m_Controls.SaveBeamforming->setEnabled(false); if (m_Controls.StepCropping->isChecked()) m_Controls.SaveCropping->setEnabled(true); else m_Controls.SaveCropping->setEnabled(false); if (m_Controls.StepBandpass->isChecked()) m_Controls.SaveBandpass->setEnabled(true); else m_Controls.SaveBandpass->setEnabled(false); if (m_Controls.StepBMode->isChecked()) m_Controls.SaveBMode->setEnabled(true); else m_Controls.SaveBMode->setEnabled(false); } void PAImageProcessing::BatchProcessing() { QFileDialog LoadDialog(nullptr, "Select Files to be processed"); LoadDialog.setFileMode(QFileDialog::FileMode::ExistingFiles); LoadDialog.setNameFilter(tr("Images (*.nrrd)")); LoadDialog.setViewMode(QFileDialog::Detail); QStringList fileNames; if (LoadDialog.exec()) fileNames = LoadDialog.selectedFiles(); QString saveDir = QFileDialog::getExistingDirectory(nullptr, tr("Select Directory To Save To"), "", QFileDialog::ShowDirsOnly | QFileDialog::DontResolveSymlinks); DisableControls(); std::set delims{ '/' }; bool doSteps[] = { m_Controls.StepBeamforming->isChecked(), m_Controls.StepCropping->isChecked() , m_Controls.StepBandpass->isChecked(), m_Controls.StepBMode->isChecked() }; bool saveSteps[] = { m_Controls.SaveBeamforming->isChecked(), m_Controls.SaveCropping->isChecked() , m_Controls.SaveBandpass->isChecked(), m_Controls.SaveBMode->isChecked() }; for (int fileNumber = 0; fileNumber < fileNames.size(); ++fileNumber) { m_Controls.progressBar->setValue(0); m_Controls.progressBar->setVisible(true); m_Controls.ProgressInfo->setVisible(true); m_Controls.ProgressInfo->setText("loading file"); QString filename = fileNames.at(fileNumber); auto split = splitpath(filename.toStdString(), delims); std::string imageName = split.at(split.size() - 1); // remove ".nrrd" imageName = imageName.substr(0, imageName.size() - 5); mitk::Image::Pointer image = mitk::IOUtil::Load(filename.toStdString().c_str()); auto BFconfig = CreateBeamformingSettings(image); // Beamforming if (doSteps[0]) { std::function progressHandle = [this](int progress, std::string progressInfo) { this->UpdateProgress(progress, progressInfo); }; m_Controls.progressBar->setValue(100); image = m_FilterBank->ApplyBeamforming(image, BFconfig, progressHandle); if (saveSteps[0]) { std::string saveFileName = saveDir.toStdString() + "/" + imageName + " beamformed" + ".nrrd"; mitk::IOUtil::Save(image, saveFileName); } } // Cropping if (doSteps[1]) { m_Controls.ProgressInfo->setText("cropping image"); image = m_FilterBank->ApplyCropping(image, m_Controls.CutoffAbove->value(), m_Controls.CutoffBelow->value(), 0, 0, 0, 0); if (saveSteps[1]) { std::string saveFileName = saveDir.toStdString() + "/" + imageName + " cropped" + ".nrrd"; mitk::IOUtil::Save(image, saveFileName); } } // Bandpass if (doSteps[2]) { m_Controls.ProgressInfo->setText("applying bandpass"); float recordTime = image->GetDimension(1)*image->GetGeometry()->GetSpacing()[1] / 1000 / m_Controls.BPSpeedOfSound->value(); // add a safeguard so the program does not chrash when applying a Bandpass that reaches out of the bounds of the image float maxFrequency = 1 / (recordTime / image->GetDimension(1)) * image->GetDimension(1) / 2 / 2 / 1000; float BPHighPass = 1000000 * m_Controls.BPhigh->value(); // [Hz] float BPLowPass = maxFrequency - 1000000 * m_Controls.BPlow->value(); // [Hz] if (BPLowPass > maxFrequency && m_Controls.UseBP->isChecked()) { QMessageBox Msgbox; Msgbox.setText("LowPass too low, disabled it."); Msgbox.exec(); BPLowPass = 0; } if (BPLowPass < 0 && m_Controls.UseBP->isChecked()) { QMessageBox Msgbox; Msgbox.setText("LowPass too high, disabled it."); Msgbox.exec(); BPLowPass = 0; } if (BPHighPass > maxFrequency && m_Controls.UseBP->isChecked()) { QMessageBox Msgbox; Msgbox.setText("HighPass too high, disabled it."); Msgbox.exec(); BPHighPass = 0; } if (BPHighPass > maxFrequency - BPLowPass) { QMessageBox Msgbox; Msgbox.setText("HighPass higher than LowPass, disabled both."); Msgbox.exec(); BPHighPass = 0; BPLowPass = 0; } image = m_FilterBank->ApplyBandpassFilter(image, BPHighPass, BPLowPass, m_Controls.BPFalloffHigh->value(), m_Controls.BPFalloffLow->value()); if (saveSteps[2]) { std::string saveFileName = saveDir.toStdString() + "/" + imageName + " bandpassed" + ".nrrd"; mitk::IOUtil::Save(image, saveFileName); } } // Bmode if (doSteps[3]) { m_Controls.ProgressInfo->setText("applying bmode filter"); if (m_Controls.BModeMethod->currentText() == "Absolute Filter") image = m_FilterBank->ApplyBmodeFilter(image, mitk::PhotoacousticFilterService::BModeMethod::Abs, m_UseLogfilter); else if (m_Controls.BModeMethod->currentText() == "Envelope Detection") image = m_FilterBank->ApplyBmodeFilter(image, mitk::PhotoacousticFilterService::BModeMethod::EnvelopeDetection, m_UseLogfilter); double desiredSpacing[2]{ image->GetGeometry()->GetSpacing()[0], m_ResampleSpacing }; image = m_FilterBank->ApplyResampling(image, desiredSpacing); if (saveSteps[3]) { std::string saveFileName = saveDir.toStdString() + "/" + imageName + " bmode" + ".nrrd"; mitk::IOUtil::Save(image, saveFileName); } } m_Controls.progressBar->setVisible(false); } EnableControls(); } void PAImageProcessing::StartBeamformingThread() { QList nodes = this->GetDataManagerSelection(); if (nodes.empty()) return; mitk::DataStorage::Pointer storage = this->GetDataStorage(); mitk::DataNode::Pointer node = nodes.front(); if (!node) { // Nothing selected. Inform the user and return QMessageBox::information(NULL, "Template", "Please load and select an image before starting image processing."); return; } mitk::BaseData* data = node->GetData(); if (data) { // test if this data item is an image or not (could also be a surface or something totally different) mitk::Image* image = dynamic_cast(data); if (image) { auto BFconfig = CreateBeamformingSettings(image); std::stringstream message; std::string name; message << "Performing beamforming for image "; if (node->GetName(name)) { // a property called "name" was found for this DataNode message << "'" << name << "'"; m_OldNodeName = name; } else m_OldNodeName = " "; message << "."; MITK_INFO << message.str(); m_Controls.progressBar->setValue(0); m_Controls.progressBar->setVisible(true); m_Controls.ProgressInfo->setVisible(true); m_Controls.ProgressInfo->setText("started"); DisableControls(); BeamformingThread *thread = new BeamformingThread(); connect(thread, &BeamformingThread::result, this, &PAImageProcessing::HandleResults); connect(thread, &BeamformingThread::updateProgress, this, &PAImageProcessing::UpdateProgress); connect(thread, &BeamformingThread::finished, thread, &QObject::deleteLater); thread->setConfig(BFconfig); thread->setInputImage(image); thread->setFilterBank(m_FilterBank); MITK_INFO << "Started new thread for Beamforming"; thread->start(); } } } void PAImageProcessing::HandleResults(mitk::Image::Pointer image, std::string nameExtension) { MITK_INFO << "Handling results..."; auto newNode = mitk::DataNode::New(); newNode->SetData(image); newNode->SetName(m_OldNodeName + nameExtension); // update level window for the current dynamic range mitk::LevelWindow levelWindow; newNode->GetLevelWindow(levelWindow); levelWindow.SetAuto(image, true, true); newNode->SetLevelWindow(levelWindow); // add new node to data storage this->GetDataStorage()->Add(newNode); // disable progress bar m_Controls.progressBar->setVisible(false); m_Controls.ProgressInfo->setVisible(false); EnableControls(); // update rendering mitk::RenderingManager::GetInstance()->InitializeViews(image->GetGeometry(), mitk::RenderingManager::REQUEST_UPDATE_ALL, true); mitk::RenderingManager::GetInstance()->RequestUpdateAll(); MITK_INFO << "Handling results...[Done]"; } void PAImageProcessing::StartBmodeThread() { QList nodes = this->GetDataManagerSelection(); if (nodes.empty()) return; mitk::DataStorage::Pointer storage = this->GetDataStorage(); mitk::DataNode::Pointer node = nodes.front(); if (!node) { // Nothing selected. Inform the user and return QMessageBox::information(NULL, "Template", "Please load and select an image before starting image processing."); return; } mitk::BaseData* data = node->GetData(); if (data) { // test if this data item is an image or not (could also be a surface or something totally different) mitk::Image* image = dynamic_cast(data); if (image) { std::stringstream message; std::string name; message << "Performing image processing for image "; if (node->GetName(name)) { // a property called "name" was found for this DataNode message << "'" << name << "'"; m_OldNodeName = name; } else m_OldNodeName = " "; message << "."; MITK_INFO << message.str(); DisableControls(); BmodeThread *thread = new BmodeThread(); connect(thread, &BmodeThread::result, this, &PAImageProcessing::HandleResults); connect(thread, &BmodeThread::finished, thread, &QObject::deleteLater); bool useGPU = m_Controls.UseGPUBmode->isChecked(); if (m_Controls.BModeMethod->currentText() == "Absolute Filter") thread->setConfig(m_UseLogfilter, m_ResampleSpacing, mitk::PhotoacousticFilterService::BModeMethod::Abs, useGPU); else if (m_Controls.BModeMethod->currentText() == "Envelope Detection") thread->setConfig(m_UseLogfilter, m_ResampleSpacing, mitk::PhotoacousticFilterService::BModeMethod::EnvelopeDetection, useGPU); thread->setInputImage(image); thread->setFilterBank(m_FilterBank); MITK_INFO << "Started new thread for Image Processing"; thread->start(); } } } void PAImageProcessing::StartCropThread() { QList nodes = this->GetDataManagerSelection(); if (nodes.empty()) return; mitk::DataStorage::Pointer storage = this->GetDataStorage(); mitk::DataNode::Pointer node = nodes.front(); if (!node) { // Nothing selected. Inform the user and return QMessageBox::information(NULL, "Template", "Please load and select an image before starting image cropping."); return; } mitk::BaseData* data = node->GetData(); if (data) { // test if this data item is an image or not (could also be a surface or something totally different) mitk::Image* image = dynamic_cast(data); if (image) { std::stringstream message; std::string name; message << "Performing image cropping for image "; if (node->GetName(name)) { // a property called "name" was found for this DataNode message << "'" << name << "'"; m_OldNodeName = name; } else m_OldNodeName = " "; message << "."; MITK_INFO << message.str(); DisableControls(); CropThread *thread = new CropThread(); connect(thread, &CropThread::result, this, &PAImageProcessing::HandleResults); connect(thread, &CropThread::finished, thread, &QObject::deleteLater); - thread->setConfig(m_Controls.CutoffAbove->value(), m_Controls.CutoffBelow->value(), 0, image->GetDimension(2) - 1); + if(m_Controls.Partial->isChecked()) + thread->setConfig(m_Controls.CutoffAbove->value(), m_Controls.CutoffBelow->value(), m_Controls.boundLow->value(), m_Controls.boundHigh->value()); + else + thread->setConfig(m_Controls.CutoffAbove->value(), m_Controls.CutoffBelow->value(), 0, image->GetDimension(2) - 1); + thread->setInputImage(image); thread->setFilterBank(m_FilterBank); MITK_INFO << "Started new thread for Image Cropping"; thread->start(); } } } void PAImageProcessing::StartBandpassThread() { QList nodes = this->GetDataManagerSelection(); if (nodes.empty()) return; mitk::DataStorage::Pointer storage = this->GetDataStorage(); mitk::DataNode::Pointer node = nodes.front(); if (!node) { // Nothing selected. Inform the user and return QMessageBox::information(NULL, "Template", "Please load and select an image before starting image cropping."); return; } mitk::BaseData* data = node->GetData(); if (data) { // test if this data item is an image or not (could also be a surface or something totally different) mitk::Image* image = dynamic_cast(data); if (image) { std::stringstream message; std::string name; message << "Performing Bandpass filter on image "; if (node->GetName(name)) { // a property called "name" was found for this DataNode message << "'" << name << "'"; m_OldNodeName = name; } else m_OldNodeName = " "; message << "."; MITK_INFO << message.str(); DisableControls(); BandpassThread *thread = new BandpassThread(); connect(thread, &BandpassThread::result, this, &PAImageProcessing::HandleResults); connect(thread, &BandpassThread::finished, thread, &QObject::deleteLater); float BPHighPass = 1000000.0f * m_Controls.BPhigh->value(); // [Now in Hz] float BPLowPass = 1000000.0f * m_Controls.BPlow->value(); // [Now in Hz] thread->setConfig(BPHighPass, BPLowPass, m_Controls.BPFalloffLow->value(), m_Controls.BPFalloffHigh->value(), 0); thread->setInputImage(image); thread->setFilterBank(m_FilterBank); MITK_INFO << "Started new thread for Bandpass filter"; thread->start(); } } } void PAImageProcessing::SliceBoundsEnabled() { if (!m_Controls.Partial->isChecked()) { m_Controls.boundLow->setEnabled(false); m_Controls.boundHigh->setEnabled(false); return; } else { m_Controls.boundLow->setEnabled(true); m_Controls.boundHigh->setEnabled(true); } } void PAImageProcessing::UpperSliceBoundChanged() { if (m_Controls.boundLow->value() > m_Controls.boundHigh->value()) { m_Controls.boundLow->setValue(m_Controls.boundHigh->value()); } } void PAImageProcessing::LowerSliceBoundChanged() { if (m_Controls.boundLow->value() > m_Controls.boundHigh->value()) { m_Controls.boundHigh->setValue(m_Controls.boundLow->value()); } } void PAImageProcessing::UpdateProgress(int progress, std::string progressInfo) { if (progress < 100) m_Controls.progressBar->setValue(progress); else m_Controls.progressBar->setValue(100); m_Controls.ProgressInfo->setText(progressInfo.c_str()); qApp->processEvents(); } void PAImageProcessing::PAMessageBox(std::string message) { if (0 != message.compare("noMessage")) { QMessageBox msgBox; msgBox.setText(message.c_str()); msgBox.exec(); } } void PAImageProcessing::UpdateImageInfo() { QList nodes = this->GetDataManagerSelection(); if (nodes.empty()) return; mitk::DataNode::Pointer node = nodes.front(); if (!node) { // Nothing selected return; } mitk::BaseData* data = node->GetData(); if (data) { // test if this data item is an image or not (could also be a surface or something totally different) mitk::Image* image = dynamic_cast(data); if (image) { // beamforming configs if (m_Controls.UseImageSpacing->isChecked()) { m_Controls.ElementCount->setValue(image->GetDimension(0)); m_Controls.Pitch->setValue(image->GetGeometry()->GetSpacing()[0]); } m_Controls.boundLow->setMaximum(image->GetDimension(2) - 1); m_Controls.boundHigh->setMaximum(image->GetDimension(2) - 1); float speedOfSound = m_Controls.SpeedOfSound->value(); // [m/s] std::stringstream frequency; float timeSpacing; if (m_Controls.UseImageSpacing->isChecked()) { timeSpacing = image->GetGeometry()->GetSpacing()[1] / 1000000.0f; MITK_INFO << "Calculated Scan Depth of " << (image->GetDimension(1)*image->GetGeometry()->GetSpacing()[1] / 1000000) * speedOfSound * 100 / 2 << "cm"; } else { timeSpacing = (2 * m_Controls.ScanDepth->value() / 1000 / speedOfSound) / image->GetDimension(1); } float maxFrequency = (1 / timeSpacing) * image->GetDimension(1) / 2 / 2 / 1000; frequency << maxFrequency / 1000000; //[MHz] frequency << "MHz"; m_Controls.BPhigh->setMaximum(maxFrequency / 1000000); m_Controls.BPlow->setMaximum(maxFrequency / 1000000); frequency << " is the maximal allowed frequency for the selected image."; m_Controls.BPhigh->setToolTip(frequency.str().c_str()); m_Controls.BPlow->setToolTip(frequency.str().c_str()); m_Controls.BPhigh->setToolTipDuration(5000); m_Controls.BPlow->setToolTipDuration(5000); } } } void PAImageProcessing::OnSelectionChanged(berry::IWorkbenchPart::Pointer /*source*/, const QList& nodes) { // iterate all selected objects, adjust warning visibility foreach(mitk::DataNode::Pointer node, nodes) { if (node.IsNotNull() && dynamic_cast(node->GetData())) { m_Controls.labelWarning->setVisible(false); m_Controls.buttonApplyBModeFilter->setEnabled(true); m_Controls.labelWarning2->setVisible(false); m_Controls.buttonApplyCropFilter->setEnabled(true); m_Controls.labelWarning3->setVisible(false); m_Controls.buttonApplyBandpass->setEnabled(true); m_Controls.labelWarning4->setVisible(false); m_Controls.buttonApplyBeamforming->setEnabled(true); UpdateImageInfo(); return; } } m_Controls.labelWarning->setVisible(true); m_Controls.buttonApplyBModeFilter->setEnabled(false); m_Controls.labelWarning2->setVisible(true); m_Controls.buttonApplyCropFilter->setEnabled(false); m_Controls.labelWarning3->setVisible(true); m_Controls.buttonApplyBandpass->setEnabled(false); m_Controls.labelWarning4->setVisible(true); m_Controls.buttonApplyBeamforming->setEnabled(false); } void PAImageProcessing::UseResampling() { if (m_Controls.DoResampling->isChecked()) { m_Controls.ResamplingValue->setEnabled(true); m_ResampleSpacing = m_Controls.ResamplingValue->value(); } else { m_Controls.ResamplingValue->setEnabled(false); m_ResampleSpacing = 0; } } void PAImageProcessing::UseLogfilter() { m_UseLogfilter = m_Controls.Logfilter->isChecked(); } void PAImageProcessing::SetResampling() { m_ResampleSpacing = m_Controls.ResamplingValue->value(); } mitk::BeamformingSettings::Pointer PAImageProcessing::CreateBeamformingSettings(mitk::Image::Pointer image) { mitk::BeamformingSettings::BeamformingAlgorithm algorithm; if ("DAS" == m_Controls.BFAlgorithm->currentText()) algorithm = mitk::BeamformingSettings::BeamformingAlgorithm::DAS; else if ("DMAS" == m_Controls.BFAlgorithm->currentText()) algorithm = mitk::BeamformingSettings::BeamformingAlgorithm::DMAS; else if ("sDMAS" == m_Controls.BFAlgorithm->currentText()) algorithm = mitk::BeamformingSettings::BeamformingAlgorithm::sDMAS; mitk::BeamformingSettings::DelayCalc delay; if ("Quad. Approx." == m_Controls.DelayCalculation->currentText()) { delay = mitk::BeamformingSettings::DelayCalc::QuadApprox; } else if ("Spherical Wave" == m_Controls.DelayCalculation->currentText()) { delay = mitk::BeamformingSettings::DelayCalc::Spherical; } mitk::BeamformingSettings::Apodization apod; if ("Von Hann" == m_Controls.Apodization->currentText()) { apod = mitk::BeamformingSettings::Apodization::Hann; } else if ("Hamming" == m_Controls.Apodization->currentText()) { apod = mitk::BeamformingSettings::Apodization::Hamm; } else if ("Box" == m_Controls.Apodization->currentText()) { apod = mitk::BeamformingSettings::Apodization::Box; } float pitchInMeters = m_Controls.Pitch->value() / 1000; // [m] float speedOfSound = m_Controls.SpeedOfSound->value(); // [m/s] unsigned int samplesPerLine = m_Controls.Samples->value(); unsigned int reconstructionLines = m_Controls.Lines->value(); unsigned int apodizatonArraySize = m_Controls.Lines->value(); float angle = m_Controls.Angle->value(); // [deg] bool useBandpass = m_Controls.UseBP->isChecked(); bool useGPU = m_Controls.UseGPUBf->isChecked(); float timeSpacing; if (m_Controls.UseImageSpacing->isChecked()) { timeSpacing = image->GetGeometry()->GetSpacing()[1] / 1000000.0f; MITK_INFO << "Calculated Scan Depth of " << (image->GetDimension(1)*image->GetGeometry()->GetSpacing()[1] / 1000000) * speedOfSound * 100 / 2 << "cm"; } else { timeSpacing = (2 * m_Controls.ScanDepth->value() / 1000 / speedOfSound) / image->GetDimension(1); } bool isPAImage; if ("US Image" == m_Controls.ImageType->currentText()) { isPAImage = false; } else if ("PA Image" == m_Controls.ImageType->currentText()) { isPAImage = true; } - bool partial = m_Controls.Partial->isChecked(); - unsigned int cropBoundLow = m_Controls.boundLow->value(); - unsigned int cropBoundHigh = m_Controls.boundHigh->value(); - - unsigned int* cropBounds = new unsigned int[2]{ cropBoundLow, cropBoundHigh }; - return mitk::BeamformingSettings::New(pitchInMeters, - speedOfSound, timeSpacing, angle, isPAImage, samplesPerLine, reconstructionLines, 0, partial, - cropBounds, image->GetDimensions(), useGPU, delay, apod, + speedOfSound, timeSpacing, angle, isPAImage, samplesPerLine, reconstructionLines, + image->GetDimensions(), useGPU, 16, delay, apod, apodizatonArraySize, algorithm, useBandpass, 50, 50); } void PAImageProcessing::EnableControls() { m_Controls.BatchProcessing->setEnabled(true); m_Controls.StepBeamforming->setEnabled(true); m_Controls.StepBandpass->setEnabled(true); m_Controls.StepCropping->setEnabled(true); m_Controls.StepBMode->setEnabled(true); UpdateSaveBoxes(); m_Controls.DoResampling->setEnabled(true); UseResampling(); m_Controls.Logfilter->setEnabled(true); m_Controls.BModeMethod->setEnabled(true); m_Controls.buttonApplyBModeFilter->setEnabled(true); m_Controls.CutoffAbove->setEnabled(true); m_Controls.CutoffBelow->setEnabled(true); m_Controls.buttonApplyCropFilter->setEnabled(true); m_Controls.BPSpeedOfSound->setEnabled(true); m_Controls.buttonApplyBandpass->setEnabled(true); m_Controls.Partial->setEnabled(true); m_Controls.boundHigh->setEnabled(true); m_Controls.boundLow->setEnabled(true); m_Controls.BFAlgorithm->setEnabled(true); m_Controls.DelayCalculation->setEnabled(true); m_Controls.ImageType->setEnabled(true); m_Controls.Apodization->setEnabled(true); m_Controls.UseBP->setEnabled(true); #ifdef PHOTOACOUSTICS_USE_GPU m_Controls.UseGPUBf->setEnabled(true); m_Controls.UseGPUBmode->setEnabled(true); #endif m_Controls.BPhigh->setEnabled(true); m_Controls.BPlow->setEnabled(true); m_Controls.BPFalloffLow->setEnabled(true); m_Controls.BPFalloffHigh->setEnabled(true); m_Controls.UseImageSpacing->setEnabled(true); UseImageSpacing(); m_Controls.Pitch->setEnabled(true); m_Controls.ElementCount->setEnabled(true); m_Controls.SpeedOfSound->setEnabled(true); m_Controls.Samples->setEnabled(true); m_Controls.Lines->setEnabled(true); m_Controls.Angle->setEnabled(true); m_Controls.buttonApplyBeamforming->setEnabled(true); } void PAImageProcessing::DisableControls() { m_Controls.BatchProcessing->setEnabled(false); m_Controls.StepBeamforming->setEnabled(false); m_Controls.StepBandpass->setEnabled(false); m_Controls.StepCropping->setEnabled(false); m_Controls.StepBMode->setEnabled(false); m_Controls.SaveBeamforming->setEnabled(false); m_Controls.SaveBandpass->setEnabled(false); m_Controls.SaveCropping->setEnabled(false); m_Controls.SaveBMode->setEnabled(false); m_Controls.DoResampling->setEnabled(false); m_Controls.ResamplingValue->setEnabled(false); m_Controls.Logfilter->setEnabled(false); m_Controls.BModeMethod->setEnabled(false); m_Controls.buttonApplyBModeFilter->setEnabled(false); m_Controls.CutoffAbove->setEnabled(false); m_Controls.CutoffBelow->setEnabled(false); m_Controls.buttonApplyCropFilter->setEnabled(false); m_Controls.BPSpeedOfSound->setEnabled(false); m_Controls.buttonApplyBandpass->setEnabled(false); m_Controls.Partial->setEnabled(false); m_Controls.boundHigh->setEnabled(false); m_Controls.boundLow->setEnabled(false); m_Controls.BFAlgorithm->setEnabled(false); m_Controls.DelayCalculation->setEnabled(false); m_Controls.ImageType->setEnabled(false); m_Controls.Apodization->setEnabled(false); m_Controls.UseBP->setEnabled(false); #ifdef PHOTOACOUSTICS_USE_GPU m_Controls.UseGPUBf->setEnabled(false); m_Controls.UseGPUBmode->setEnabled(false); #endif m_Controls.BPhigh->setEnabled(false); m_Controls.BPlow->setEnabled(false); m_Controls.BPFalloffLow->setEnabled(false); m_Controls.BPFalloffHigh->setEnabled(false); m_Controls.UseImageSpacing->setEnabled(false); m_Controls.ScanDepth->setEnabled(false); m_Controls.Pitch->setEnabled(false); m_Controls.ElementCount->setEnabled(false); m_Controls.SpeedOfSound->setEnabled(false); m_Controls.Samples->setEnabled(false); m_Controls.Lines->setEnabled(false); m_Controls.Angle->setEnabled(false); m_Controls.buttonApplyBeamforming->setEnabled(false); } void PAImageProcessing::UseImageSpacing() { if (m_Controls.UseImageSpacing->isChecked()) { m_Controls.ScanDepth->setDisabled(true); } else { m_Controls.ScanDepth->setEnabled(true); } } #include void BeamformingThread::run() { mitk::Image::Pointer resultImage = mitk::Image::New(); mitk::Image::Pointer resultImageBuffer; std::function progressHandle = [this](int progress, std::string progressInfo) { emit updateProgress(progress, progressInfo); }; resultImageBuffer = m_FilterBank->ApplyBeamforming(m_InputImage, m_BFconfig, progressHandle); emit result(resultImageBuffer, "_bf"); } void BeamformingThread::setConfig(mitk::BeamformingSettings::Pointer BFconfig) { m_BFconfig = BFconfig; } void BeamformingThread::setInputImage(mitk::Image::Pointer image) { m_InputImage = image; } void BmodeThread::run() { mitk::Image::Pointer resultImage = m_FilterBank->ApplyBmodeFilter(m_InputImage, m_Method, m_UseLogfilter); if (m_ResampleSpacing != 0) { double desiredSpacing[2]{ m_InputImage->GetGeometry()->GetSpacing()[0], m_ResampleSpacing }; resultImage = m_FilterBank->ApplyResampling(resultImage, desiredSpacing); } emit result(resultImage, "_bmode"); } void BmodeThread::setConfig(bool useLogfilter, double resampleSpacing, mitk::PhotoacousticFilterService::BModeMethod method, bool useGPU) { m_UseLogfilter = useLogfilter; m_ResampleSpacing = resampleSpacing; m_Method = method; m_UseGPU = useGPU; } void BmodeThread::setInputImage(mitk::Image::Pointer image) { m_InputImage = image; } void CropThread::run() { mitk::Image::Pointer resultImage; - resultImage = m_FilterBank->ApplyCropping(m_InputImage, m_CutAbove, m_CutBelow, 0, 0, 0, 0); + resultImage = m_FilterBank->ApplyCropping(m_InputImage, m_CutAbove, m_CutBelow, 0, 0, m_CutSliceFirst, m_InputImage->GetDimension(2) - m_CutSliceLast); emit result(resultImage, "_cropped"); } void CropThread::setConfig(unsigned int CutAbove, unsigned int CutBelow, unsigned int CutSliceFirst, unsigned int CutSliceLast) { m_CutAbove = CutAbove; m_CutBelow = CutBelow; m_CutSliceLast = CutSliceLast; m_CutSliceFirst = CutSliceFirst; } void CropThread::setInputImage(mitk::Image::Pointer image) { m_InputImage = image; } void BandpassThread::run() { mitk::Image::Pointer resultImage = m_FilterBank->ApplyBandpassFilter(m_InputImage, m_BPHighPass, m_BPLowPass, m_TukeyAlphaHighPass, m_TukeyAlphaLowPass); emit result(resultImage, "_bandpassed"); } void BandpassThread::setConfig(float BPHighPass, float BPLowPass, float TukeyAlphaHighPass, float TukeyAlphaLowPass, float recordTime) { m_BPHighPass = BPHighPass; m_BPLowPass = BPLowPass; m_TukeyAlphaHighPass = TukeyAlphaHighPass; m_TukeyAlphaLowPass = TukeyAlphaLowPass; m_RecordTime = recordTime; } void BandpassThread::setInputImage(mitk::Image::Pointer image) { m_InputImage = image; } diff --git a/Plugins/org.mitk.gui.qt.photoacoustics.imageprocessing/src/internal/PAImageProcessingControls.ui b/Plugins/org.mitk.gui.qt.photoacoustics.imageprocessing/src/internal/PAImageProcessingControls.ui index 07818f493a..94a156f405 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,1001 +1,1005 @@ PAImageProcessingControls 0 0 385 1278 0 0 QmitkTemplate <html><head/><body><p><span style=" font-weight:600;">Batch Processing</span></p></body></html> Start Batch Processing Bandpass true Crop true Save false Save false Save true Beamform true BMode true Save true <html><head/><body><p><span style=" font-weight:600;">B-mode Filter Settings</span></p></body></html> - Envelope Detection + Absolute Filter Absolute Filter Envelope Detection Do Resampling true 0 0 13 0 11 3 0.010000000000000 1.000000000000000 0.010000000000000 0.150000000000000 [mm] Resampled Depth Spacing Logarithmic Compression Use GPU QLabel { color: rgb(255, 0, 0) } <html><head/><body><p align="center"><span style=" font-size:10pt; font-weight:600;">Please select an image!</span></p></body></html> 0 0 Do image processing Apply B-mode Filter Qt::Horizontal <html><head/><body><p><span style=" font-weight:600;">Bandpass Filter Settings</span></p></body></html> QLayout::SetDefaultConstraint 0 0 0 3 0.010000000000000 200.000000000000000 + + 0.100000000000000 + 1.000000000000000 [MHz] f High Pass [MHz] f Low Pass 0 0 3 200.000000000000000 + + 0.100000000000000 + 1 200.000000000000000 3000.000000000000000 5.000000000000000 1540.000000000000000 [m/s] Speed of Sound 1.000000000000000 0.100000000000000 0.500000000000000 Tukey Window α High Pass 2 1.000000000000000 0.100000000000000 0.000000000000000 Tukey Window α Low Pass <html><head/><body><p align="center"><span style=" font-size:10pt; font-weight:600; color:#ff0000;">Please select an image!</span></p></body></html> Apply Bandpass Qt::Horizontal <html><head/><body><p><span style=" font-weight:600;">Crop Filter Settings</span></p></body></html> 99999 300 Cut Top Cut Bottom 99999 5 800 + + + + false + + + 99999 + + + + + + + Maximal beamformed slice + + + max + + + + + + + false + + + 99999 + + + 10 + + + + + + + minimal beamformed slice + + + min + + + + + + + select slices + + + <html><head/><body><p align="center"><span style=" font-size:10pt; font-weight:600; color:#ff0000;">Please select an image!</span></p></body></html> Apply Crop Filer Qt::Horizontal <html><head/><body><p><span style=" font-weight:600;">Beamforming Filter Settings</span></p></body></html> 5 2 Delay Calculation Auto Get Depth true Apply Beamforming Beamforming Method [mm] Scan Depth 0 0 3 0.010000000000000 9.000000000000000 0.050000000000000 0.300000000000000 Transducer Elements 0 0 4 300.000000000000000 0.100000000000000 50.000000000000000 [mm] Transducer Pitch 0 0 64 1024 128 128 0 0 256 16384 256 2048 0 0 64 2048 128 256 Samples Reconstruction Lines true 0 0 100 0 0 0 DAS DMAS sDMAS 0 0 - Spherical Wave + Quad. Approx. Quad. Approx. Spherical Wave 0 0 PA Image US Image Image Type 0 0 Von Hann Hamming Box Apodization 0 0 1 200.000000000000000 3000.000000000000000 5.000000000000000 1540.000000000000000 [m/s] Speed of Sound - - - - - false - - - 99999 - - - - - - - minimal beamformed slice - - - min - - - - + - - - - - false - - - 99999 - - - 10 - - - - - - - Maximal beamformed slice - - - max - - - - + - - - - select slices - - - Compute On GPU true true Auto Use Bandpass 0 0 1 1.000000000000000 180.000000000000000 27.000000000000000 [°] Element Angle <html><head/><body><p align="center"><span style=" font-size:10pt; font-weight:600; color:#ff0000;">Please select an image!</span></p></body></html> Qt::Vertical 20 40