diff --git a/Modules/PhotoacousticsAlgorithms/Resources/DelayCalculation.cl b/Modules/PhotoacousticsAlgorithms/Resources/DelayCalculation.cl index c0ee5224e4..30f9cdabba 100644 --- a/Modules/PhotoacousticsAlgorithms/Resources/DelayCalculation.cl +++ b/Modules/PhotoacousticsAlgorithms/Resources/DelayCalculation.cl @@ -1,70 +1,85 @@ /*=================================================================== 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, - float percentOfImageReconstructed // parameters + float totalSamples_i, + float probeRadius, + float mult1, + float mult2 // 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-isPAImage) * percentOfImageReconstructed; + float s_i = (float)globalPosY / (float)outputS * totalSamples_i; float l_s = (float)globalPosX / (float)outputL * (float)inputL; // the currently calculated line + + float elementHeight = 0; + if (concave) + elementHeight = (1 - cos(abs(inputL / 2.f - l_s) * mult1)) * mult2; - float delayMultiplicator = delayMultiplicatorRaw / s_i; - gDest[globalPosY * (outputL / 2) + globalPosX] = delayMultiplicator * pow((l_s - l_i), 2) + s_i + (1-isPAImage)*s_i; + float delayMultiplicator = delayMultiplicatorRaw / (s_i - elementHeight); + + gDest[globalPosY * (outputL / 2) + globalPosX] = delayMultiplicator * pow((l_s - l_i), 2) + (s_i - elementHeight) + (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, - float percentOfImageReconstructed // parameters + float totalSamples_i, + float probeRadius, + float mult1, + float mult2 // 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-isPAImage) * percentOfImageReconstructed; + float s_i = (float)globalPosY / (float)outputS * totalSamples_i; float l_s = (float)globalPosX / (float)outputL * (float)inputL; // the currently calculated line + float elementHeight = 0; + if (concave) + elementHeight = (1 - cos(abs(inputL / 2.f - l_s) * mult1)) * mult2; + gDest[globalPosY * (outputL / 2) + globalPosX] = sqrt( - pow(s_i, 2) + pow(s_i - elementHeight, 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/mitkPhotoacousticOCLDelayCalculation.cpp b/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLDelayCalculation.cpp index cc236090bc..d9d2d727cc 100644 --- a/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLDelayCalculation.cpp +++ b/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLDelayCalculation.cpp @@ -1,136 +1,140 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #define _USE_MATH_DEFINES #include #include "./OpenCLFilter/mitkPhotoacousticOCLDelayCalculation.h" #include "usServiceReference.h" #include "mitkImageReadAccessor.h" mitk::OCLDelayCalculation::OCLDelayCalculation(mitk::BeamformingSettings::Pointer settings) : m_PixelCalculation(NULL), m_Conf(settings) { this->AddSourceFile("DelayCalculation.cl"); this->m_FilterID = "DelayCalculation"; m_ChunkSize[0] = 128; m_ChunkSize[1] = 128; m_ChunkSize[2] = 8; this->Initialize(); } mitk::OCLDelayCalculation::~OCLDelayCalculation() { if (this->m_PixelCalculation) { clReleaseKernel(m_PixelCalculation); } } void mitk::OCLDelayCalculation::Update() { //Check if context & program available if (!this->Initialize()) { us::ServiceReference ref = GetModuleContext()->GetServiceReference(); OclResourceService* resources = GetModuleContext()->GetService(ref); // clean-up also the resources resources->InvalidateStorage(); mitkThrow() << "Filter is not initialized. Cannot update."; } else { // Execute this->Execute(); } } void mitk::OCLDelayCalculation::Execute() { cl_int clErr = 0; unsigned int gridDim[3] = { m_Conf->GetReconstructionLines() / 2, m_Conf->GetSamplesPerLine(), 1 }; m_BufferSize = gridDim[0] * gridDim[1] * 1; try { this->InitExecNoInput(this->m_PixelCalculation, gridDim, m_BufferSize, sizeof(unsigned short)); } catch (const mitk::Exception& e) { MITK_ERROR << "Caught exception while initializing Delay Calculation filter: " << e.what(); return; } // This calculation is the same for all kernels, so for performance reasons simply perform it here instead of within the kernels if (m_Conf->GetDelayCalculationMethod() == BeamformingSettings::DelayCalc::QuadApprox) m_DelayMultiplicatorRaw = pow(1 / (m_Conf->GetTimeSpacing()*m_Conf->GetSpeedOfSound()) * m_Conf->GetPitchInMeters() * (float)m_Conf->GetTransducerElements() / (float)m_Conf->GetInputDim()[0], 2) / 2; else if (m_Conf->GetDelayCalculationMethod() == BeamformingSettings::DelayCalc::Spherical) m_DelayMultiplicatorRaw = 1 / (m_Conf->GetTimeSpacing()*m_Conf->GetSpeedOfSound()) * (m_Conf->GetPitchInMeters()*(float)m_Conf->GetTransducerElements()); // as openCL does not support bool as a kernel argument, we need to buffer this value in a char... m_IsPAImage = m_Conf->GetIsPhotoacousticImage(); unsigned int reconstructionLines = this->m_Conf->GetReconstructionLines(); unsigned int samplesperLine = this->m_Conf->GetSamplesPerLine(); - float percentOfImageReconstructed = (float)(m_Conf->GetReconstructionDepth()) / - (float)(m_Conf->GetInputDim()[1] * m_Conf->GetSpeedOfSound() * m_Conf->GetTimeSpacing() / (float)(2 - (int)m_Conf->GetIsPhotoacousticImage())); - percentOfImageReconstructed = percentOfImageReconstructed <= 1 ? percentOfImageReconstructed : 1; + float totalSamples_i = (float)(m_Conf->GetReconstructionDepth()) / (float)(m_Conf->GetSpeedOfSound() * m_Conf->GetTimeSpacing()); + totalSamples_i = totalSamples_i <= m_Conf->GetInputDim()[1] ? totalSamples_i : m_Conf->GetInputDim()[1]; + + float probeRadius = m_Conf->GetGeometry() == mitk::BeamformingSettings::ProbeGeometry::Concave ? 0 : m_Conf->GetProbeRadius(); + float mult1 = m_Conf->GetPitchInMeters() / (2 * itk::Math::pi*probeRadius) * 2 * itk::Math::pi; + float mult2 = probeRadius / (m_Conf->GetSpeedOfSound()*m_Conf->GetTimeSpacing()); clErr = clSetKernelArg(this->m_PixelCalculation, 1, sizeof(cl_mem), &(this->m_UsedLines)); clErr |= clSetKernelArg(this->m_PixelCalculation, 2, sizeof(cl_uint), &(this->m_Conf->GetInputDim()[0])); clErr |= clSetKernelArg(this->m_PixelCalculation, 3, sizeof(cl_uint), &(this->m_Conf->GetInputDim()[1])); clErr |= clSetKernelArg(this->m_PixelCalculation, 4, sizeof(cl_uint), &(reconstructionLines)); clErr |= clSetKernelArg(this->m_PixelCalculation, 5, sizeof(cl_uint), &(samplesperLine)); clErr |= clSetKernelArg(this->m_PixelCalculation, 6, sizeof(cl_char), &(this->m_IsPAImage)); clErr |= clSetKernelArg(this->m_PixelCalculation, 7, sizeof(cl_float), &(this->m_DelayMultiplicatorRaw)); - clErr |= clSetKernelArg(this->m_PixelCalculation, 8, sizeof(cl_float), &(percentOfImageReconstructed)); + clErr |= clSetKernelArg(this->m_PixelCalculation, 8, sizeof(cl_float), &(totalSamples_i)); + clErr |= clSetKernelArg(this->m_PixelCalculation, 9, sizeof(cl_float), &(probeRadius)); CHECK_OCL_ERR(clErr); // execute the filter on a 3D NDRange if (!this->ExecuteKernelChunksInBatches(m_PixelCalculation, 2, m_ChunkSize, 16, 50)) mitkThrow() << "openCL Error when executing Kernel"; // signalize the GPU-side data changed m_Output->Modified(GPU_DATA); } us::Module *mitk::OCLDelayCalculation::GetModule() { return us::GetModuleContext()->GetModule(); } bool mitk::OCLDelayCalculation::Initialize() { bool buildErr = true; cl_int clErr = 0; if (OclFilter::Initialize()) { if (m_Conf->GetDelayCalculationMethod() == BeamformingSettings::DelayCalc::QuadApprox) this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "ckDelayCalculationQuad", &clErr); if (m_Conf->GetDelayCalculationMethod() == BeamformingSettings::DelayCalc::Spherical) this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "ckDelayCalculationSphe", &clErr); buildErr |= CHECK_OCL_ERR(clErr); } return (OclFilter::IsInitialized() && buildErr); } diff --git a/Modules/PhotoacousticsAlgorithms/source/utils/mitkBeamformingUtils.cpp b/Modules/PhotoacousticsAlgorithms/source/utils/mitkBeamformingUtils.cpp index d3c087beda..884bbe7f95 100644 --- a/Modules/PhotoacousticsAlgorithms/source/utils/mitkBeamformingUtils.cpp +++ b/Modules/PhotoacousticsAlgorithms/source/utils/mitkBeamformingUtils.cpp @@ -1,609 +1,615 @@ /*=================================================================== 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; float probeRadius = config->GetProbeRadius(); bool concave = config->GetGeometry() == mitk::BeamformingSettings::ProbeGeometry::Concave; float elementHeight = 0; short usedLines = (maxLine - minLine); float totalSamples_i = (float)(config->GetReconstructionDepth()) / (float)(config->GetSpeedOfSound() * config->GetTimeSpacing()); totalSamples_i = totalSamples_i <= inputS ? totalSamples_i : inputS; l_i = (float)line / outputL * inputL; for (short sample = 0; sample < outputS; ++sample) { s_i = (float)sample / outputS * totalSamples_i; part = part_multiplicator*s_i; if (part < 1) part = 1; maxLine = std::round(std::min((l_i + part) + 1, inputL)); minLine = std::round(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) { if (concave) { elementHeight = (probeRadius - std::cos(((abs(inputL / 2.f - l_s)*config->GetPitchInMeters()) / (2 * itk::Math::pi*probeRadius)) * 2 * itk::Math::pi)*probeRadius) / (config->GetSpeedOfSound()*config->GetTimeSpacing()); + delayMultiplicator = pow((1 / (config->GetTimeSpacing()*config->GetSpeedOfSound()) * + (config->GetPitchInMeters()*config->GetTransducerElements()) / inputL), 2) / (s_i - elementHeight) / 2; } AddSample = delayMultiplicator * pow((l_s - l_i), 2) + (s_i - elementHeight) + (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; // just a default value 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 probeRadius = config->GetProbeRadius(); bool concave = config->GetGeometry() == mitk::BeamformingSettings::ProbeGeometry::Concave; float elementHeight = 0; short usedLines = (maxLine - minLine); float totalSamples_i = (float)(config->GetReconstructionDepth()) / (float)(config->GetSpeedOfSound() * config->GetTimeSpacing()); totalSamples_i = totalSamples_i <= inputS ? totalSamples_i : inputS; l_i = (float)line / outputL * inputL; for (short sample = 0; sample < outputS; ++sample) { s_i = (float)sample / outputS * totalSamples_i; part = part_multiplicator*s_i; if (part < 1) part = 1; maxLine = std::round(std::min((l_i + part) + 1, inputL)); minLine = std::round(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) { if (concave) { elementHeight = (probeRadius - std::cos(((abs(inputL / 2.f - l_s)*config->GetPitchInMeters()) / (2 * itk::Math::pi*probeRadius)) * 2 * itk::Math::pi)*probeRadius) / (config->GetSpeedOfSound()*config->GetTimeSpacing()); } AddSample = (int)sqrt( pow(s_i-elementHeight, 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 probeRadius = config->GetProbeRadius(); bool concave = config->GetGeometry() == mitk::BeamformingSettings::ProbeGeometry::Concave; float elementHeight = 0; float mult = 0; short usedLines = (maxLine - minLine); float totalSamples_i = (float)(config->GetReconstructionDepth()) / (float)(config->GetSpeedOfSound() * config->GetTimeSpacing()); totalSamples_i = totalSamples_i <= inputS ? totalSamples_i : inputS; l_i = (float)line / outputL * inputL; for (short sample = 0; sample < outputS; ++sample) { s_i = (float)sample / outputS * totalSamples_i; part = part_multiplicator*s_i; if (part < 1) part = 1; maxLine = std::round(std::min((l_i + part) + 1, inputL)); minLine = std::round(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) { if (concave) { elementHeight = (probeRadius - std::cos(((abs(inputL / 2.f - l_s)*config->GetPitchInMeters()) / (2 * itk::Math::pi*probeRadius)) * 2 * itk::Math::pi)*probeRadius) / (config->GetSpeedOfSound()*config->GetTimeSpacing()); + delayMultiplicator = pow((1 / (config->GetTimeSpacing()*config->GetSpeedOfSound()) * + (config->GetPitchInMeters()*config->GetTransducerElements()) / inputL), 2) / (s_i - elementHeight) / 2; } AddSample[l_s] = (short)(delayMultiplicator * pow((minLine + l_s - l_i), 2) + (s_i - elementHeight)) + (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 probeRadius = config->GetProbeRadius(); bool concave = config->GetGeometry() == mitk::BeamformingSettings::ProbeGeometry::Concave; float elementHeight = 0; float mult = 0; short usedLines = (maxLine - minLine); float totalSamples_i = (float)(config->GetReconstructionDepth()) / (float)(config->GetSpeedOfSound() * config->GetTimeSpacing()); totalSamples_i = totalSamples_i <= inputS ? totalSamples_i : inputS; l_i = (float)line / outputL * inputL; for (short sample = 0; sample < outputS; ++sample) { s_i = (float)sample / outputS * totalSamples_i; part = part_multiplicator*s_i; if (part < 1) part = 1; maxLine = std::round(std::min((l_i + part) + 1, inputL)); minLine = std::round(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) { if (concave) { elementHeight = (probeRadius - std::cos(((abs(inputL / 2.f - l_s)*config->GetPitchInMeters()) / (2 * itk::Math::pi*probeRadius)) * 2 * itk::Math::pi)*probeRadius) / (config->GetSpeedOfSound()*config->GetTimeSpacing()); } AddSample[l_s] = (short)sqrt( pow(s_i - elementHeight, 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 probeRadius = config->GetProbeRadius(); bool concave = config->GetGeometry() == mitk::BeamformingSettings::ProbeGeometry::Concave; float elementHeight = 0; float mult = 0; short usedLines = (maxLine - minLine); float totalSamples_i = (float)(config->GetReconstructionDepth()) / (float)(config->GetSpeedOfSound() * config->GetTimeSpacing()); totalSamples_i = totalSamples_i <= inputS ? totalSamples_i : inputS; l_i = (float)line / outputL * inputL; for (short sample = 0; sample < outputS; ++sample) { s_i = (float)sample / outputS * totalSamples_i; part = part_multiplicator*s_i; if (part < 1) part = 1; maxLine = std::round(std::min((l_i + part) + 1, inputL)); minLine = std::round(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) { if (concave) { elementHeight = (probeRadius - std::cos(((abs(inputL / 2.f - l_s)*config->GetPitchInMeters()) / (2 * itk::Math::pi*probeRadius)) * 2 * itk::Math::pi)*probeRadius) / (config->GetSpeedOfSound()*config->GetTimeSpacing()); + delayMultiplicator = pow((1 / (config->GetTimeSpacing()*config->GetSpeedOfSound()) * + (config->GetPitchInMeters()*config->GetTransducerElements()) / inputL), 2) / (s_i - elementHeight) / 2; } AddSample[l_s] = (short)(delayMultiplicator * pow((minLine + l_s - l_i), 2) + (s_i - elementHeight)) + (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 probeRadius = config->GetProbeRadius(); bool concave = config->GetGeometry() == mitk::BeamformingSettings::ProbeGeometry::Concave; float elementHeight = 0; float mult = 0; short usedLines = (maxLine - minLine); float totalSamples_i = (float)(config->GetReconstructionDepth()) / (float)(config->GetSpeedOfSound() * config->GetTimeSpacing()); totalSamples_i = totalSamples_i <= inputS ? totalSamples_i : inputS; l_i = (float)line / outputL * inputL; for (short sample = 0; sample < outputS; ++sample) { s_i = (float)sample / outputS * totalSamples_i; part = part_multiplicator*s_i; if (part < 1) part = 1; maxLine = std::round(std::min((l_i + part) + 1, inputL)); minLine = std::round(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) { if (concave) { elementHeight = (probeRadius - std::cos(((abs(inputL / 2.f - l_s)*config->GetPitchInMeters()) / (2 * itk::Math::pi*probeRadius)) * 2 * itk::Math::pi)*probeRadius) / (config->GetSpeedOfSound()*config->GetTimeSpacing()); } AddSample[l_s] = (short)sqrt( pow(s_i - elementHeight, 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; } }