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/mitkBandpassFilter.cpp b/Modules/PhotoacousticsAlgorithms/source/filters/mitkBandpassFilter.cpp index d9c711177f..438aa798d7 100644 --- a/Modules/PhotoacousticsAlgorithms/source/filters/mitkBandpassFilter.cpp +++ b/Modules/PhotoacousticsAlgorithms/source/filters/mitkBandpassFilter.cpp @@ -1,267 +1,264 @@ /*=================================================================== mitkCastToFloatImageFilter The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #define _USE_MATH_DEFINES #include #include "mitkBandpassFilter.h" #include "../ITKFilter/ITKUltrasound/itkFFT1DComplexConjugateToRealImageFilter.h" #include "../ITKFilter/ITKUltrasound/itkFFT1DRealToComplexConjugateImageFilter.h" #include "itkComplexToModulusImageFilter.h" #include "itkMultiplyImageFilter.h" #include "mitkIOUtil.h" #include "mitkITKImageImport.h" #include "mitkImageCast.h" mitk::BandpassFilter::BandpassFilter() : m_HighPass(0), m_LowPass(50), m_HighPassAlpha(1), m_LowPassAlpha(1), m_FilterService(mitk::PhotoacousticFilterService::New()) { MITK_INFO << "Instantiating BandpassFilter..."; SetNumberOfIndexedInputs(1); SetNumberOfIndexedOutputs(1); MITK_INFO << "Instantiating BandpassFilter...[Done]"; } mitk::BandpassFilter::~BandpassFilter() { MITK_INFO << "Destructed BandpassFilter."; } void mitk::BandpassFilter::SanityCheckPreconditions() { auto input = GetInput(); std::string type = input->GetPixelType().GetTypeAsString(); if (!(type == "scalar (float)" || type == " (float)")) { MITK_ERROR << "This filter can currently only handle float type images."; mitkThrow() << "This filter can currently only handle float type images."; } if (m_HighPass > m_LowPass) { MITK_ERROR << "High Pass is higher than Low Pass; aborting."; mitkThrow() << "High Pass is higher than Low Pass; aborting."; } } itk::Image::Pointer BPFunction(mitk::Image::Pointer reference, float cutoffFrequencyPixelHighPass, float cutoffFrequencyPixelLowPass, float alphaHighPass, float alphaLowPass) { float *imageData = new float[reference->GetDimension(0) * reference->GetDimension(1)]; float width = cutoffFrequencyPixelLowPass - cutoffFrequencyPixelHighPass; float center = cutoffFrequencyPixelHighPass + width / 2.f; for (unsigned int n = 0; n < reference->GetDimension(1); ++n) { imageData[reference->GetDimension(0) * n] = 0; } for (int n = 0; n < width; ++n) { imageData[reference->GetDimension(0) * (int)(n + center - (width / 2.f))] = 1; if (n <= (alphaHighPass * (width - 1)) / 2.f) { if (alphaHighPass > 0.00001f) { imageData[reference->GetDimension(0) * (int)(n + center - (width / 2.f))] = (1 + cos(itk::Math::pi * (2 * n / (alphaHighPass * (width - 1)) - 1))) / 2; } else { imageData[reference->GetDimension(0) * (int)(n + center - (width / 2.f))] = 1; } } else if (n >= (width - 1) * (1 - alphaLowPass / 2.f)) { if (alphaLowPass > 0.00001f) { imageData[reference->GetDimension(0) * (int)(n + center - (width / 2.f))] = (1 + cos(itk::Math::pi * (2 * n / (alphaLowPass * (width - 1)) + 1 - 2 / alphaLowPass))) / 2; } else { imageData[reference->GetDimension(0) * (int)(n + center - (width / 2.f))] = 1; } } } for (unsigned int n = reference->GetDimension(1) / 2; n < reference->GetDimension(1); ++n) { imageData[reference->GetDimension(0) * n] = imageData[(reference->GetDimension(1) - (n + 1)) * reference->GetDimension(0)]; } for (unsigned int line = 1; line < reference->GetDimension(0); ++line) { for (unsigned int sample = 0; sample < reference->GetDimension(1); ++sample) { imageData[reference->GetDimension(0) * sample + line] = imageData[reference->GetDimension(0) * sample]; } } typedef itk::Image ImageType; ImageType::RegionType region; ImageType::IndexType start; start.Fill(0); region.SetIndex(start); ImageType::SizeType size; size[0] = reference->GetDimension(0); size[1] = reference->GetDimension(1); size[2] = reference->GetDimension(2); region.SetSize(size); ImageType::SpacingType SpacingItk; SpacingItk[0] = reference->GetGeometry()->GetSpacing()[0]; SpacingItk[1] = reference->GetGeometry()->GetSpacing()[1]; SpacingItk[2] = reference->GetGeometry()->GetSpacing()[2]; ImageType::Pointer image = ImageType::New(); image->SetRegions(region); image->Allocate(); image->FillBuffer(itk::NumericTraits::Zero); image->SetSpacing(SpacingItk); ImageType::IndexType pixelIndex; for (unsigned int slice = 0; slice < reference->GetDimension(2); ++slice) { for (unsigned int line = 0; line < reference->GetDimension(0); ++line) { for (unsigned int sample = 0; sample < reference->GetDimension(1); ++sample) { pixelIndex[0] = line; pixelIndex[1] = sample; pixelIndex[2] = slice; image->SetPixel(pixelIndex, imageData[line + sample * reference->GetDimension(0)]); } } } delete[] imageData; return image; } void mitk::BandpassFilter::GenerateData() { SanityCheckPreconditions(); auto input = GetInput(); auto output = GetOutput(); mitk::Image::Pointer resampledInput = input; double powerOfTwo = std::log2(input->GetDimension(1)); int finalSize = 0; double spacingResize = 1; // check if this is a power of two by checking that log2 is int if (std::fmod(powerOfTwo, 1.0) >= std::numeric_limits::epsilon()) { finalSize = (int)pow(2, std::ceil(powerOfTwo)); double dim[2] = {(double)input->GetDimension(0), (double)finalSize }; resampledInput = m_FilterService->ApplyResamplingToDim(input, dim); spacingResize = (double)input->GetDimension(1) / finalSize; } // do a fourier transform, multiply with an appropriate window for the filter, and transform back typedef itk::Image RealImageType; - MITK_INFO << resampledInput->GetDimension(0) << " " << resampledInput->GetDimension(1) << " " << resampledInput->GetDimension(2) << "asd" << sizeof(long double); RealImageType::Pointer image; mitk::CastToItkImage(resampledInput, image); typedef itk::FFT1DRealToComplexConjugateImageFilter ForwardFFTFilterType; typedef ForwardFFTFilterType::OutputImageType ComplexImageType; ForwardFFTFilterType::Pointer forwardFFTFilter = ForwardFFTFilterType::New(); forwardFFTFilter->SetInput(image); forwardFFTFilter->SetDirection(1); try { forwardFFTFilter->UpdateOutputInformation(); } catch (itk::ExceptionObject &error) { std::cerr << "Error: " << error << std::endl; MITK_ERROR << "Bandpass could not be applied"; mitkThrow() << "bandpass error"; } if (m_HighPass > m_LowPass) mitkThrow() << "High pass frequency higher than low pass frequency, abort"; float singleVoxel = spacingResize / (m_TimeSpacing * resampledInput->GetDimension(1)); // [Hz] if(m_IsBFImage) singleVoxel = spacingResize / (resampledInput->GetGeometry()->GetSpacing()[1] / 1e3 / m_SpeedOfSound * resampledInput->GetDimension(1)); // [Hz] float cutoffPixelHighPass = std::min((m_HighPass / singleVoxel), (float)resampledInput->GetDimension(1) / 2.0f); float cutoffPixelLowPass = std::min((m_LowPass / singleVoxel), (float)resampledInput->GetDimension(1) / 2.0f); MITK_INFO << "SingleVoxel: " << singleVoxel; MITK_INFO << "HighPass: " << m_HighPass; MITK_INFO << "LowPass: " << m_LowPass; MITK_INFO << "cutoffPixelHighPass: " << cutoffPixelHighPass; MITK_INFO << "cutoffPixelLowPass: " << cutoffPixelLowPass; RealImageType::Pointer fftMultiplicator = BPFunction(resampledInput, cutoffPixelHighPass, cutoffPixelLowPass, m_HighPassAlpha, m_LowPassAlpha); typedef itk::MultiplyImageFilter MultiplyFilterType; MultiplyFilterType::Pointer multiplyFilter = MultiplyFilterType::New(); multiplyFilter->SetInput1(forwardFFTFilter->GetOutput()); multiplyFilter->SetInput2(fftMultiplicator); /* GrabItkImageMemory(fftMultiplicator, output); mitk::IOUtil::Save(output, "D:/fftImage.nrrd"); typedef itk::ComplexToModulusImageFilter< ComplexImageType, RealImageType > modulusType; modulusType::Pointer modul = modulusType::New(); modul->SetInput(multiplyFilter->GetOutput()); GrabItkImageMemory(modul->GetOutput(), output); mitk::IOUtil::Save(output, "D:/fftout.nrrd"); modul->SetInput(forwardFFTFilter->GetOutput()); GrabItkImageMemory(modul->GetOutput(), output); mitk::IOUtil::Save(output, "D:/fftin.nrrd");*/ typedef itk::FFT1DComplexConjugateToRealImageFilter InverseFilterType; InverseFilterType::Pointer inverseFFTFilter = InverseFilterType::New(); inverseFFTFilter->SetInput(multiplyFilter->GetOutput()); inverseFFTFilter->SetDirection(1); GrabItkImageMemory(inverseFFTFilter->GetOutput(), output); - double spacing_x = input->GetGeometry()->GetSpacing()[0]; - double spacing_y = input->GetGeometry()->GetSpacing()[1]; - double dim[2] = {spacing_x, spacing_y}; - auto resampledOutput = m_FilterService->ApplyResampling(output, dim); + double dim[2] = { (double)input->GetDimension(0), (double)input->GetDimension(1) }; + auto resampledOutput = m_FilterService->ApplyResamplingToDim(output, dim); output->Initialize(mitk::MakeScalarPixelType(), 3, input->GetDimensions()); output->SetSpacing(resampledOutput->GetGeometry()->GetSpacing()); ImageReadAccessor copy(resampledOutput); output->SetImportVolume(copy.GetData()); } 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; }