diff --git a/Modules/OpenCL/mitkOclDataSetToDataSetFilter.cpp b/Modules/OpenCL/mitkOclDataSetToDataSetFilter.cpp index be0075557e..1a7e821ab1 100644 --- a/Modules/OpenCL/mitkOclDataSetToDataSetFilter.cpp +++ b/Modules/OpenCL/mitkOclDataSetToDataSetFilter.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. ===================================================================*/ #include "mitkOclDataSetToDataSetFilter.h" #include "mitkOclDataSet.h" #include "mitkException.h" mitk::OclDataSetToDataSetFilter::OclDataSetToDataSetFilter() { m_Output = mitk::OclDataSet::New(); } mitk::OclDataSetToDataSetFilter::~OclDataSetToDataSetFilter() { } mitk::OclDataSet::Pointer mitk::OclDataSetToDataSetFilter::GetGPUOutput() { return this->m_Output; } void* mitk::OclDataSetToDataSetFilter::GetOutput() { void* pData = m_Output->TransferDataToCPU(m_CommandQue); return pData; } int mitk::OclDataSetToDataSetFilter::GetBytesPerElem() { return this->m_CurrentSizeOutput; } bool mitk::OclDataSetToDataSetFilter::InitExec(cl_kernel ckKernel, unsigned int* dimensions, size_t outputDataSize, unsigned int outputBpE) { cl_int clErr = 0; if (m_Input.IsNull()) mitkThrow() << "Input DataSet is null."; // get DataSet size once const unsigned int uiDataSetWidth = dimensions[0]; const unsigned int uiDataSetHeight = dimensions[1]; const unsigned int uiDataSetDepth = dimensions[2]; // compute work sizes this->SetWorkingSize(8, uiDataSetWidth, 8, uiDataSetHeight, 8, uiDataSetDepth); cl_mem clBuffIn = m_Input->GetGPUBuffer(); cl_mem clBuffOut = m_Output->GetGPUBuffer(); if (!clBuffIn) { if (m_Input->TransferDataToGPU(m_CommandQue) != CL_SUCCESS) { mitkThrow() << "Could not create / initialize gpu DataSet."; } clBuffIn = m_Input->GetGPUBuffer(); } // output DataSet not initialized or output buffer size changed if (!clBuffOut || (size_t)m_Output->GetBufferSize() != outputDataSize) { MITK_DEBUG << "Create GPU DataSet call " << uiDataSetWidth << "x" << uiDataSetHeight << "x" << uiDataSetDepth; + MITK_INFO << "Create GPU Buffer of size " << outputDataSize * outputBpE / 1024.f / 1024.f << "MB"; m_Output->SetBpE(outputBpE); m_Output->SetBufferSize(outputDataSize); clBuffOut = m_Output->CreateGPUBuffer(); m_CurrentSizeOutput = outputBpE; } clErr = 0; clErr = clSetKernelArg(ckKernel, 0, sizeof(cl_mem), &clBuffIn); clErr |= clSetKernelArg(ckKernel, 1, sizeof(cl_mem), &clBuffOut); CHECK_OCL_ERR(clErr); if (clErr != CL_SUCCESS) mitkThrow() << "OpenCL Part initialization failed with " << GetOclErrorAsString(clErr); return(clErr == CL_SUCCESS); } bool mitk::OclDataSetToDataSetFilter::InitExecNoInput(cl_kernel ckKernel, unsigned int* dimensions, size_t outputDataSize, unsigned int outputBpE) { cl_int clErr = 0; // get DataSet size once const unsigned int uiDataSetWidth = dimensions[0]; const unsigned int uiDataSetHeight = dimensions[1]; const unsigned int uiDataSetDepth = dimensions[2]; // compute work sizes this->SetWorkingSize(8, uiDataSetWidth, 8, uiDataSetHeight, 8, uiDataSetDepth); cl_mem clBuffOut = m_Output->GetGPUBuffer(); // output DataSet not initialized or output buffer size changed if (!clBuffOut || (size_t)m_Output->GetBufferSize() != outputDataSize) { MITK_DEBUG << "Create GPU DataSet call " << uiDataSetWidth << "x" << uiDataSetHeight << "x" << uiDataSetDepth; m_Output->SetBpE(outputBpE); m_Output->SetBufferSize(outputDataSize); clBuffOut = m_Output->CreateGPUBuffer(); m_CurrentSizeOutput = outputBpE; } clErr = clSetKernelArg(ckKernel, 0, sizeof(cl_mem), &clBuffOut); CHECK_OCL_ERR(clErr); if (clErr != CL_SUCCESS) mitkThrow() << "OpenCL Part initialization failed with " << GetOclErrorAsString(clErr); return(clErr == CL_SUCCESS); } diff --git a/Modules/PhotoacousticsAlgorithms/include/OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.h b/Modules/PhotoacousticsAlgorithms/include/OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.h index 7490c71f57..b60cd2315d 100644 --- a/Modules/PhotoacousticsAlgorithms/include/OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.h +++ b/Modules/PhotoacousticsAlgorithms/include/OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.h @@ -1,145 +1,145 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef _MITKPHOTOACOUSTICSOCLBEAMFORMER_H_ #define _MITKPHOTOACOUSTICSOCLBEAMFORMER_H_ #include #if defined(PHOTOACOUSTICS_USE_GPU) || DOXYGEN #include "mitkOclDataSetToDataSetFilter.h" #include "mitkBeamformingSettings.h" #include "mitkPhotoacousticOCLDelayCalculation.h" #include "mitkPhotoacousticOCLUsedLinesCalculation.h" #include namespace mitk { /*! * \brief Class implementing a mitk::OclDataSetToDataSetFilter for beamforming on GPU * * The class must be given a configuration class instance of mitk::BeamformingSettings for beamforming parameters through mitk::PhotoacousticOCLBeamformingFilter::SetConfig(BeamformingSettings settings) * Additional configuration of the apodisation function is needed. */ class PhotoacousticOCLBeamformingFilter : public OclDataSetToDataSetFilter, public itk::Object { public: mitkClassMacroItkParent(PhotoacousticOCLBeamformingFilter, itk::Object); mitkNewMacro1Param(Self, BeamformingSettings::Pointer); /** * @brief SetInput Set the input data through an image. Arbitrary images are supported */ - void SetInput(Image::Pointer image, unsigned int inputSlices); + 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_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/source/OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.cpp b/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.cpp index 3c3759fe24..9bf7ace7cc 100644 --- a/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.cpp +++ b/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.cpp @@ -1,259 +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_inputSlices(1), m_Conf(settings), m_InputImage(mitk::Image::New()), m_ApodizationBuffer(nullptr), m_MemoryLocationsBuffer(nullptr), m_DelaysBuffer(nullptr), m_UsedLinesBuffer(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); } 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_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; } //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), &(m_OutputDim[2])); + 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)); // 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)) + 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, 16, 50)) + 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()) { 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); return (OclFilter::IsInitialized() && buildErr); } -void mitk::PhotoacousticOCLBeamformingFilter::SetInput(mitk::Image::Pointer image, unsigned int inputSlices) +void mitk::PhotoacousticOCLBeamformingFilter::SetInput(mitk::Image::Pointer image) { OclDataSetToDataSetFilter::SetInput(image); m_InputImage = image; - m_inputSlices = inputSlices; + 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/filters/mitkBeamformingFilter.cpp b/Modules/PhotoacousticsAlgorithms/source/filters/mitkBeamformingFilter.cpp index 203d1bd952..24a92884df 100644 --- a/Modules/PhotoacousticsAlgorithms/source/filters/mitkBeamformingFilter.cpp +++ b/Modules/PhotoacousticsAlgorithms/source/filters/mitkBeamformingFilter.cpp @@ -1,298 +1,298 @@ /*=================================================================== 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) {}; #if defined(PHOTOACOUSTICS_USE_GPU) m_BeamformingOclFilter = mitk::PhotoacousticOCLBeamformingFilter::New(m_Conf); #else m_BeamformingOclFilter = mitk::PhotoacousticOCLBeamformingFilter::New(); #endif 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; mitk::Vector3D spacing; spacing[0] = m_Conf->GetPitchInMeters() * m_Conf->GetTransducerElements() * 1000 / m_Conf->GetReconstructionLines(); float desiredYSpacing = m_Conf->GetReconstructionDepth() * 1000 / m_Conf->GetSamplesPerLine(); float maxYSpacing = m_Conf->GetSpeedOfSound() * m_Conf->GetTimeSpacing() * input->GetDimension(1) / m_Conf->GetSamplesPerLine() * 1000; spacing[1] = desiredYSpacing < maxYSpacing ? desiredYSpacing : maxYSpacing; spacing[2] = 1; unsigned int dim[] = { m_Conf->GetReconstructionLines(), m_Conf->GetSamplesPerLine(), input->GetDimension(2)}; output->Initialize(mitk::MakeScalarPixelType(), 3, dim); 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 = 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 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 { 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(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, num_Slices); + m_BeamformingOclFilter->SetInput(inputBatch); 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); } free(out); } } 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; }