diff --git a/Modules/PhotoacousticsAlgorithms/Resources/DelayCalculation.cl b/Modules/PhotoacousticsAlgorithms/Resources/DelayCalculation.cl index 0be471b2b4..acb4aafb9c 100644 --- a/Modules/PhotoacousticsAlgorithms/Resources/DelayCalculation.cl +++ b/Modules/PhotoacousticsAlgorithms/Resources/DelayCalculation.cl @@ -1,69 +1,69 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ __kernel void ckDelayCalculationQuad( __global unsigned short *gDest, __global unsigned short *usedLines, unsigned int inputL, unsigned int inputS, unsigned int outputL, unsigned int outputS, char isPAImage, float delayMultiplicatorRaw // parameters ) { uint globalPosX = get_global_id(0); uint globalPosY = get_global_id(1); if (globalPosX * 2 < outputL && globalPosY < outputS) { float l_i = 0; // we calculate the delays relative to line zero - float s_i = (float)globalPosY / (float)outputS * (float)inputS / 2; + float s_i = (float)globalPosY / (float)outputS * (float)inputS / (2-isPAImage); float l_s = (float)globalPosX / (float)outputL * (float)inputL; // the currently calculated line float delayMultiplicator = delayMultiplicatorRaw / s_i; gDest[globalPosY * (outputL / 2) + globalPosX] = delayMultiplicator * pow((l_s - l_i), 2) + s_i + (1-isPAImage)*s_i; } } __kernel void ckDelayCalculationSphe( __global unsigned short *gDest, __global unsigned short *usedLines, unsigned int inputL, unsigned int inputS, unsigned int outputL, unsigned int outputS, char isPAImage, float delayMultiplicatorRaw // parameters ) { uint globalPosX = get_global_id(0); uint globalPosY = get_global_id(1); if (globalPosX * 2 < outputL && globalPosY < outputS) { float l_i = 0; // we calculate the delays relative to line zero - float s_i = (float)globalPosY / (float)outputS * (float)inputS / 2; + float s_i = (float)globalPosY / (float)outputS * (float)inputS / (2-isPAImage); float l_s = (float)globalPosX / (float)outputL * (float)inputL; // the currently calculated line gDest[globalPosY * (outputL / 2) + globalPosX] = sqrt( pow(s_i, 2) + pow((delayMultiplicatorRaw * ((l_s - l_i)) / inputL), 2) ) + (1-isPAImage)*s_i; } } \ No newline at end of file diff --git a/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.cpp b/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.cpp index 588047ab63..16960e2e47 100644 --- a/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.cpp +++ b/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.cpp @@ -1,261 +1,260 @@ /*=================================================================== 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_Conf(settings) { - MITK_INFO << "Initializing OCL beamforming Filter..."; + 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 << "Initializing OCL beamforming Filter...[Done]"; + 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 { - MITK_INFO << "Updating Workgroup size for new dimensions"; size_t outputSize = (size_t)m_Conf->GetReconstructionLines() * (size_t)m_Conf->GetSamplesPerLine() * (size_t)m_Conf->GetInputDim()[2]; m_OutputDim[0] = m_Conf->GetReconstructionLines(); m_OutputDim[1] = m_Conf->GetSamplesPerLine(); m_OutputDim[2] = m_Conf->GetInputDim()[2]; 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, 9, sizeof(cl_uint), &(reconstructionLines)); clErr |= clSetKernelArg(this->m_PixelCalculation, 10, sizeof(cl_uint), &(samplesPerLine)); // execute the filter on a 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) { OclDataSetToDataSetFilter::SetInput(image); m_InputImage = image; } 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 31b30a8f50..ad4ee3b1a8 100644 --- a/Modules/PhotoacousticsAlgorithms/source/filters/mitkBeamformingFilter.cpp +++ b/Modules/PhotoacousticsAlgorithms/source/filters/mitkBeamformingFilter.cpp @@ -1,298 +1,299 @@ /*=================================================================== 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 * m_Conf->GetSpeedOfSound() * 1000 / m_Conf->GetSamplesPerLine(); + 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 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 * ((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) { } if (batchSize < 1) { MITK_ERROR << "device memory too small for GPU beamforming"; return; } mitk::ImageReadAccessor copy(input); for (unsigned int i = 0; i < batches; ++i) { m_ProgressHandle(input->GetDimension(2) / batches * i, "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->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/utils/mitkBeamformingUtils.cpp b/Modules/PhotoacousticsAlgorithms/source/utils/mitkBeamformingUtils.cpp index f805227198..86dd243c0f 100644 --- a/Modules/PhotoacousticsAlgorithms/source/utils/mitkBeamformingUtils.cpp +++ b/Modules/PhotoacousticsAlgorithms/source/utils/mitkBeamformingUtils.cpp @@ -1,548 +1,548 @@ /*=================================================================== mitkPhotoacousticBeamformingFilter The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkProperties.h" #include "mitkImageReadAccessor.h" #include #include #include #include #include #include "mitkImageCast.h" #include "mitkBeamformingUtils.h" mitk::BeamformingUtils::BeamformingUtils() { } mitk::BeamformingUtils::~BeamformingUtils() { } float* mitk::BeamformingUtils::VonHannFunction(int samples) { float* ApodWindow = new float[samples]; for (int n = 0; n < samples; ++n) { ApodWindow[n] = (1 - cos(2 * itk::Math::pi * n / (samples - 1))) / 2; } return ApodWindow; } float* mitk::BeamformingUtils::HammFunction(int samples) { float* ApodWindow = new float[samples]; for (int n = 0; n < samples; ++n) { ApodWindow[n] = 0.54 - 0.46*cos(2 * itk::Math::pi*n / (samples - 1)); } return ApodWindow; } float* mitk::BeamformingUtils::BoxFunction(int samples) { float* ApodWindow = new float[samples]; for (int n = 0; n < samples; ++n) { ApodWindow[n] = 1; } return ApodWindow; } void mitk::BeamformingUtils::DASQuadraticLine( float* input, float* output, float inputDim[2], float outputDim[2], const short& line, const mitk::BeamformingSettings::Pointer config) { const float* apodisation = config->GetApodizationFunction(); const short apodArraySize = config->GetApodizationArraySize(); float& inputS = inputDim[1]; float& inputL = inputDim[0]; float& outputS = outputDim[1]; float& outputL = outputDim[0]; short AddSample = 0; short maxLine = 0; short minLine = 0; float delayMultiplicator = 0; float l_i = 0; float s_i = 0; float part = 0; float tan_phi = std::tan(config->GetAngle() / 360 * 2 * itk::Math::pi); float part_multiplicator = tan_phi * config->GetTimeSpacing() * config->GetSpeedOfSound() / config->GetPitchInMeters() * inputL / config->GetTransducerElements(); float apod_mult = 1; short usedLines = (maxLine - minLine); //quadratic delay l_i = line / outputL * inputL; for (short sample = 0; sample < outputS; ++sample) { - s_i = (float)sample / outputS * inputS / 2; + s_i = (float)sample / outputS * inputS / (2-(int)config->GetIsPhotoacousticImage()); part = part_multiplicator*s_i; if (part < 1) part = 1; maxLine = (short)std::min((l_i + part) + 1, inputL); minLine = (short)std::max((l_i - part), 0.0f); usedLines = (maxLine - minLine); apod_mult = (float)apodArraySize / (float)usedLines; delayMultiplicator = pow((1 / (config->GetTimeSpacing()*config->GetSpeedOfSound()) * (config->GetPitchInMeters()*config->GetTransducerElements()) / inputL), 2) / s_i / 2; for (short l_s = minLine; l_s < maxLine; ++l_s) { AddSample = delayMultiplicator * pow((l_s - l_i), 2) + s_i + (1 - config->GetIsPhotoacousticImage())*s_i; if (AddSample < inputS && AddSample >= 0) output[sample*(short)outputL + line] += input[l_s + AddSample*(short)inputL] * apodisation[(short)((l_s - minLine)*apod_mult)]; else --usedLines; } output[sample*(short)outputL + line] = output[sample*(short)outputL + line] / usedLines; } } void mitk::BeamformingUtils::DASSphericalLine( float* input, float* output, float inputDim[2], float outputDim[2], const short& line, const mitk::BeamformingSettings::Pointer config) { const float* apodisation = config->GetApodizationFunction(); const short apodArraySize = config->GetApodizationArraySize(); float& inputS = inputDim[1]; float& inputL = inputDim[0]; float& outputS = outputDim[1]; float& outputL = outputDim[0]; short AddSample = 0; short maxLine = 0; short minLine = 0; float l_i = 0; float s_i = 0; float part = 0.07 * inputL; float tan_phi = std::tan(config->GetAngle() / 360 * 2 * itk::Math::pi); float part_multiplicator = tan_phi * config->GetTimeSpacing() * config->GetSpeedOfSound() / config->GetPitchInMeters() * inputL / (float)config->GetTransducerElements(); float apod_mult = 1; short usedLines = (maxLine - minLine); //exact delay l_i = (float)line / outputL * inputL; for (short sample = 0; sample < outputS; ++sample) { - s_i = (float)sample / outputS * inputS / 2; + s_i = (float)sample / outputS * inputS / (2 - (int)config->GetIsPhotoacousticImage()); part = part_multiplicator*s_i; if (part < 1) part = 1; maxLine = (short)std::min((l_i + part) + 1, inputL); minLine = (short)std::max((l_i - part), 0.0f); usedLines = (maxLine - minLine); apod_mult = (float)apodArraySize / (float)usedLines; for (short l_s = minLine; l_s < maxLine; ++l_s) { AddSample = (int)sqrt( pow(s_i, 2) + pow((1 / (config->GetTimeSpacing()*config->GetSpeedOfSound()) * (((float)l_s - l_i)*config->GetPitchInMeters()*(float)config->GetTransducerElements()) / inputL), 2) ) + (1 - config->GetIsPhotoacousticImage())*s_i; if (AddSample < inputS && AddSample >= 0) output[sample*(short)outputL + line] += input[l_s + AddSample*(short)inputL] * apodisation[(short)((l_s - minLine)*apod_mult)]; else --usedLines; } output[sample*(short)outputL + line] = output[sample*(short)outputL + line] / usedLines; } } void mitk::BeamformingUtils::DMASQuadraticLine( float* input, float* output, float inputDim[2], float outputDim[2], const short& line, const mitk::BeamformingSettings::Pointer config) { const float* apodisation = config->GetApodizationFunction(); const short apodArraySize = config->GetApodizationArraySize(); float& inputS = inputDim[1]; float& inputL = inputDim[0]; float& outputS = outputDim[1]; float& outputL = outputDim[0]; short maxLine = 0; short minLine = 0; float delayMultiplicator = 0; float l_i = 0; float s_i = 0; float part = 0.07 * inputL; float tan_phi = std::tan(config->GetAngle() / 360 * 2 * itk::Math::pi); float part_multiplicator = tan_phi * config->GetTimeSpacing() * config->GetSpeedOfSound() / config->GetPitchInMeters() * inputL / (float)config->GetTransducerElements(); float apod_mult = 1; float mult = 0; short usedLines = (maxLine - minLine); //quadratic delay l_i = line / outputL * inputL; for (short sample = 0; sample < outputS; ++sample) { - s_i = sample / outputS * inputS / 2; + s_i = sample / outputS * inputS / (2 - (int)config->GetIsPhotoacousticImage()); part = part_multiplicator*s_i; if (part < 1) part = 1; maxLine = (short)std::min((l_i + part) + 1, inputL); minLine = (short)std::max((l_i - part), 0.0f); usedLines = (maxLine - minLine); apod_mult = (float)apodArraySize / (float)usedLines; delayMultiplicator = pow((1 / (config->GetTimeSpacing()*config->GetSpeedOfSound()) * (config->GetPitchInMeters()*config->GetTransducerElements()) / inputL), 2) / s_i / 2; //calculate the AddSamples beforehand to save some time short* AddSample = new short[maxLine - minLine]; for (short l_s = 0; l_s < maxLine - minLine; ++l_s) { AddSample[l_s] = (short)(delayMultiplicator * pow((minLine + l_s - l_i), 2) + s_i) + (1 - config->GetIsPhotoacousticImage())*s_i; } float s_1 = 0; float s_2 = 0; for (short l_s1 = minLine; l_s1 < maxLine - 1; ++l_s1) { if (AddSample[l_s1 - minLine] < inputS && AddSample[l_s1 - minLine] >= 0) { for (short l_s2 = l_s1 + 1; l_s2 < maxLine; ++l_s2) { if (AddSample[l_s2 - minLine] < inputS && AddSample[l_s2 - minLine] >= 0) { s_2 = input[l_s2 + AddSample[l_s2 - minLine] * (short)inputL]; s_1 = input[l_s1 + AddSample[l_s1 - minLine] * (short)inputL]; mult = s_2 * apodisation[(int)((l_s2 - minLine)*apod_mult)] * s_1 * apodisation[(int)((l_s1 - minLine)*apod_mult)]; output[sample*(short)outputL + line] += sqrt(fabs(mult)) * ((mult > 0) - (mult < 0)); } } } else --usedLines; } output[sample*(short)outputL + line] = output[sample*(short)outputL + line] / (float)(pow(usedLines, 2) - (usedLines - 1)); delete[] AddSample; } } void mitk::BeamformingUtils::DMASSphericalLine( float* input, float* output, float inputDim[2], float outputDim[2], const short& line, const mitk::BeamformingSettings::Pointer config) { const float* apodisation = config->GetApodizationFunction(); const short apodArraySize = config->GetApodizationArraySize(); float& inputS = inputDim[1]; float& inputL = inputDim[0]; float& outputS = outputDim[1]; float& outputL = outputDim[0]; short maxLine = 0; short minLine = 0; float l_i = 0; float s_i = 0; float part = 0.07 * inputL; float tan_phi = std::tan(config->GetAngle() / 360 * 2 * itk::Math::pi); float part_multiplicator = tan_phi * config->GetTimeSpacing() * config->GetSpeedOfSound() / config->GetPitchInMeters() * inputL / (float)config->GetTransducerElements(); float apod_mult = 1; float mult = 0; short usedLines = (maxLine - minLine); //exact delay l_i = (float)line / outputL * inputL; for (short sample = 0; sample < outputS; ++sample) { - s_i = (float)sample / outputS * inputS / 2; + s_i = (float)sample / outputS * inputS / (2 - (int)config->GetIsPhotoacousticImage()); part = part_multiplicator*s_i; if (part < 1) part = 1; maxLine = (short)std::min((l_i + part) + 1, inputL); minLine = (short)std::max((l_i - part), 0.0f); usedLines = (maxLine - minLine); apod_mult = (float)apodArraySize / (float)usedLines; //calculate the AddSamples beforehand to save some time short* AddSample = new short[maxLine - minLine]; for (short l_s = 0; l_s < maxLine - minLine; ++l_s) { AddSample[l_s] = (short)sqrt( pow(s_i, 2) + pow((1 / (config->GetTimeSpacing()*config->GetSpeedOfSound()) * (((float)minLine + (float)l_s - l_i)*config->GetPitchInMeters()*(float)config->GetTransducerElements()) / inputL), 2) ) + (1 - config->GetIsPhotoacousticImage())*s_i; } float s_1 = 0; float s_2 = 0; for (short l_s1 = minLine; l_s1 < maxLine - 1; ++l_s1) { if (AddSample[l_s1 - minLine] < inputS && AddSample[l_s1 - minLine] >= 0) { for (short l_s2 = l_s1 + 1; l_s2 < maxLine; ++l_s2) { if (AddSample[l_s2 - minLine] < inputS && AddSample[l_s2 - minLine] >= 0) { s_2 = input[l_s2 + AddSample[l_s2 - minLine] * (short)inputL]; s_1 = input[l_s1 + AddSample[l_s1 - minLine] * (short)inputL]; mult = s_2 * apodisation[(int)((l_s2 - minLine)*apod_mult)] * s_1 * apodisation[(int)((l_s1 - minLine)*apod_mult)]; output[sample*(short)outputL + line] += sqrt(fabs(mult)) * ((mult > 0) - (mult < 0)); } } } else --usedLines; } output[sample*(short)outputL + line] = output[sample*(short)outputL + line] / (float)(pow(usedLines, 2) - (usedLines - 1)); delete[] AddSample; } } void mitk::BeamformingUtils::sDMASQuadraticLine( float* input, float* output, float inputDim[2], float outputDim[2], const short& line, const mitk::BeamformingSettings::Pointer config) { const float* apodisation = config->GetApodizationFunction(); const short apodArraySize = config->GetApodizationArraySize(); float& inputS = inputDim[1]; float& inputL = inputDim[0]; float& outputS = outputDim[1]; float& outputL = outputDim[0]; short maxLine = 0; short minLine = 0; float delayMultiplicator = 0; float l_i = 0; float s_i = 0; float part = 0.07 * inputL; float tan_phi = std::tan(config->GetAngle() / 360 * 2 * itk::Math::pi); float part_multiplicator = tan_phi * config->GetTimeSpacing() * config->GetSpeedOfSound() / config->GetPitchInMeters() * inputL / (float)config->GetTransducerElements(); float apod_mult = 1; float mult = 0; short usedLines = (maxLine - minLine); //quadratic delay l_i = line / outputL * inputL; for (short sample = 0; sample < outputS; ++sample) { - s_i = sample / outputS * inputS / 2; + s_i = sample / outputS * inputS / (2 - (int)config->GetIsPhotoacousticImage()); part = part_multiplicator*s_i; if (part < 1) part = 1; maxLine = (short)std::min((l_i + part) + 1, inputL); minLine = (short)std::max((l_i - part), 0.0f); usedLines = (maxLine - minLine); apod_mult = (float)apodArraySize / (float)usedLines; delayMultiplicator = pow((1 / (config->GetTimeSpacing()*config->GetSpeedOfSound()) * (config->GetPitchInMeters()*config->GetTransducerElements()) / inputL), 2) / s_i / 2; //calculate the AddSamples beforehand to save some time short* AddSample = new short[maxLine - minLine]; for (short l_s = 0; l_s < maxLine - minLine; ++l_s) { AddSample[l_s] = (short)(delayMultiplicator * pow((minLine + l_s - l_i), 2) + s_i) + (1 - config->GetIsPhotoacousticImage())*s_i; } float s_1 = 0; float s_2 = 0; float sign = 0; for (short l_s1 = minLine; l_s1 < maxLine - 1; ++l_s1) { if (AddSample[l_s1 - minLine] < inputS && AddSample[l_s1 - minLine] >= 0) { s_1 = input[l_s1 + AddSample[l_s1 - minLine] * (short)inputL]; sign += s_1; for (short l_s2 = l_s1 + 1; l_s2 < maxLine; ++l_s2) { if (AddSample[l_s2 - minLine] < inputS && AddSample[l_s2 - minLine] >= 0) { s_2 = input[l_s2 + AddSample[l_s2 - minLine] * (short)inputL]; mult = s_2 * apodisation[(int)((l_s2 - minLine)*apod_mult)] * s_1 * apodisation[(int)((l_s1 - minLine)*apod_mult)]; output[sample*(short)outputL + line] += sqrt(fabs(mult)) * ((mult > 0) - (mult < 0)); } } } else --usedLines; } output[sample*(short)outputL + line] = output[sample*(short)outputL + line] / (float)(pow(usedLines, 2) - (usedLines - 1)) * ((sign > 0) - (sign < 0)); delete[] AddSample; } } void mitk::BeamformingUtils::sDMASSphericalLine( float* input, float* output, float inputDim[2], float outputDim[2], const short& line, const mitk::BeamformingSettings::Pointer config) { const float* apodisation = config->GetApodizationFunction(); const short apodArraySize = config->GetApodizationArraySize(); float& inputS = inputDim[1]; float& inputL = inputDim[0]; float& outputS = outputDim[1]; float& outputL = outputDim[0]; short maxLine = 0; short minLine = 0; float l_i = 0; float s_i = 0; float part = 0.07 * inputL; float tan_phi = std::tan(config->GetAngle() / 360 * 2 * itk::Math::pi); float part_multiplicator = tan_phi * config->GetTimeSpacing() * config->GetSpeedOfSound() / config->GetPitchInMeters() * inputL / (float)config->GetTransducerElements(); float apod_mult = 1; float mult = 0; short usedLines = (maxLine - minLine); //exact delay l_i = (float)line / outputL * inputL; for (short sample = 0; sample < outputS; ++sample) { - s_i = (float)sample / outputS * inputS / 2; + s_i = (float)sample / outputS * inputS / (2 - (int)config->GetIsPhotoacousticImage()); part = part_multiplicator*s_i; if (part < 1) part = 1; maxLine = (short)std::min((l_i + part) + 1, inputL); minLine = (short)std::max((l_i - part), 0.0f); usedLines = (maxLine - minLine); apod_mult = (float)apodArraySize / (float)usedLines; //calculate the AddSamples beforehand to save some time short* AddSample = new short[maxLine - minLine]; for (short l_s = 0; l_s < maxLine - minLine; ++l_s) { AddSample[l_s] = (short)sqrt( pow(s_i, 2) + pow((1 / (config->GetTimeSpacing()*config->GetSpeedOfSound()) * (((float)minLine + (float)l_s - l_i)*config->GetPitchInMeters()*(float)config->GetTransducerElements()) / inputL), 2) ) + (1 - config->GetIsPhotoacousticImage())*s_i; } float s_1 = 0; float s_2 = 0; float sign = 0; for (short l_s1 = minLine; l_s1 < maxLine - 1; ++l_s1) { if (AddSample[l_s1 - minLine] < inputS && AddSample[l_s1 - minLine] >= 0) { s_1 = input[l_s1 + AddSample[l_s1 - minLine] * (short)inputL]; sign += s_1; for (short l_s2 = l_s1 + 1; l_s2 < maxLine; ++l_s2) { if (AddSample[l_s2 - minLine] < inputS && AddSample[l_s2 - minLine] >= 0) { s_2 = input[l_s2 + AddSample[l_s2 - minLine] * (short)inputL]; mult = s_2 * apodisation[(int)((l_s2 - minLine)*apod_mult)] * s_1 * apodisation[(int)((l_s1 - minLine)*apod_mult)]; output[sample*(short)outputL + line] += sqrt(fabs(mult)) * ((mult > 0) - (mult < 0)); } } } else --usedLines; } output[sample*(short)outputL + line] = output[sample*(short)outputL + line] / (float)(pow(usedLines, 2) - (usedLines - 1)) * ((sign > 0) - (sign < 0)); delete[] AddSample; } }