diff --git a/Modules/PhotoacousticsAlgorithms/Resources/UsedLinesCalculation.cl b/Modules/PhotoacousticsAlgorithms/Resources/UsedLinesCalculation.cl index 132c7d88fd..866b5f5d21 100644 --- a/Modules/PhotoacousticsAlgorithms/Resources/UsedLinesCalculation.cl +++ b/Modules/PhotoacousticsAlgorithms/Resources/UsedLinesCalculation.cl @@ -1,170 +1,132 @@ /*=================================================================== 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 ckUsedLines( __global unsigned short* dDest, // output buffer float partMult, unsigned int inputL, unsigned int inputS, unsigned int outputL, unsigned int outputS, float totalSamples_i // parameters ) { // get thread identifier unsigned int globalPosX = get_global_id(0); unsigned int globalPosY = get_global_id(1); // terminate non-valid threads if ( globalPosX < outputL && globalPosY < outputS) { float l_i = (float)globalPosX / outputL * inputL; float s_i = (float)globalPosY / (float)outputS * totalSamples_i; float part = partMult * s_i; if (part < 1) part = 1; unsigned short maxLine = min((l_i + part) + 1, (float)inputL); unsigned short minLine = max((l_i - part), 0.0f); dDest[globalPosY * 3 * outputL + 3 * globalPosX] = (maxLine - minLine); //usedLines dDest[globalPosY * 3 * outputL + 3 * globalPosX + 1] = minLine; //minLine dDest[globalPosY * 3 * outputL + 3 * globalPosX + 2] = maxLine; //maxLine } } __kernel void ckUsedLines_g( __global unsigned short* dDest, // output buffer __global float* elementHeights, __global float* elementPositions, - float cos_deg, + float sin_deg, float probeRadius, unsigned int inputL, unsigned int inputS, unsigned int outputL, unsigned int outputS, float horizontalExtent, float verticalExtent ) { // get thread identifier unsigned int globalPosX = get_global_id(0); unsigned int globalPosY = get_global_id(1); + float x_sensor_pos = 0; + float y_sensor_pos = 0; + float center_to_sensor_a = 0; + float center_to_sensor_b = 0; + float center_to_sensor_c = 0; + float distance_to_sensor_direction = 0; + float distance_sensor_target = 0; + // terminate non-valid threads if ( globalPosX < outputL && globalPosY < outputS) { - float cos = 0; - float a = 0; - float d = 0; - - float l_p = (float)globalPosX / outputL * horizontalExtent; - float s_p = (float)globalPosY / (float)outputS * verticalExtent; - + float x_cm = (float)globalPosX / outputL * horizontalExtent; // x*Xspacing + float y_cm = (float)globalPosY / (float)outputS * verticalExtent; // y*Yspacing + int maxLine = inputL; int minLine = 0; - - for(int l_s = 0; l_s < inputL; l_s+=32) - { - a = sqrt((probeRadius - s_p)*(probeRadius - s_p) + (l_p - horizontalExtent/2)*(l_p - horizontalExtent/2)); - d = sqrt((s_p - elementHeights[l_s])*(s_p - elementHeights[l_s]) + (l_p - elementPositions[l_s])*(l_p - elementPositions[l_s])); - cos = (d*d + probeRadius*probeRadius - a*a)/(2*probeRadius*d); - - if(cos > cos_deg) - { - minLine = l_s-32; - if(minLine < 0) - minLine = 0; - break; - } - } - for(int l_s = minLine; l_s < inputL; l_s+=8) - { - a = sqrt((probeRadius - s_p)*(probeRadius - s_p) + (l_p - horizontalExtent/2)*(l_p - horizontalExtent/2)); - d = sqrt((s_p - elementHeights[l_s])*(s_p - elementHeights[l_s]) + (l_p - elementPositions[l_s])*(l_p - elementPositions[l_s])); - cos = (d*d + probeRadius*probeRadius - a*a)/(2*probeRadius*d); - - if(cos > cos_deg) - { - minLine = l_s-8; - if(minLine < 0) - minLine = 0; - break; - } - } - for(int l_s = minLine; l_s < inputL; l_s+=1) + + float x_center_pos = horizontalExtent / 2.0; + float y_center_pos = probeRadius; + + for (int l_s = minLine; l_s <= inputL; l_s += 1) { - a = sqrt((probeRadius - s_p)*(probeRadius - s_p) + (l_p - horizontalExtent/2)*(l_p - horizontalExtent/2)); - d = sqrt((s_p - elementHeights[l_s])*(s_p - elementHeights[l_s]) + (l_p - elementPositions[l_s])*(l_p - elementPositions[l_s])); - cos = (d*d + probeRadius*probeRadius - a*a)/(2*probeRadius*d); - - if(cos > cos_deg) + x_sensor_pos = elementPositions[l_s]; + y_sensor_pos = elementHeights[l_s]; + + distance_sensor_target = sqrt((x_cm - x_sensor_pos)*(x_cm - x_sensor_pos) + (y_cm - y_sensor_pos)*(y_cm - y_sensor_pos)); + + // solving line equation + center_to_sensor_a = y_sensor_pos - y_center_pos; + center_to_sensor_b = x_center_pos - x_sensor_pos; + center_to_sensor_c = -(center_to_sensor_a * x_center_pos + center_to_sensor_b * y_center_pos); + distance_to_sensor_direction = fabs((center_to_sensor_a * x_cm + center_to_sensor_b * y_cm + center_to_sensor_c)) / (sqrt(center_to_sensor_a*center_to_sensor_a + center_to_sensor_b*center_to_sensor_b)); + + if (distance_to_sensor_direction < (sin_deg*distance_sensor_target)) { minLine = l_s; break; } } - - for(int l_s = inputL; l_s >= 0 ; l_s-=32) - { - a = sqrt((probeRadius - s_p)*(probeRadius - s_p) + (l_p - horizontalExtent/2)*(l_p - horizontalExtent/2)); - d = sqrt((s_p - elementHeights[l_s])*(s_p - elementHeights[l_s]) + (l_p - elementPositions[l_s])*(l_p - elementPositions[l_s])); - cos = (d*d + probeRadius*probeRadius - a*a)/(2*probeRadius*d); - cos = 0; - - if(cos > cos_deg) - { - maxLine = l_s+32; - if(maxLine > inputL) - minLine = inputL; - break; - } - } - for(int l_s = maxLine; l_s >= 0 ; l_s-=8) - { - a = sqrt((probeRadius - s_p)*(probeRadius - s_p) + (l_p - horizontalExtent/2)*(l_p - horizontalExtent/2)); - d = sqrt((s_p - elementHeights[l_s])*(s_p - elementHeights[l_s]) + (l_p - elementPositions[l_s])*(l_p - elementPositions[l_s])); - cos = (d*d + probeRadius*probeRadius - a*a)/(2*probeRadius*d); - cos = 0; - - if(cos > cos_deg) - { - maxLine = l_s+8; - if(maxLine > inputL) - minLine = inputL; - break; - } - } - for(int l_s = maxLine; l_s >= 0 ; l_s-=1) - { - a = sqrt((probeRadius - s_p)*(probeRadius - s_p) + (l_p - horizontalExtent/2)*(l_p - horizontalExtent/2)); - d = sqrt((s_p - elementHeights[l_s])*(s_p - elementHeights[l_s]) + (l_p - elementPositions[l_s])*(l_p - elementPositions[l_s])); - cos = (d*d + probeRadius*probeRadius - a*a)/(2*probeRadius*d); - cos = 0; - - if(cos > cos_deg) + + for (int l_s = maxLine - 1; l_s > minLine; l_s -= 1) // start with maxline-1 otherwise elementPositions[] will go out of range + { + x_sensor_pos = elementPositions[l_s]; + y_sensor_pos = elementHeights[l_s]; + + distance_sensor_target = sqrt((x_cm - x_sensor_pos)*(x_cm - x_sensor_pos) + (y_cm - y_sensor_pos)*(y_cm - y_sensor_pos)); + + // solving line equation + center_to_sensor_a = y_sensor_pos - y_center_pos; + center_to_sensor_b = x_center_pos - x_sensor_pos; + center_to_sensor_c = -(center_to_sensor_a * x_center_pos + center_to_sensor_b * y_center_pos); + distance_to_sensor_direction = fabs((center_to_sensor_a * x_cm + center_to_sensor_b * y_cm + center_to_sensor_c)) / (sqrt(center_to_sensor_a*center_to_sensor_a + center_to_sensor_b*center_to_sensor_b)); + + if (distance_to_sensor_direction < sin_deg*distance_sensor_target) { maxLine = l_s; break; } } - + dDest[globalPosY * 3 * outputL + 3 * globalPosX] = (maxLine - minLine); //usedLines dDest[globalPosY * 3 * outputL + 3 * globalPosX + 1] = minLine; //minLine dDest[globalPosY * 3 * outputL + 3 * globalPosX + 2] = maxLine; //maxLine } -} \ No newline at end of file +} diff --git a/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLUsedLinesCalculation.cpp b/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLUsedLinesCalculation.cpp index 8da3c4c197..d671f523cf 100644 --- a/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLUsedLinesCalculation.cpp +++ b/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLUsedLinesCalculation.cpp @@ -1,167 +1,165 @@ /*=================================================================== 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/mitkPhotoacousticOCLUsedLinesCalculation.h" #include "usServiceReference.h" #include "mitkImageReadAccessor.h" mitk::OCLUsedLinesCalculation::OCLUsedLinesCalculation(mitk::BeamformingSettings::Pointer settings) : m_PixelCalculation(NULL), m_Conf(settings) { this->AddSourceFile("UsedLinesCalculation.cl"); this->m_FilterID = "UsedLinesCalculation"; m_ChunkSize[0] = 128; m_ChunkSize[1] = 128; m_ChunkSize[2] = 8; this->Initialize(); } mitk::OCLUsedLinesCalculation::~OCLUsedLinesCalculation() { if (this->m_PixelCalculation) { clReleaseKernel(m_PixelCalculation); } } void mitk::OCLUsedLinesCalculation::SetElementHeightsBuffer(cl_mem elementHeightsBuffer) { m_ElementHeightsBuffer = elementHeightsBuffer; } void mitk::OCLUsedLinesCalculation::SetElementPositionsBuffer(cl_mem elementPositionsBuffer) { m_ElementPositionsBuffer = elementPositionsBuffer; } void mitk::OCLUsedLinesCalculation::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::OCLUsedLinesCalculation::Execute() { cl_int clErr = 0; unsigned int gridDim[3] = { m_Conf->GetReconstructionLines(), m_Conf->GetSamplesPerLine(), 1 }; size_t outputSize = gridDim[0] * gridDim[1] * 3; try { this->InitExecNoInput(this->m_PixelCalculation, gridDim, outputSize, sizeof(unsigned short)); } catch (const mitk::Exception& e) { MITK_ERROR << "Caught exception while initializing UsedLines filter: " << e.what(); return; } if (m_Conf->GetGeometry() == mitk::BeamformingSettings::ProbeGeometry::Linear) { // This calculation is the same for all kernels, so for performance reasons simply perform it here instead of within the kernels m_part = (tan(m_Conf->GetAngle() / 360 * 2 * itk::Math::pi) * ((m_Conf->GetSpeedOfSound() * m_Conf->GetTimeSpacing())) / (m_Conf->GetPitchInMeters() * m_Conf->GetTransducerElements())) * m_Conf->GetInputDim()[0]; - unsigned int reconLines = this->m_Conf->GetReconstructionLines(); unsigned int samplesPerLine = this->m_Conf->GetSamplesPerLine(); 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]; - clErr = clSetKernelArg(this->m_PixelCalculation, 1, sizeof(cl_float), &(this->m_part)); 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), &(reconLines)); clErr |= clSetKernelArg(this->m_PixelCalculation, 5, sizeof(cl_uint), &(samplesPerLine)); clErr |= clSetKernelArg(this->m_PixelCalculation, 6, sizeof(cl_float), &(totalSamples_i)); } else { unsigned int reconLines = this->m_Conf->GetReconstructionLines(); unsigned int samplesPerLine = this->m_Conf->GetSamplesPerLine(); float probeRadius = m_Conf->GetProbeRadius(); - float cos_deg = std::cos(m_Conf->GetAngle()/2.f / 360 * 2 * itk::Math::pi); + float sin_deg = std::sin(m_Conf->GetAngle() / 360 * 2 * itk::Math::pi); float horizontalExtent = m_Conf->GetHorizontalExtent(); float verticalExtent = m_Conf->GetReconstructionDepth(); clErr = clSetKernelArg(this->m_PixelCalculation, 1, sizeof(cl_mem), &(this->m_ElementHeightsBuffer)); clErr |= clSetKernelArg(this->m_PixelCalculation, 2, sizeof(cl_mem), &(this->m_ElementPositionsBuffer)); - clErr |= clSetKernelArg(this->m_PixelCalculation, 3, sizeof(cl_float), &(cos_deg)); + clErr |= clSetKernelArg(this->m_PixelCalculation, 3, sizeof(cl_float), &(sin_deg)); clErr |= clSetKernelArg(this->m_PixelCalculation, 4, sizeof(cl_float), &(probeRadius)); clErr |= clSetKernelArg(this->m_PixelCalculation, 5, sizeof(cl_uint), &(this->m_Conf->GetInputDim()[0])); clErr |= clSetKernelArg(this->m_PixelCalculation, 6, sizeof(cl_uint), &(this->m_Conf->GetInputDim()[1])); clErr |= clSetKernelArg(this->m_PixelCalculation, 7, sizeof(cl_uint), &(reconLines)); clErr |= clSetKernelArg(this->m_PixelCalculation, 8, sizeof(cl_uint), &(samplesPerLine)); clErr |= clSetKernelArg(this->m_PixelCalculation, 9, sizeof(cl_float), &(horizontalExtent)); clErr |= clSetKernelArg(this->m_PixelCalculation, 10, sizeof(cl_float), &(verticalExtent)); } CHECK_OCL_ERR(clErr); // execute the filter on a 2D 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::OCLUsedLinesCalculation::GetModule() { return us::GetModuleContext()->GetModule(); } bool mitk::OCLUsedLinesCalculation::Initialize() { bool buildErr = true; cl_int clErr = 0; if (OclFilter::Initialize()) { if (m_Conf->GetGeometry() == mitk::BeamformingSettings::ProbeGeometry::Linear) { this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "ckUsedLines", &clErr); buildErr |= CHECK_OCL_ERR(clErr); } else { this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "ckUsedLines_g", &clErr); buildErr |= CHECK_OCL_ERR(clErr); } } return (OclFilter::IsInitialized() && buildErr); } #endif diff --git a/Modules/PhotoacousticsAlgorithms/source/utils/mitkBeamformingUtils.cpp b/Modules/PhotoacousticsAlgorithms/source/utils/mitkBeamformingUtils.cpp index 901b7f344d..8cfd04f92b 100644 --- a/Modules/PhotoacousticsAlgorithms/source/utils/mitkBeamformingUtils.cpp +++ b/Modules/PhotoacousticsAlgorithms/source/utils/mitkBeamformingUtils.cpp @@ -1,414 +1,412 @@ /*=================================================================== 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; } unsigned short* mitk::BeamformingUtils::MinMaxLines(const mitk::BeamformingSettings::Pointer config) { unsigned int outputL = (unsigned int)config->GetReconstructionLines(); unsigned int outputS = (unsigned int)config->GetSamplesPerLine(); unsigned short* dDest = new unsigned short[outputL * outputS * 2]; unsigned int inputL = (unsigned int)config->GetInputDim()[0]; float horizontalExtent = config->GetHorizontalExtent(); float verticalExtent = config->GetReconstructionDepth(); float partMult = (tan(config->GetAngle() / 360 * 2 * itk::Math::pi) * ((config->GetSpeedOfSound() * config->GetTimeSpacing())) / (config->GetPitchInMeters() * config->GetTransducerElements())) * inputL; float totalSamples_i = (float)(config->GetReconstructionDepth()) / (float)(config->GetSpeedOfSound() * config->GetTimeSpacing()); totalSamples_i = totalSamples_i <= config->GetInputDim()[1] ? totalSamples_i : config->GetInputDim()[1]; if ((int)config->GetGeometry() == 0) // if this is raw data from a linear probe geometry { for (unsigned int x = 0; x < outputL; ++x) { for (unsigned int y = 0; y < outputS; ++y) { float l_i = (float)x / outputL * inputL; float s_i = (float)y / (float)outputS * totalSamples_i; float part = partMult * s_i; if (part < 1) part = 1; unsigned short maxLine = std::min((l_i + part) + 1, (float)inputL); unsigned short minLine = std::max((l_i - part), 0.0f); dDest[y * 2 * outputL + 2 * x] = (unsigned short)minLine; //minLine dDest[y * 2 * outputL + 2 * x + 1] = (unsigned short)maxLine; //maxLine } } } else // if this is *not* raw data from a linear probe geometry (currently meaning its a concave geometry) { float probeRadius = config->GetProbeRadius(); float* elementHeights = config->GetElementHeights(); float* elementPositions = config->GetElementPositions(); float sin_deg = std::sin(config->GetAngle() / 360 * 2 * itk::Math::pi); float x_center_pos = horizontalExtent / 2.0; float y_center_pos = probeRadius; for (unsigned int x = 0; x < outputL; ++x) { for (unsigned int y = 0; y < outputS; ++y) { float x_cm = (float)x / outputL * horizontalExtent; // x*Xspacing float y_cm = (float)y / (float)outputS * verticalExtent; // y*Yspacing unsigned int maxLine = inputL; unsigned int minLine = 0; for (unsigned int l_s = minLine; l_s <= inputL; l_s += 1) { float x_sensor_pos = elementPositions[l_s]; float y_sensor_pos = elementHeights[l_s]; float distance_sensor_target = sqrt((x_cm - x_sensor_pos)*(x_cm - x_sensor_pos) + (y_cm - y_sensor_pos)*(y_cm - y_sensor_pos)); // solving line equation float center_to_sensor_a = y_sensor_pos - y_center_pos; float center_to_sensor_b = x_center_pos - x_sensor_pos; float center_to_sensor_c = -(center_to_sensor_a * x_center_pos + center_to_sensor_b * y_center_pos); float distance_to_sensor_direction = std::fabs((center_to_sensor_a * x_cm + center_to_sensor_b * y_cm + center_to_sensor_c)) / (sqrt(center_to_sensor_a*center_to_sensor_a + center_to_sensor_b*center_to_sensor_b)); if (distance_to_sensor_direction < sin_deg*distance_sensor_target) { minLine = l_s; break; } } for (unsigned int l_s = maxLine - 1; l_s > minLine; l_s -= 1) // start with maxline-1 otherwise elementPositions[] will go out of range { float x_sensor_pos = elementPositions[l_s]; float y_sensor_pos = elementHeights[l_s]; float distance_sensor_target = sqrt((x_cm - x_sensor_pos)*(x_cm - x_sensor_pos) + (y_cm - y_sensor_pos)*(y_cm - y_sensor_pos)); // solving line equation - float center_pos_x = horizontalExtent / 2.0; - float center_pos_y = probeRadius; - float center_to_sensor_a = y_sensor_pos - center_pos_y; - float center_to_sensor_b = center_pos_x - x_sensor_pos; - float center_to_sensor_c = -(center_to_sensor_a * center_pos_x + center_to_sensor_b * center_pos_y); + float center_to_sensor_a = y_sensor_pos - y_center_pos; + float center_to_sensor_b = x_center_pos - x_sensor_pos; + float center_to_sensor_c = -(center_to_sensor_a * x_center_pos + center_to_sensor_b * y_center_pos); float distance_to_sensor_direction = std::fabs((center_to_sensor_a * x_cm + center_to_sensor_b * y_cm + center_to_sensor_c)) / (sqrt(center_to_sensor_a*center_to_sensor_a + center_to_sensor_b*center_to_sensor_b)); if (distance_to_sensor_direction < sin_deg*distance_sensor_target) { maxLine = l_s; break; } } dDest[y * 2 * outputL + 2 * x] = (unsigned short)minLine; //minLine dDest[y * 2 * outputL + 2 * x + 1] = (unsigned short)maxLine; //maxLine } } } return dDest; } 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(); const float* elementHeights = config->GetElementHeights(); const float* elementPositions = config->GetElementPositions(); 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_p = 0; float s_i = 0; float apod_mult = 1; short usedLines = (maxLine - minLine); float totalSamples_i = (float)(config->GetReconstructionDepth()) / (float)(config->GetSpeedOfSound() * config->GetTimeSpacing()); totalSamples_i = totalSamples_i <= inputS ? totalSamples_i : inputS; l_p = (float)line / outputL * config->GetHorizontalExtent(); for (short sample = 0; sample < outputS; ++sample) { s_i = (float)sample / outputS * totalSamples_i; minLine = config->GetMinMaxLines()[2*sample*(short)outputL + 2*line]; maxLine = config->GetMinMaxLines()[2*sample*(short)outputL + 2*line + 1]; 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-elementHeights[l_s]/(config->GetSpeedOfSound()*config->GetTimeSpacing()), 2) + pow((1 / (config->GetTimeSpacing()*config->GetSpeedOfSound())) * (l_p - elementPositions[l_s]), 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::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(); const float* elementHeights = config->GetElementHeights(); const float* elementPositions = config->GetElementPositions(); float& inputS = inputDim[1]; float& inputL = inputDim[0]; float& outputS = outputDim[1]; float& outputL = outputDim[0]; short maxLine = 0; short minLine = 0; float l_p = 0; float s_i = 0; float apod_mult = 1; 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_p = (float)line / outputL * config->GetHorizontalExtent(); for (short sample = 0; sample < outputS; ++sample) { s_i = (float)sample / outputS * totalSamples_i; minLine = config->GetMinMaxLines()[2 * sample*(short)outputL + 2 * line]; maxLine = config->GetMinMaxLines()[2 * sample*(short)outputL + 2 * line + 1]; 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] = (int)sqrt( pow(s_i - elementHeights[l_s + minLine] / (config->GetSpeedOfSound()*config->GetTimeSpacing()), 2) + pow((1 / (config->GetTimeSpacing()*config->GetSpeedOfSound())) * (l_p - elementPositions[l_s + minLine]), 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::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(); const float* elementHeights = config->GetElementHeights(); const float* elementPositions = config->GetElementPositions(); float& inputS = inputDim[1]; float& inputL = inputDim[0]; float& outputS = outputDim[1]; float& outputL = outputDim[0]; short maxLine = 0; short minLine = 0; float l_p = 0; float s_i = 0; float apod_mult = 1; 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_p = (float)line / outputL * config->GetHorizontalExtent(); for (short sample = 0; sample < outputS; ++sample) { s_i = (float)sample / outputS * totalSamples_i; minLine = config->GetMinMaxLines()[2 * sample*(short)outputL + 2 * line]; maxLine = config->GetMinMaxLines()[2 * sample*(short)outputL + 2 * line + 1]; 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] = (int)sqrt( pow(s_i - elementHeights[l_s + minLine] / (config->GetSpeedOfSound()*config->GetTimeSpacing()), 2) + pow((1 / (config->GetTimeSpacing()*config->GetSpeedOfSound())) * (l_p - elementPositions[l_s + minLine]), 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; } }