diff --git a/Modules/PhotoacousticsAlgorithms/MitkPABeamformingTool/PABeamformingTool.cpp b/Modules/PhotoacousticsAlgorithms/MitkPABeamformingTool/PABeamformingTool.cpp index 04663d930e..83eb8d5540 100644 --- a/Modules/PhotoacousticsAlgorithms/MitkPABeamformingTool/PABeamformingTool.cpp +++ b/Modules/PhotoacousticsAlgorithms/MitkPABeamformingTool/PABeamformingTool.cpp @@ -1,378 +1,388 @@ /*=================================================================== 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 #include #include #include #include #include #include #include #include #include struct InputParameters { mitk::Image::Pointer inputImage; std::string outputFilename; bool verbose; std::string settingsFile; std::string imageType; }; struct BandpassSettings { float highPass; float lowPass; float alphaLow; float alphaHigh; float speedOfSound; }; struct CropSettings { int above; int below; int right; int left; int zStart; int zEnd; }; struct ResampleSettings { double spacing; int dimX; }; struct BModeSettings { mitk::PhotoacousticFilterService::BModeMethod method; bool UseLogFilter; }; struct ProcessSettings { bool DoBeamforming; bool DoBandpass; bool DoCropping; bool DoResampling; bool DoBmode; }; InputParameters parseInput(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setCategory("MITK-Photoacoustics"); parser.setTitle("Mitk Photoacoustics Beamforming Tool"); parser.setDescription("Reads a nrrd file as an input and applies a beamforming method as set with the parameters defined in an additionally provided xml file."); parser.setContributor("Computer Assisted Medical Interventions, DKFZ"); parser.setArgumentPrefix("--", "-"); parser.beginGroup("Required parameters"); parser.addArgument( "inputImage", "i", mitkCommandLineParser::Image, "Input image (mitk::Image)", "input image (.nrrd file)", us::Any(), false, false, false, mitkCommandLineParser::Input); parser.addArgument( "output", "o", mitkCommandLineParser::File, "Output filename", "output image (.nrrd file)", us::Any(), false, false, false, mitkCommandLineParser::Output); parser.addArgument( "settings", "s", mitkCommandLineParser::String, "settings file", "file containing beamforming and other specifications(.xml file)", us::Any(), false); parser.addArgument( "type", "t", mitkCommandLineParser::String, "image type", "Specifies whether to use the PA or US subsection of the xml file must be 'PA' or 'US'", us::Any(), false); parser.endGroup(); parser.beginGroup("Optional parameters"); parser.addArgument( "verbose", "v", mitkCommandLineParser::Bool, "Verbose Output", "Whether to produce verbose, or rather debug output. (default: false)"); parser.endGroup(); InputParameters input; std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size() == 0) exit(-1); input.verbose = (bool)parsedArgs.count("verbose"); MITK_INFO(input.verbose) << "### VERBOSE OUTPUT ENABLED ###"; if (parsedArgs.count("inputImage")) { MITK_INFO(input.verbose) << "Reading input image..."; input.inputImage = mitk::IOUtil::Load(us::any_cast(parsedArgs["inputImage"])); MITK_INFO(input.verbose) << "Reading input image...[Done]"; } else mitkThrow() << "No input image given."; if (parsedArgs.count("output")) input.outputFilename = us::any_cast(parsedArgs["output"]); else mitkThrow() << "No output image path given.."; if (parsedArgs.count("settings")) input.settingsFile = us::any_cast(parsedArgs["settings"]); else mitkThrow() << "No settings image path given.."; if (parsedArgs.count("type")) input.imageType = us::any_cast(parsedArgs["type"]); else mitkThrow() << "No settings image type given.."; return input; } TiXmlElement* getRootChild(TiXmlElement* root, std::string childName) { for (TiXmlElement* elem = root->FirstChildElement(); elem != NULL; elem = elem->NextSiblingElement()) { if (elem->Value() == childName) return elem; } return nullptr; } void ParseXML(std::string xmlFile, InputParameters input, mitk::BeamformingSettings::Pointer *bfSet, BandpassSettings& bandpassSet, CropSettings& cropSet, ResampleSettings& resSet, BModeSettings& bmodeSet, ProcessSettings& processSet) { MITK_INFO << "Loading configuration File \"" << xmlFile << "\""; TiXmlDocument doc(xmlFile); if (!doc.LoadFile()) mitkThrow() << "Failed to load settings file \"" << xmlFile << "\" Error: " << doc.ErrorDesc(); TiXmlElement* root = doc.FirstChildElement(); if (root == NULL) { mitkThrow() << "Failed to load file: No root element."; doc.Clear(); } TiXmlElement* elemtop = getRootChild(root, input.imageType); if(elemtop == nullptr) mitkThrow() << "The xml file is wrongly formatted, no element \"" << input.imageType << "\" found"; for (TiXmlElement* elem = elemtop->FirstChildElement(); elem != NULL; elem = elem->NextSiblingElement()) { std::string elemName = elem->Value(); if (elemName == "Beamforming") { processSet.DoBeamforming = std::stoi(elem->Attribute("do")); if (!processSet.DoBeamforming) continue; float PitchInMeters = std::stof(elem->Attribute("pitchInMeters")); float SpeedOfSound = std::stof(elem->Attribute("speedOfSound")); float Angle = std::stof(elem->Attribute("angle")); bool IsPhotoacousticImage = std::stoi(elem->Attribute("isPhotoacousticImage")); unsigned int SamplesPerLine = std::stoi(elem->Attribute("samplesPerLine")); unsigned int ReconstructionLines = std::stoi(elem->Attribute("reconstructionLines")); float ReconstructionDepth = std::stof(elem->Attribute("reconstructionDepth")); bool UseGPU = std::stoi(elem->Attribute("useGPU")); unsigned int GPUBatchSize = std::stoi(elem->Attribute("GPUBatchSize")); std::string apodizationStr = elem->Attribute("apodization"); mitk::BeamformingSettings::Apodization Apodization = mitk::BeamformingSettings::Apodization::Box; if (apodizationStr == "Box") Apodization = mitk::BeamformingSettings::Apodization::Box; else if (apodizationStr == "Hann") Apodization = mitk::BeamformingSettings::Apodization::Hann; else if (apodizationStr == "Hamm") Apodization = mitk::BeamformingSettings::Apodization::Hamm; else mitkThrow() << "Apodization incorrectly defined in settings"; + std::string geomStr = elem->Attribute("geometry"); + mitk::BeamformingSettings::ProbeGeometry ProbeGeometry = mitk::BeamformingSettings::ProbeGeometry::Linear; + if (geomStr == "linear") + ProbeGeometry = mitk::BeamformingSettings::ProbeGeometry::Linear; + else if(geomStr == "concave") + ProbeGeometry = mitk::BeamformingSettings::ProbeGeometry::Concave; + else + mitkThrow() << "geometry incorrectly defined in settings"; + + float radius = std::stof(elem->Attribute("radiusInMm")); + unsigned int ApodizationArraySize = std::stoi(elem->Attribute("apodizationArraySize")); std::string algorithmStr = elem->Attribute("algorithm"); mitk::BeamformingSettings::BeamformingAlgorithm Algorithm = mitk::BeamformingSettings::BeamformingAlgorithm::DAS; if (algorithmStr == "DAS") Algorithm = mitk::BeamformingSettings::BeamformingAlgorithm::DAS; else if (algorithmStr == "DMAS") Algorithm = mitk::BeamformingSettings::BeamformingAlgorithm::DMAS; else if (algorithmStr == "sDMAS") Algorithm = mitk::BeamformingSettings::BeamformingAlgorithm::sDMAS; else - { mitkThrow() << "Beamforming algorithm incorrectly defined in settings"; - } *bfSet = mitk::BeamformingSettings::New( (float)(input.inputImage->GetGeometry()->GetSpacing()[0] / 1000), SpeedOfSound, (float)(input.inputImage->GetGeometry()->GetSpacing()[1] / 1000000), Angle, IsPhotoacousticImage, SamplesPerLine, ReconstructionLines, input.inputImage->GetDimensions(), ReconstructionDepth, UseGPU, GPUBatchSize, - mitk::BeamformingSettings::DelayCalc::Spherical, Apodization, ApodizationArraySize, - Algorithm + Algorithm, + ProbeGeometry, + radius ); } if (elemName == "Bandpass") { processSet.DoBandpass = std::stoi(elem->Attribute("do")); if (!processSet.DoBandpass) continue; bandpassSet.highPass = std::stof(elem->Attribute("highPass")); bandpassSet.lowPass = std::stof(elem->Attribute("lowPass")); bandpassSet.alphaLow = std::stof(elem->Attribute("alphaLow")); bandpassSet.alphaHigh = std::stof(elem->Attribute("alphaHigh")); bandpassSet.speedOfSound = std::stof(elem->Attribute("speedOfSound")); } if (elemName == "Cropping") { processSet.DoCropping = std::stoi(elem->Attribute("do")); if (!processSet.DoCropping) continue; cropSet.above = std::stoi(elem->Attribute("cutAbove")); cropSet.below = std::stoi(elem->Attribute("cutBelow")); cropSet.right = std::stoi(elem->Attribute("cutRight")); cropSet.left = std::stoi(elem->Attribute("cutLeft")); cropSet.zStart = std::stoi(elem->Attribute("firstSlice")); cropSet.zEnd = std::stoi(elem->Attribute("cutSlices")); } if (elemName == "Resampling") { processSet.DoResampling = std::stoi(elem->Attribute("do")); if (!processSet.DoResampling) continue; resSet.spacing = std::stod(elem->Attribute("spacing")); resSet.dimX = std::stoi(elem->Attribute("dimX")); } if (elemName == "BMode") { processSet.DoBmode = std::stoi(elem->Attribute("do")); if (!processSet.DoBmode) continue; std::string methodStr = elem->Attribute("method"); if (methodStr == "EnvelopeDetection") bmodeSet.method = mitk::PhotoacousticFilterService::BModeMethod::EnvelopeDetection; else if (elem->Attribute("method") == "Abs") bmodeSet.method = mitk::PhotoacousticFilterService::BModeMethod::Abs; else mitkThrow() << "BMode method incorrectly set in configuration file"; bmodeSet.UseLogFilter = (bool)std::stoi(elem->Attribute("useLogFilter")); } } } int main(int argc, char * argv[]) { auto input = parseInput(argc, argv); mitk::BeamformingSettings::Pointer bfSettings; BandpassSettings bandpassSettings{5,10,1,1,1540}; BModeSettings bmodeSettings{ mitk::PhotoacousticFilterService::BModeMethod::EnvelopeDetection, false }; CropSettings cropSettings{ 0,0,0,0,0,0 }; ResampleSettings resSettings{ 0.15 }; - ProcessSettings processSettings{ true, false, false }; + ProcessSettings processSettings{ false, false, false, false, false }; MITK_INFO << "Parsing settings XML..."; try { ParseXML(input.settingsFile, input, &bfSettings, bandpassSettings, cropSettings, resSettings, bmodeSettings, processSettings); } catch (mitk::Exception e) { MITK_INFO << e; return -1; } MITK_INFO << "Parsing settings XML...[Done]"; MITK_INFO(input.verbose) << "Beamforming input image..."; mitk::Image::Pointer inputImage = input.inputImage; if (!(inputImage->GetPixelType().GetTypeAsString() == "scalar (float)" || inputImage->GetPixelType().GetTypeAsString() == " (float)")) { // we need a float image, so cast it here MITK_INFO(input.verbose) << "Casting input image to float..."; mitk::CastToFloatImageFilter::Pointer castFilter = mitk::CastToFloatImageFilter::New(); castFilter->SetInput(inputImage); castFilter->Update(); inputImage = castFilter->GetOutput(); MITK_INFO << inputImage->GetPixelType().GetPixelTypeAsString(); MITK_INFO(input.verbose) << "Casting input image to float...[Done]"; } mitk::PhotoacousticFilterService::Pointer m_FilterService = mitk::PhotoacousticFilterService::New(); mitk::Image::Pointer output = inputImage; if (processSettings.DoBeamforming) { MITK_INFO(input.verbose) << "Beamforming input image..."; output = m_FilterService->ApplyBeamforming(output, bfSettings); MITK_INFO(input.verbose) << "Beamforming input image...[Done]"; } + if (processSettings.DoCropping) + { + int err; + MITK_INFO(input.verbose) << "Applying Crop filter to image..."; + output = m_FilterService->ApplyCropping(output, + cropSettings.above, cropSettings.below, cropSettings.right, cropSettings.left, cropSettings.zStart, cropSettings.zEnd, &err); + MITK_INFO(input.verbose) << "Applying Crop filter to image...[Done]"; + } if (processSettings.DoBandpass) { MITK_INFO(input.verbose) << "Bandpassing input image..."; output = m_FilterService->ApplyBandpassFilter(output, bandpassSettings.highPass*1e6, bandpassSettings.lowPass*1e6, bandpassSettings.alphaHigh, bandpassSettings.alphaLow, 1, bandpassSettings.speedOfSound, true); MITK_INFO(input.verbose) << "Bandpassing input image...[Done]"; } - if (processSettings.DoCropping) + if (processSettings.DoBmode) { - int err; - MITK_INFO(input.verbose) << "Applying Crop filter to image..."; - output = m_FilterService->ApplyCropping(output, - cropSettings.above, cropSettings.below, cropSettings.right, cropSettings.left, cropSettings.zStart, cropSettings.zEnd, &err); - MITK_INFO(input.verbose) << "Applying Crop filter to image...[Done]"; + MITK_INFO(input.verbose) << "Applying BModeFilter to image..."; + output = m_FilterService->ApplyBmodeFilter(output, bmodeSettings.method, bmodeSettings.UseLogFilter); + MITK_INFO(input.verbose) << "Applying BModeFilter to image...[Done]"; } if (processSettings.DoResampling) { - double spacing[3] = {output->GetGeometry()->GetSpacing()[0]*output->GetDimension(0)/resSettings.dimX, resSettings.spacing, output->GetGeometry()->GetSpacing()[2]}; + double spacing[3] = { output->GetGeometry()->GetSpacing()[0] * output->GetDimension(0) / resSettings.dimX, resSettings.spacing, output->GetGeometry()->GetSpacing()[2] }; MITK_INFO(input.verbose) << "Applying Resample filter to image..."; output = m_FilterService->ApplyResampling(output, spacing); if (output->GetDimension(0) != resSettings.dimX) { - double dim[3] = {(double)resSettings.dimX, (double)output->GetDimension(1), (double)output->GetDimension(2)}; + double dim[3] = { (double)resSettings.dimX, (double)output->GetDimension(1), (double)output->GetDimension(2) }; output = m_FilterService->ApplyResamplingToDim(output, dim); } MITK_INFO(input.verbose) << "Applying Resample filter to image...[Done]"; } - if (processSettings.DoBmode) - { - MITK_INFO(input.verbose) << "Applying BModeFilter to image..."; - output = m_FilterService->ApplyBmodeFilter(output, bmodeSettings.method, bmodeSettings.UseLogFilter); - MITK_INFO(input.verbose) << "Applying BModeFilter to image...[Done]"; - } MITK_INFO(input.verbose) << "Saving image..."; mitk::IOUtil::Save(output, input.outputFilename); MITK_INFO(input.verbose) << "Saving image...[Done]"; MITK_INFO(input.verbose) << "Beamforming input image...[Done]"; } diff --git a/Modules/PhotoacousticsAlgorithms/MitkPABeamformingTool/name.options.xml b/Modules/PhotoacousticsAlgorithms/MitkPABeamformingTool/name.options.xml new file mode 100644 index 0000000000..ec367585b5 --- /dev/null +++ b/Modules/PhotoacousticsAlgorithms/MitkPABeamformingTool/name.options.xml @@ -0,0 +1,131 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/Modules/PhotoacousticsAlgorithms/Resources/DAS.cl b/Modules/PhotoacousticsAlgorithms/Resources/DAS.cl index 2540ec58db..56492ff8fb 100644 --- a/Modules/PhotoacousticsAlgorithms/Resources/DAS.cl +++ b/Modules/PhotoacousticsAlgorithms/Resources/DAS.cl @@ -1,64 +1,129 @@ /*=================================================================== 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 ckDAS( __global float* dSource, // input image __global float* dDest, // output buffer __global unsigned short* usedLines, __global unsigned short* delays, __constant float* apodArray, unsigned short apodArraySize, unsigned int inputL, unsigned int inputS, unsigned int Slices, unsigned int outputL, unsigned int outputS // parameters ) { // get thread identifier unsigned int globalPosX = get_global_id(0); unsigned int globalPosY = get_global_id(1); unsigned int globalPosZ = get_global_id(2); // terminate non-valid threads if ( globalPosX < outputL && globalPosY < outputS && globalPosZ < Slices ) { float l_i = (float)globalPosX / (float)outputL * (float)inputL; unsigned short curUsedLines = usedLines[globalPosY * 3 * outputL + 3 * globalPosX]; unsigned short minLine = usedLines[globalPosY * 3 * outputL + 3 * globalPosX + 1]; unsigned short maxLine = usedLines[globalPosY * 3 *outputL + 3 * globalPosX + 2]; float apod_mult = (float)apodArraySize / (float)curUsedLines; unsigned short Delay = 0; float output = 0; float mult = 0; for (short l_s = minLine; l_s < maxLine; ++l_s) { Delay = delays[globalPosY * (outputL / 2) + (int)(fabs(l_s - l_i)/(float)inputL * (float)outputL)]; if (Delay < inputS && Delay >= 0) { output += apodArray[(int)((l_s - minLine)*apod_mult)] * dSource[(int)(globalPosZ * inputL * inputS + Delay * inputL + l_s)]; } else --curUsedLines; } dDest[ globalPosZ * outputL * outputS + globalPosY * outputL + globalPosX ] = output / (float)curUsedLines; } +} + +__kernel void ckDAS_g( + __global float* dSource, // input image + __global float* dDest, // output buffer + __global float* elementHeights, + __global float* elementPositions, + __constant float* apodArray, + unsigned short apodArraySize, + unsigned int inputL, + unsigned int inputS, + int Slices, + int outputL, + int outputS, + float angle, + float probeRadius, + float totalSamples_i, + float horizontalExtent, + float mult, + char isPAImage, + __global unsigned short* usedLines // parameters +) +{ + // get thread identifier + int globalPosX = get_global_id(0); + int globalPosY = get_global_id(1); + int globalPosZ = get_global_id(2); + + // terminate non-valid threads + if ( globalPosX < outputL && globalPosY < outputS && globalPosZ < Slices ) + { + int AddSample = 0; + float l_i = 0; + float l_p = 0; + float s_i = 0; + + float apod_mult = 1; + + float output = 0; + + l_i = (float)globalPosX / outputL * inputL; + l_p = (float)globalPosX / outputL * horizontalExtent; + s_i = (float)globalPosY / outputS * totalSamples_i; + + unsigned short curUsedLines = usedLines[globalPosY * 3 * outputL + 3 * globalPosX]; + unsigned short minLine = usedLines[globalPosY * 3 * outputL + 3 * globalPosX + 1]; + unsigned short maxLine = usedLines[globalPosY * 3 *outputL + 3 * globalPosX + 2]; + + apod_mult = (float)apodArraySize / curUsedLines; + + for (int l_s = minLine; l_s < maxLine; ++l_s) + { + AddSample = (int)sqrt( + pow(s_i-elementHeights[l_s]*mult, 2) + + + pow(mult * (l_p - elementPositions[l_s]), 2) + ) + (1 - isPAImage)*s_i; + if (AddSample < inputS && AddSample >= 0) + output += dSource[l_s + AddSample*inputL] * + apodArray[(int)((l_s - minLine)*apod_mult)]; + else + --curUsedLines; + } + dDest[ globalPosZ * outputL * outputS + globalPosY * outputL + globalPosX ] = output / curUsedLines; + } } \ No newline at end of file diff --git a/Modules/PhotoacousticsAlgorithms/Resources/DMAS.cl b/Modules/PhotoacousticsAlgorithms/Resources/DMAS.cl index 89cc58bf12..0a925ddadc 100644 --- a/Modules/PhotoacousticsAlgorithms/Resources/DMAS.cl +++ b/Modules/PhotoacousticsAlgorithms/Resources/DMAS.cl @@ -1,85 +1,170 @@ /*=================================================================== 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 ckDMAS( __global float* dSource, // input image __global float* dDest, // output buffer __global unsigned short* usedLines, __global unsigned short* AddSamples, __constant float* apodArray, unsigned short apodArraySize, unsigned int inputL, unsigned int inputS, unsigned int Slices, unsigned int outputL, unsigned int outputS // parameters ) { // get thread identifier unsigned int globalPosX = get_global_id(0); unsigned int globalPosY = get_global_id(1); unsigned int globalPosZ = get_global_id(2); // terminate non-valid threads if ( globalPosX < outputL && globalPosY < outputS && globalPosZ < Slices ) { float l_i = (float)globalPosX / (float)outputL * (float)inputL; unsigned short curUsedLines = usedLines[globalPosY * 3 * outputL + 3 * globalPosX]; unsigned short minLine = usedLines[globalPosY * 3 * outputL + 3 * globalPosX + 1]; unsigned short maxLine = usedLines[globalPosY * 3 *outputL + 3 * globalPosX + 2]; float apod_mult = (float)apodArraySize / (float)curUsedLines; unsigned short Delay1 = 0; unsigned short Delay2 = 0; float output = 0; float mult = 0; float s_1 = 0; float s_2 = 0; float apod_1 = 0; for (short l_s1 = minLine; l_s1 < maxLine; ++l_s1) { Delay1 = AddSamples[globalPosY * (outputL / 2) + (int)(fabs(l_s1 - l_i)/(float)inputL * (float)outputL)]; if (Delay1 < inputS && Delay1 >= 0) { s_1 = dSource[(int)(globalPosZ * inputL * inputS + Delay1 * inputL + l_s1)]; apod_1 = apodArray[(int)((l_s1 - minLine)*apod_mult)]; for (short l_s2 = l_s1 + 1; l_s2 < maxLine; ++l_s2) { Delay2 = AddSamples[globalPosY * (outputL / 2) + (int)(fabs(l_s2 - l_i)/(float)inputL * (float)outputL)]; if (Delay2 < inputS && Delay2 >= 0) { s_2 = dSource[(int)(globalPosZ * inputL * inputS + Delay2 * inputL + l_s2)]; - mult = apodArray[(int)((l_s2 - minLine)*apod_mult)] * s_2 - * apod_1 * s_1; + mult = apodArray[(int)((l_s2 - minLine)*apod_mult)] * s_2 * apod_1 * s_1; - output += sqrt(fabs(mult)) * ((mult > 0) - (mult < 0)); + output += sqrt(fabs(mult)) * sign(mult); } } } else --curUsedLines; } dDest[ globalPosZ * outputL * outputS + globalPosY * outputL + globalPosX ] = output / (float)(curUsedLines * curUsedLines - (curUsedLines - 1)); } } + +__kernel void ckDMAS_g( + __global float* dSource, // input image + __global float* dDest, // output buffer + __global float* elementHeights, + __global float* elementPositions, + __constant float* apodArray, + unsigned short apodArraySize, + unsigned int inputL, + unsigned int inputS, + int Slices, + int outputL, + int outputS, + float angle, + float probeRadius, + float totalSamples_i, + float horizontalExtent, + float mult, + char isPAImage, + __global unsigned short* usedLines // parameters +) +{ + // get thread identifier + int globalPosX = get_global_id(0); + int globalPosY = get_global_id(1); + int globalPosZ = get_global_id(2); + + // terminate non-valid threads + if ( globalPosX < outputL && globalPosY < outputS && globalPosZ < Slices ) + { + int AddSample1 = 0; + int AddSample2 = 0; + + float output = 0; + + float s_1 = 0; + float s_2 = 0; + float apod_1 = 0; + + float l_i = (float)globalPosX / outputL * inputL; + float l_p = (float)globalPosX / outputL * horizontalExtent; + float s_i = (float)globalPosY / outputS * totalSamples_i; + + unsigned short curUsedLines = usedLines[globalPosY * 3 * outputL + 3 * globalPosX]; + unsigned short minLine = usedLines[globalPosY * 3 * outputL + 3 * globalPosX + 1]; + unsigned short maxLine = usedLines[globalPosY * 3 *outputL + 3 * globalPosX + 2]; + + float apod_mult = (float)apodArraySize / curUsedLines; + + float multiplication = 0; + + for (int l_s1 = minLine; l_s1 < maxLine; ++l_s1) + { + AddSample1 = (int)sqrt( + pow(s_i-elementHeights[l_s1]*mult, 2) + + + pow(mult * (l_p - elementPositions[l_s1]), 2) + ) + (1 - isPAImage)*s_i; + + if (AddSample1 < inputS && AddSample1 >= 0) + { + s_1 = dSource[(int)(globalPosZ * inputL * inputS + AddSample1 * inputL + l_s1)]; + apod_1 = apodArray[(int)((l_s1 - minLine)*apod_mult)]; + + for (int l_s2 = minLine; l_s2 < maxLine; ++l_s2) + { + AddSample2 = (int)sqrt( + pow(s_i-elementHeights[l_s2]*mult, 2) + + + pow(mult * (l_p - elementPositions[l_s2]), 2) + ) + (1 - isPAImage)*s_i; + if (AddSample2 < inputS && AddSample2 >= 0) + { + s_2 = dSource[(int)(globalPosZ * inputL * inputS + AddSample2 * inputL + l_s2)]; + multiplication = apodArray[(int)((l_s2 - minLine)*apod_mult)] * s_2 * apod_1 * s_1; + + output += sqrt(fabs(multiplication)) * sign(multiplication); + } + } + } + else + --curUsedLines; + } + dDest[ globalPosZ * outputL * outputS + globalPosY * outputL + globalPosX ] = output / (float)(pow((float)curUsedLines, 2) - (curUsedLines - 1)); + } +} \ No newline at end of file diff --git a/Modules/PhotoacousticsAlgorithms/Resources/DelayCalculation.cl b/Modules/PhotoacousticsAlgorithms/Resources/DelayCalculation.cl index c0ee5224e4..b71208e6e8 100644 --- a/Modules/PhotoacousticsAlgorithms/Resources/DelayCalculation.cl +++ b/Modules/PhotoacousticsAlgorithms/Resources/DelayCalculation.cl @@ -1,70 +1,44 @@ /*=================================================================== 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 - ) -{ - 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 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, - float percentOfImageReconstructed // parameters + float totalSamples_i // 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 - + 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/Resources/UsedLinesCalculation.cl b/Modules/PhotoacousticsAlgorithms/Resources/UsedLinesCalculation.cl index aa63d1562e..132c7d88fd 100644 --- a/Modules/PhotoacousticsAlgorithms/Resources/UsedLinesCalculation.cl +++ b/Modules/PhotoacousticsAlgorithms/Resources/UsedLinesCalculation.cl @@ -1,52 +1,170 @@ /*=================================================================== 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, - char isPAImage, - float percentOfImageReconstructed // parameters + float totalSamples_i // parameters ) { // get thread identifier unsigned int globalPosX = get_global_id(0); unsigned int globalPosY = get_global_id(1); - //unsigned short outputS = get_global_size(1); - //unsigned short outputL = get_global_size(0); - // terminate non-valid threads if ( globalPosX < outputL && globalPosY < outputS) { float l_i = (float)globalPosX / outputL * inputL; - float s_i = (float)globalPosY / (float)outputS * (float)inputS / (2-isPAImage) * percentOfImageReconstructed; + 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 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); + + // 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; + + 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) + { + 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; + 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) + { + 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/Resources/sDMAS.cl b/Modules/PhotoacousticsAlgorithms/Resources/sDMAS.cl index b362047846..bb6096b3c6 100644 --- a/Modules/PhotoacousticsAlgorithms/Resources/sDMAS.cl +++ b/Modules/PhotoacousticsAlgorithms/Resources/sDMAS.cl @@ -1,87 +1,174 @@ /*=================================================================== 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 cksDMAS( __global float* dSource, // input image __global float* dDest, // output buffer __global unsigned short* usedLines, __global unsigned short* AddSamples, __constant float* apodArray, unsigned short apodArraySize, unsigned int inputL, unsigned int inputS, unsigned int Slices, unsigned int outputL, unsigned int outputS // parameters ) { // get thread identifier unsigned int globalPosX = get_global_id(0); unsigned int globalPosY = get_global_id(1); unsigned int globalPosZ = get_global_id(2); // terminate non-valid threads if ( globalPosX < outputL && globalPosY < outputS && globalPosZ < Slices ) { float l_i = (float)globalPosX / (float)outputL * (float)inputL; unsigned short curUsedLines = usedLines[globalPosY * 3 * outputL + 3 * globalPosX]; unsigned short minLine = usedLines[globalPosY * 3 * outputL + 3 * globalPosX + 1]; unsigned short maxLine = usedLines[globalPosY * 3 *outputL + 3 * globalPosX + 2]; float apod_mult = (float)apodArraySize / (float)curUsedLines; unsigned short Delay1 = 0; unsigned short Delay2 = 0; float output = 0; float mult = 0; float s_1 = 0; float s_2 = 0; - float sign = 0; + float dSign = 0; float apod_1 = 0; for (short l_s1 = minLine; l_s1 < maxLine; ++l_s1) { Delay1 = AddSamples[globalPosY * (outputL / 2) + (int)(fabs(l_s1 - l_i)/(float)inputL * (float)outputL)]; if (Delay1 < inputS && Delay1 >= 0) { s_1 = dSource[(int)(globalPosZ * inputL * inputS + Delay1 * inputL + l_s1)]; apod_1 = apodArray[(int)((l_s1 - minLine)*apod_mult)]; - sign += s_1; + dSign += s_1; for (short l_s2 = l_s1 + 1; l_s2 < maxLine; ++l_s2) { Delay2 = AddSamples[globalPosY * (outputL / 2) + (int)(fabs(l_s2 - l_i)/(float)inputL * (float)outputL)]; if (Delay2 < inputS && Delay2 >= 0) { s_2 = dSource[(int)(globalPosZ * inputL * inputS + Delay2 * inputL + l_s2)]; - mult = apodArray[(int)((l_s2 - minLine)*apod_mult)] * s_2 - * apod_1 * s_1; + mult = apodArray[(int)((l_s2 - minLine)*apod_mult)] * s_2 * apod_1 * s_1; output += sqrt(fabs(mult)) * ((mult > 0) - (mult < 0)); } } } else --curUsedLines; } - dDest[ globalPosZ * outputL * outputS + globalPosY * outputL + globalPosX ] = output / (float)(curUsedLines * curUsedLines - (curUsedLines - 1)) * ((sign > 0) - (sign < 0)); + dDest[ globalPosZ * outputL * outputS + globalPosY * outputL + globalPosX ] = output / (float)(curUsedLines * curUsedLines - (curUsedLines - 1)) * sign(dSign); } } + +__kernel void cksDMAS_g( + __global float* dSource, // input image + __global float* dDest, // output buffer + __global float* elementHeights, + __global float* elementPositions, + __constant float* apodArray, + unsigned short apodArraySize, + unsigned int inputL, + unsigned int inputS, + int Slices, + int outputL, + int outputS, + float angle, + float probeRadius, + float totalSamples_i, + float horizontalExtent, + float mult, + char isPAImage, + __global unsigned short* usedLines // parameters +) +{ + // get thread identifier + int globalPosX = get_global_id(0); + int globalPosY = get_global_id(1); + int globalPosZ = get_global_id(2); + + // terminate non-valid threads + if ( globalPosX < outputL && globalPosY < outputS && globalPosZ < Slices ) + { + int AddSample1 = 0; + int AddSample2 = 0; + + float output = 0; + + float s_1 = 0; + float s_2 = 0; + float apod_1 = 0; + + float l_i = (float)globalPosX / outputL * inputL; + float l_p = (float)globalPosX / outputL * horizontalExtent; + float s_i = (float)globalPosY / outputS * totalSamples_i; + + unsigned short curUsedLines = usedLines[globalPosY * 3 * outputL + 3 * globalPosX]; + unsigned short minLine = usedLines[globalPosY * 3 * outputL + 3 * globalPosX + 1]; + unsigned short maxLine = usedLines[globalPosY * 3 *outputL + 3 * globalPosX + 2]; + + float apod_mult = (float)apodArraySize / curUsedLines; + float dSign = 0; + + float multiplication = 0; + + for (int l_s1 = minLine; l_s1 < maxLine; ++l_s1) + { + AddSample1 = (int)sqrt( + pow(s_i-elementHeights[l_s1]*mult, 2) + + + pow(mult * (l_p - elementPositions[l_s1]), 2) + ) + (1 - isPAImage)*s_i; + + if (AddSample1 < inputS && AddSample1 >= 0) + { + s_1 = dSource[(int)(globalPosZ * inputL * inputS + AddSample1 * inputL + l_s1)]; + apod_1 = apodArray[(int)((l_s1 - minLine)*apod_mult)]; + dSign += s_1; + + for (int l_s2 = minLine; l_s2 < maxLine; ++l_s2) + { + AddSample2 = (int)sqrt( + pow(s_i-elementHeights[l_s2]*mult, 2) + + + pow(mult * (l_p - elementPositions[l_s2]), 2) + ) + (1 - isPAImage)*s_i; + if (AddSample2 < inputS && AddSample2 >= 0) + { + s_2 = dSource[(int)(globalPosZ * inputL * inputS + AddSample2 * inputL + l_s2)]; + multiplication = apodArray[(int)((l_s2 - minLine)*apod_mult)] * s_2 * apod_1 * s_1; + + output += sqrt(fabs(multiplication)) * sign(multiplication); + } + } + } + else + --curUsedLines; + } + dDest[ globalPosZ * outputL * outputS + globalPosY * outputL + globalPosX ] = output / (float)(pow((float)curUsedLines, 2) - (curUsedLines - 1)) * sign(dSign); + } +} \ No newline at end of file diff --git a/Modules/PhotoacousticsAlgorithms/include/OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.h b/Modules/PhotoacousticsAlgorithms/include/OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.h index b60cd2315d..bd9bf91913 100644 --- a/Modules/PhotoacousticsAlgorithms/include/OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.h +++ b/Modules/PhotoacousticsAlgorithms/include/OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.h @@ -1,145 +1,147 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef _MITKPHOTOACOUSTICSOCLBEAMFORMER_H_ #define _MITKPHOTOACOUSTICSOCLBEAMFORMER_H_ #include #if defined(PHOTOACOUSTICS_USE_GPU) || DOXYGEN #include "mitkOclDataSetToDataSetFilter.h" #include "mitkBeamformingSettings.h" #include "mitkPhotoacousticOCLDelayCalculation.h" #include "mitkPhotoacousticOCLUsedLinesCalculation.h" #include namespace mitk { /*! * \brief Class implementing a mitk::OclDataSetToDataSetFilter for beamforming on GPU * * The class must be given a configuration class instance of mitk::BeamformingSettings for beamforming parameters through mitk::PhotoacousticOCLBeamformingFilter::SetConfig(BeamformingSettings settings) * Additional configuration of the apodisation function is needed. */ class PhotoacousticOCLBeamformingFilter : public OclDataSetToDataSetFilter, public itk::Object { public: mitkClassMacroItkParent(PhotoacousticOCLBeamformingFilter, itk::Object); mitkNewMacro1Param(Self, BeamformingSettings::Pointer); /** * @brief SetInput Set the input data through an image. Arbitrary images are supported */ void SetInput(Image::Pointer image); /** * brief SetInput Manually set the input data while providing 3 dimensions and memory size of the input data (Bytes per element). */ void SetInput(void* data, unsigned int* dimensions, unsigned int BpE); /** * @brief GetOutput Get a pointer to the processed data. The standard datatype is float. */ void* GetOutput(); /** * @brief GetOutputAsImage Returns an mitk::Image constructed from the processed data */ mitk::Image::Pointer GetOutputAsImage(); /** \brief Update the filter */ void Update(); /** \brief Set the Apodisation function to apply when beamforming */ void SetApodisation(const float* apodisation, unsigned short apodArraySize) { m_ApodArraySize = apodArraySize; m_Apodisation = apodisation; } protected: PhotoacousticOCLBeamformingFilter(BeamformingSettings::Pointer settings); virtual ~PhotoacousticOCLBeamformingFilter(); /** \brief Initialize the filter */ bool Initialize(); /** \brief Updated the used data for beamforming depending on whether the configuration has significantly changed */ void UpdateDataBuffers(); /** \brief Execute the filter */ void Execute(); mitk::PixelType GetOutputType() { return mitk::MakeScalarPixelType(); } int GetBytesPerElem() { return sizeof(float); } virtual us::Module* GetModule(); private: /** The OpenCL kernel for the filter */ cl_kernel m_PixelCalculation; unsigned int m_OutputDim[3]; const float* m_Apodisation; unsigned short m_ApodArraySize; unsigned int m_inputSlices; unsigned short m_PAImage; BeamformingSettings::Pointer m_Conf; mitk::Image::Pointer m_InputImage; size_t m_ChunkSize[3]; mitk::OCLUsedLinesCalculation::Pointer m_UsedLinesCalculation; mitk::OCLDelayCalculation::Pointer m_DelayCalculation; - cl_mem m_ApodizationBuffer; - cl_mem m_MemoryLocationsBuffer; + cl_mem m_ApodizationBuffer; cl_mem m_DelaysBuffer; cl_mem m_UsedLinesBuffer; + cl_mem m_ElementHeightsBuffer; + cl_mem m_ElementPositionsBuffer; + }; } #else namespace mitk { class PhotoacousticOCLBeamformingFilter : public itk::Object { public: mitkClassMacroItkParent(mitk::PhotoacousticOCLBeamformingFilter, itk::Object); itkNewMacro(Self); protected: /** Constructor */ PhotoacousticOCLBeamformingFilter() {} /** Destructor */ ~PhotoacousticOCLBeamformingFilter() override {} }; } #endif #endif diff --git a/Modules/PhotoacousticsAlgorithms/include/OpenCLFilter/mitkPhotoacousticOCLUsedLinesCalculation.h b/Modules/PhotoacousticsAlgorithms/include/OpenCLFilter/mitkPhotoacousticOCLUsedLinesCalculation.h index 4394764e23..52021ea7c5 100644 --- a/Modules/PhotoacousticsAlgorithms/include/OpenCLFilter/mitkPhotoacousticOCLUsedLinesCalculation.h +++ b/Modules/PhotoacousticsAlgorithms/include/OpenCLFilter/mitkPhotoacousticOCLUsedLinesCalculation.h @@ -1,79 +1,84 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef _MITKPHOTOACOUSTICSOCLUSEDLINESCALCULATION_H_ #define _MITKPHOTOACOUSTICSOCLUSEDLINESCALCULATION_H_ #if defined(PHOTOACOUSTICS_USE_GPU) || DOXYGEN #include "mitkOclDataSetToDataSetFilter.h" #include #include "mitkBeamformingSettings.h" namespace mitk { /*! * \brief Class implementing a mitk::OclDataSetToDataSetFilter to calculate which lines each sample should use when beamforming. * * The class must be given a configuration class instance of mitk::BeamformingSettings for beamforming parameters through mitk::OCLDelayCalculation::SetConfig(BeamformingSettings conf) */ class OCLUsedLinesCalculation : public OclDataSetToDataSetFilter, public itk::Object { public: mitkClassMacroItkParent(OCLUsedLinesCalculation, itk::Object); mitkNewMacro1Param(Self, mitk::BeamformingSettings::Pointer); void Update(); + void SetElementHeightsBuffer(cl_mem elementHeightsBuffer); + void SetElementPositionsBuffer(cl_mem elementPositionsBuffer); + protected: /** Constructor */ OCLUsedLinesCalculation(mitk::BeamformingSettings::Pointer settings); /** Destructor */ virtual ~OCLUsedLinesCalculation(); /** Initialize the filter */ bool Initialize(); void Execute(); mitk::PixelType GetOutputType() { return mitk::MakeScalarPixelType(); } int GetBytesPerElem() { return sizeof(unsigned short); } virtual us::Module* GetModule(); int m_sizeThis; private: /** The OpenCL kernel for the filter */ cl_kernel m_PixelCalculation; BeamformingSettings::Pointer m_Conf; float m_part; size_t m_ChunkSize[3]; + cl_mem m_ElementHeightsBuffer; + cl_mem m_ElementPositionsBuffer; }; } #endif #endif diff --git a/Modules/PhotoacousticsAlgorithms/include/mitkBeamformingFilter.h b/Modules/PhotoacousticsAlgorithms/include/mitkBeamformingFilter.h index 8724808d36..6ee498e470 100644 --- a/Modules/PhotoacousticsAlgorithms/include/mitkBeamformingFilter.h +++ b/Modules/PhotoacousticsAlgorithms/include/mitkBeamformingFilter.h @@ -1,89 +1,90 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef MITK_PHOTOACOUSTICS_BEAMFORMING_FILTER #define MITK_PHOTOACOUSTICS_BEAMFORMING_FILTER #include "mitkImageToImageFilter.h" #include #include "./OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.h" #include "mitkBeamformingSettings.h" #include "mitkBeamformingUtils.h" #include "MitkPhotoacousticsAlgorithmsExports.h" namespace mitk { /*! * \brief Class implementing an mitk::ImageToImageFilter for beamforming on both CPU and GPU * * The class must be given a configuration class instance of mitk::BeamformingSettings for beamforming parameters through mitk::BeamformingFilter::Configure(BeamformingSettings settings) * Whether the GPU is used can be set in the configuration. * For significant problems or important messages a string is written, which can be accessed via GetMessageString(). */ class MITKPHOTOACOUSTICSALGORITHMS_EXPORT BeamformingFilter : public ImageToImageFilter { public: mitkClassMacro(BeamformingFilter, ImageToImageFilter); mitkNewMacro1Param(Self, mitk::BeamformingSettings::Pointer); itkCloneMacro(Self); /** \brief Sets a callback for progress checking * * An std::function can be set, through which progress of the currently updating filter is reported. * The integer argument is a number between 0 an 100 to indicate how far completion has been achieved, the std::string argument indicates what the filter is currently doing. */ void SetProgressHandle(std::function progressHandle); protected: BeamformingFilter(mitk::BeamformingSettings::Pointer settings); ~BeamformingFilter() override; void GenerateInputRequestedRegion() override; void GenerateOutputInformation() override; void GenerateData() override; //##Description //## @brief Time when Header was last initialized itk::TimeStamp m_TimeOfHeaderInitialization; /** \brief The std::function, through which progress of the currently updating filter is reported. */ std::function m_ProgressHandle; float* m_OutputData; float* m_InputData; float* m_InputDataPuffer; + unsigned short* m_MinMaxLines; /** \brief Current configuration set */ BeamformingSettings::Pointer m_Conf; /** * The size of the apodization array when it last changed. */ int m_LastApodizationArraySize; /** \brief Pointer to the GPU beamforming filter class; for performance reasons the filter is initialized within the constructor and kept for all later computations. */ mitk::PhotoacousticOCLBeamformingFilter::Pointer m_BeamformingOclFilter; }; } // namespace mitk #endif //MITK_PHOTOACOUSTICS_BEAMFORMING_FILTER diff --git a/Modules/PhotoacousticsAlgorithms/include/mitkBeamformingSettings.h b/Modules/PhotoacousticsAlgorithms/include/mitkBeamformingSettings.h index 1b36181e79..b50654b3f0 100644 --- a/Modules/PhotoacousticsAlgorithms/include/mitkBeamformingSettings.h +++ b/Modules/PhotoacousticsAlgorithms/include/mitkBeamformingSettings.h @@ -1,218 +1,248 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef MITK_BEAMFORMING_SETTINGS #define MITK_BEAMFORMING_SETTINGS #include #include #include #include namespace mitk { /*! * \brief Class holding the configuration data for the beamforming filters mitk::BeamformingFilter and mitk::PhotoacousticOCLBeamformingFilter * * A detailed description can be seen below. All parameters should be set manually for successfull beamforming. */ class MITKPHOTOACOUSTICSALGORITHMS_EXPORT BeamformingSettings : public itk::Object { public: mitkClassMacroItkParent(BeamformingSettings, itk::Object); itkCloneMacro(Self); - /** \brief Available delay calculation methods: - * - Spherical delay for best results. - * - DEPRECATED quadratic Taylor approximation for slightly faster results with hardly any quality loss. - */ - enum DelayCalc { QuadApprox, Spherical }; - /** \brief Available apodization functions: * - Hamming function. * - Von-Hann function. * - Box function. */ enum Apodization { Hamm, Hann, Box }; /** \brief Available beamforming algorithms: * - DAS (Delay and sum). * - DMAS (Delay multiply and sum). */ enum BeamformingAlgorithm { DMAS, DAS, sDMAS }; + /** \brief Available geometries for Probes: + * - Linear + * - Concave + */ + enum ProbeGeometry { Linear, Concave}; + itkGetConstMacro(PitchInMeters, float); itkGetConstMacro(SpeedOfSound, float); itkGetConstMacro(TimeSpacing, float); itkGetConstMacro(Angle, float); itkGetConstMacro(IsPhotoacousticImage, bool); itkGetConstMacro(TransducerElements, unsigned int); itkGetConstMacro(SamplesPerLine, unsigned int); itkGetConstMacro(ReconstructionLines, unsigned int); itkGetConstMacro(InputDim, const unsigned int*); itkGetConstMacro(UseGPU, bool); itkGetConstMacro(GPUBatchSize, unsigned int); - itkGetConstMacro(DelayCalculationMethod, DelayCalc); itkGetConstMacro(ApodizationFunction, const float*); itkGetConstMacro(Apod, Apodization); itkGetConstMacro(ApodizationArraySize, int); itkGetConstMacro(Algorithm, BeamformingAlgorithm); itkGetConstMacro(ReconstructionDepth, float); + itkGetConstMacro(Geometry, ProbeGeometry); + itkGetConstMacro(ProbeRadius, float); + itkGetConstMacro(ElementHeights, float*); + itkGetConstMacro(ElementPositions, float*); + itkGetConstMacro(HorizontalExtent, float); /** \brief function for mitk::PhotoacousticOCLBeamformingFilter to check whether buffers need to be updated * this method only checks parameters relevant for the openCL implementation */ static bool SettingsChangedOpenCL(const BeamformingSettings::Pointer lhs, const BeamformingSettings::Pointer rhs) { return !((std::abs(lhs->GetAngle() - rhs->GetAngle()) < 0.01f) && // 0.01 degree error margin (lhs->GetApod() == rhs->GetApod()) && - (lhs->GetDelayCalculationMethod() == rhs->GetDelayCalculationMethod()) && + (lhs->GetGeometry() == rhs->GetGeometry()) && + (abs(lhs->GetProbeRadius() - rhs->GetProbeRadius()) < 0.001f) && (lhs->GetIsPhotoacousticImage() == rhs->GetIsPhotoacousticImage()) && (std::abs(lhs->GetPitchInMeters() - rhs->GetPitchInMeters()) < 0.000001f) && // 0.0001 mm error margin (lhs->GetReconstructionLines() == rhs->GetReconstructionLines()) && (lhs->GetSamplesPerLine() == rhs->GetSamplesPerLine()) && (lhs->GetReconstructionDepth() == rhs->GetReconstructionDepth()) && (std::abs(lhs->GetSpeedOfSound() - rhs->GetSpeedOfSound()) < 0.01f) && (std::abs(lhs->GetTimeSpacing() - rhs->GetTimeSpacing()) < 0.00000000001f) && //0.01 ns error margin (lhs->GetTransducerElements() == rhs->GetTransducerElements())); } static Pointer New(float pitchInMeters, float speedOfSound, float timeSpacing, float angle, bool isPhotoacousticImage, unsigned int samplesPerLine, unsigned int reconstructionLines, unsigned int* inputDim, float reconstructionDepth, bool useGPU, unsigned int GPUBatchSize, - DelayCalc delayCalculationMethod, Apodization apod, unsigned int apodizationArraySize, - BeamformingAlgorithm algorithm) + BeamformingAlgorithm algorithm, + ProbeGeometry geometry, + float probeRadius) { Pointer smartPtr = new BeamformingSettings(pitchInMeters, speedOfSound, timeSpacing, angle, isPhotoacousticImage, samplesPerLine, reconstructionLines, inputDim, reconstructionDepth, useGPU, GPUBatchSize, - delayCalculationMethod, apod, apodizationArraySize, - algorithm); + algorithm, + geometry, + probeRadius); smartPtr->UnRegister(); return smartPtr; } + unsigned short* GetMinMaxLines(); + protected: /** */ BeamformingSettings(float pitchInMeters, float speedOfSound, float timeSpacing, float angle, bool isPhotoacousticImage, unsigned int samplesPerLine, unsigned int reconstructionLines, unsigned int* inputDim, float reconstructionDepth, bool useGPU, unsigned int GPUBatchSize, - DelayCalc delayCalculationMethod, Apodization apod, unsigned int apodizationArraySize, - BeamformingAlgorithm algorithm + BeamformingAlgorithm algorithm, + ProbeGeometry geometry, + float probeRadius ); ~BeamformingSettings(); /** \brief Pitch of the used transducer in [m]. */ float m_PitchInMeters; /** \brief Speed of sound in the used medium in [m/s]. */ float m_SpeedOfSound; /** \brief The time spacing of the input image */ float m_TimeSpacing; // [s] /** \brief The angle of the transducer elements */ float m_Angle; /** \brief Flag whether processed image is a photoacoustic image or an ultrasound image */ bool m_IsPhotoacousticImage; /** \brief How many transducer elements the used transducer had. */ unsigned int m_TransducerElements; /** \brief How many vertical samples should be used in the final image. */ unsigned int m_SamplesPerLine; /** \brief How many lines should be reconstructed in the final image. */ unsigned int m_ReconstructionLines; /** \brief Sets the dimensions of the inputImage. */ const unsigned int* m_InputDim; /** \brief The Depth up to which the filter should reconstruct the image [m] */ float m_ReconstructionDepth; /** \brief Decides whether GPU computing should be used */ bool m_UseGPU; unsigned int m_GPUBatchSize; /** \brief Sets the amount of image slices in batches when GPU is used */ - /** \brief Sets how the delays for beamforming should be calculated. - */ - DelayCalc m_DelayCalculationMethod; - const float* m_ApodizationFunction; /** \brief Sets the used apodization function. */ Apodization m_Apod; /** \brief Sets the resolution of the apodization array (must be greater than 0). */ int m_ApodizationArraySize; /** \brief Sets the used beamforming algorithm. */ BeamformingAlgorithm m_Algorithm; + + /** \brief Sets the used probe geometry + */ + ProbeGeometry m_Geometry; + + /** \brief Sets the radius of the curved probe [m] + */ + float m_ProbeRadius; + + /** + */ + float *m_ElementHeights; + + /** + */ + float *m_ElementPositions; + + /** + */ + float m_HorizontalExtent; + + /** + */ + unsigned short* m_MinMaxLines; }; } #endif //MITK_BEAMFORMING_SETTINGS diff --git a/Modules/PhotoacousticsAlgorithms/include/mitkBeamformingUtils.h b/Modules/PhotoacousticsAlgorithms/include/mitkBeamformingUtils.h index 9060284cfd..fff7126538 100644 --- a/Modules/PhotoacousticsAlgorithms/include/mitkBeamformingUtils.h +++ b/Modules/PhotoacousticsAlgorithms/include/mitkBeamformingUtils.h @@ -1,80 +1,72 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef MITK_BEAMFORMING_FILTER_UTILS #define MITK_BEAMFORMING_FILTER_UTILS #include "mitkImageToImageFilter.h" #include #include "./OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.h" #include "mitkBeamformingSettings.h" namespace mitk { /*! * \brief Class implementing util functionality for beamforming on CPU * */ class BeamformingUtils final { public: - /** \brief Function to perform beamforming on CPU for a single line, using DAS and quadratic delay - */ - static void DASQuadraticLine(float* input, float* output, float inputDim[2], float outputDim[2], const short& line, const mitk::BeamformingSettings::Pointer config); - /** \brief Function to perform beamforming on CPU for a single line, using DAS and spherical delay */ static void DASSphericalLine(float* input, float* output, float inputDim[2], float outputDim[2], const short& line, const mitk::BeamformingSettings::Pointer config); - /** \brief Function to perform beamforming on CPU for a single line, using DMAS and quadratic delay - */ - static void DMASQuadraticLine(float* input, float* output, float inputDim[2], float outputDim[2], const short& line, const mitk::BeamformingSettings::Pointer config); - /** \brief Function to perform beamforming on CPU for a single line, using DMAS and spherical delay */ static void DMASSphericalLine(float* input, float* output, float inputDim[2], float outputDim[2], const short& line, const mitk::BeamformingSettings::Pointer config); - /** \brief Function to perform beamforming on CPU for a single line, using signed DMAS and quadratic delay - */ - static void sDMASQuadraticLine(float* input, float* output, float inputDim[2], float outputDim[2], const short& line, const mitk::BeamformingSettings::Pointer config); - /** \brief Function to perform beamforming on CPU for a single line, using signed DMAS and spherical delay */ static void sDMASSphericalLine(float* input, float* output, float inputDim[2], float outputDim[2], const short& line, const mitk::BeamformingSettings::Pointer config); /** \brief Pointer holding the Von-Hann apodization window for beamforming * @param samples the resolution at which the window is created */ static float* VonHannFunction(int samples); /** \brief Function to create a Hamming apodization window * @param samples the resolution at which the window is created */ static float* HammFunction(int samples); /** \brief Function to create a Box apodization window * @param samples the resolution at which the window is created */ static float* BoxFunction(int samples); + /** \brief + */ + static unsigned short* MinMaxLines(const mitk::BeamformingSettings::Pointer config); + protected: BeamformingUtils(); ~BeamformingUtils(); }; } // namespace mitk #endif //MITK_BEAMFORMING_FILTER_UTILS diff --git a/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.cpp b/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.cpp index 9bf7ace7cc..8721ea4e24 100644 --- a/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.cpp +++ b/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.cpp @@ -1,259 +1,329 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #if defined(PHOTOACOUSTICS_USE_GPU) || DOXYGEN #include "./OpenCLFilter/mitkPhotoacousticOCLBeamformingFilter.h" #include "usServiceReference.h" mitk::PhotoacousticOCLBeamformingFilter::PhotoacousticOCLBeamformingFilter(BeamformingSettings::Pointer settings) : m_PixelCalculation(NULL), m_inputSlices(1), m_Conf(settings), m_InputImage(mitk::Image::New()), m_ApodizationBuffer(nullptr), - m_MemoryLocationsBuffer(nullptr), m_DelaysBuffer(nullptr), - m_UsedLinesBuffer(nullptr) + m_UsedLinesBuffer(nullptr), + m_ElementHeightsBuffer(nullptr), + m_ElementPositionsBuffer(nullptr) { MITK_INFO << "Instantiating OCL beamforming Filter..."; this->AddSourceFile("DAS.cl"); this->AddSourceFile("DMAS.cl"); this->AddSourceFile("sDMAS.cl"); this->m_FilterID = "OpenCLBeamformingFilter"; this->Initialize(); unsigned int dim[] = { 128, 2048, 2 }; m_InputImage->Initialize(mitk::MakeScalarPixelType(), 3, dim); m_ChunkSize[0] = 128; m_ChunkSize[1] = 128; m_ChunkSize[2] = 8; m_UsedLinesCalculation = mitk::OCLUsedLinesCalculation::New(m_Conf); m_DelayCalculation = mitk::OCLDelayCalculation::New(m_Conf); MITK_INFO << "Instantiating OCL beamforming Filter...[Done]"; } mitk::PhotoacousticOCLBeamformingFilter::~PhotoacousticOCLBeamformingFilter() { if (this->m_PixelCalculation) { clReleaseKernel(m_PixelCalculation); } if (m_ApodizationBuffer) clReleaseMemObject(m_ApodizationBuffer); + if (m_ElementHeightsBuffer) clReleaseMemObject(m_ElementHeightsBuffer); + if (m_ElementPositionsBuffer) clReleaseMemObject(m_ElementPositionsBuffer); } 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(); + us::ServiceReference ref = GetModuleContext()->GetServiceReference(); OclResourceService* resources = GetModuleContext()->GetService(ref); - cl_ulong globalMemSize = oclGetGlobalMemSize(resources->GetCurrentDevice());*/ - //Initialize the Output + cl_context gpuContext = resources->GetContext(); + //Initialize the Output try { size_t outputSize = (size_t)m_Conf->GetReconstructionLines() * (size_t)m_Conf->GetSamplesPerLine() * (size_t)m_inputSlices; m_OutputDim[0] = m_Conf->GetReconstructionLines(); m_OutputDim[1] = m_Conf->GetSamplesPerLine(); m_OutputDim[2] = m_inputSlices; this->InitExec(this->m_PixelCalculation, m_OutputDim, outputSize, sizeof(float)); } catch (const mitk::Exception& e) { MITK_ERROR << "Caught exception while initializing filter: " << e.what(); return; } - //TODO FIXME cl_int clErr = 0; MITK_DEBUG << "Updating GPU Buffers for new configuration"; // create the apodisation buffer - if (m_Apodisation == nullptr) { MITK_INFO << "No apodisation function set; Beamforming will be done without any apodisation."; m_Apodisation = new float[1]{ 1 }; m_ApodArraySize = 1; } - us::ServiceReference ref = GetModuleContext()->GetServiceReference(); - OclResourceService* resources = GetModuleContext()->GetService(ref); - cl_context gpuContext = resources->GetContext(); - if (m_ApodizationBuffer) clReleaseMemObject(m_ApodizationBuffer); this->m_ApodizationBuffer = clCreateBuffer(gpuContext, CL_MEM_READ_ONLY | CL_MEM_USE_HOST_PTR, sizeof(float) * m_ApodArraySize, const_cast(m_Apodisation), &clErr); CHECK_OCL_ERR(clErr); + if (m_ElementHeightsBuffer) clReleaseMemObject(m_ElementHeightsBuffer); + this->m_ElementHeightsBuffer = clCreateBuffer(gpuContext, CL_MEM_READ_ONLY | CL_MEM_USE_HOST_PTR, sizeof(float) * m_Conf->GetTransducerElements(), const_cast(m_Conf->GetElementHeights()), &clErr); + CHECK_OCL_ERR(clErr); + + if (m_ElementPositionsBuffer) clReleaseMemObject(m_ElementPositionsBuffer); + this->m_ElementPositionsBuffer = clCreateBuffer(gpuContext, CL_MEM_READ_ONLY | CL_MEM_USE_HOST_PTR, sizeof(float) * m_Conf->GetTransducerElements(), const_cast(m_Conf->GetElementPositions()), &clErr); + CHECK_OCL_ERR(clErr); // calculate used lines + m_UsedLinesCalculation->SetElementPositionsBuffer(m_ElementPositionsBuffer); + m_UsedLinesCalculation->SetElementHeightsBuffer(m_ElementHeightsBuffer); m_UsedLinesCalculation->Update(); m_UsedLinesBuffer = m_UsedLinesCalculation->GetGPUOutput()->GetGPUBuffer(); // calculate the Delays m_DelayCalculation->SetInputs(m_UsedLinesBuffer); m_DelayCalculation->Update(); m_DelaysBuffer = m_DelayCalculation->GetGPUOutput()->GetGPUBuffer(); - - //m_ConfOld = m_Conf; } void mitk::PhotoacousticOCLBeamformingFilter::Execute() { cl_int clErr = 0; UpdateDataBuffers(); - - unsigned int reconstructionLines = this->m_Conf->GetReconstructionLines(); - unsigned int samplesPerLine = this->m_Conf->GetSamplesPerLine(); - - clErr = clSetKernelArg(this->m_PixelCalculation, 2, sizeof(cl_mem), &(this->m_UsedLinesBuffer)); - clErr |= clSetKernelArg(this->m_PixelCalculation, 3, sizeof(cl_mem), &(this->m_DelaysBuffer)); - clErr |= clSetKernelArg(this->m_PixelCalculation, 4, sizeof(cl_mem), &(this->m_ApodizationBuffer)); - clErr |= clSetKernelArg(this->m_PixelCalculation, 5, sizeof(cl_ushort), &(this->m_ApodArraySize)); - clErr |= clSetKernelArg(this->m_PixelCalculation, 6, sizeof(cl_uint), &(this->m_Conf->GetInputDim()[0])); - clErr |= clSetKernelArg(this->m_PixelCalculation, 7, sizeof(cl_uint), &(this->m_Conf->GetInputDim()[1])); - clErr |= clSetKernelArg(this->m_PixelCalculation, 8, sizeof(cl_uint), &(m_inputSlices)); - clErr |= clSetKernelArg(this->m_PixelCalculation, 9, sizeof(cl_uint), &(reconstructionLines)); - clErr |= clSetKernelArg(this->m_PixelCalculation, 10, sizeof(cl_uint), &(samplesPerLine)); - + if (m_Conf->GetGeometry() == mitk::BeamformingSettings::ProbeGeometry::Linear) + { + unsigned int reconstructionLines = this->m_Conf->GetReconstructionLines(); + unsigned int samplesPerLine = this->m_Conf->GetSamplesPerLine(); + + clErr = clSetKernelArg(this->m_PixelCalculation, 2, sizeof(cl_mem), &(this->m_UsedLinesBuffer)); + clErr |= clSetKernelArg(this->m_PixelCalculation, 3, sizeof(cl_mem), &(this->m_DelaysBuffer)); + clErr |= clSetKernelArg(this->m_PixelCalculation, 4, sizeof(cl_mem), &(this->m_ApodizationBuffer)); + clErr |= clSetKernelArg(this->m_PixelCalculation, 5, sizeof(cl_ushort), &(this->m_ApodArraySize)); + clErr |= clSetKernelArg(this->m_PixelCalculation, 6, sizeof(cl_uint), &(this->m_Conf->GetInputDim()[0])); + clErr |= clSetKernelArg(this->m_PixelCalculation, 7, sizeof(cl_uint), &(this->m_Conf->GetInputDim()[1])); + clErr |= clSetKernelArg(this->m_PixelCalculation, 8, sizeof(cl_uint), &(m_inputSlices)); + clErr |= clSetKernelArg(this->m_PixelCalculation, 9, sizeof(cl_uint), &(reconstructionLines)); + clErr |= clSetKernelArg(this->m_PixelCalculation, 10, sizeof(cl_uint), &(samplesPerLine)); + } + else + { + int reconstructionLines = this->m_Conf->GetReconstructionLines(); + int samplesPerLine = this->m_Conf->GetSamplesPerLine(); + float angle = this->m_Conf->GetAngle(); + float probeRadius = this->m_Conf->GetProbeRadius(); + 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 horizontalExtent = m_Conf->GetHorizontalExtent(); + float mult = 1 / (this->m_Conf->GetTimeSpacing() * this->m_Conf->GetSpeedOfSound()); + char isPAImage = (char)m_Conf->GetIsPhotoacousticImage(); + + clErr = clSetKernelArg(this->m_PixelCalculation, 2, sizeof(cl_mem), &(this->m_ElementHeightsBuffer)); + clErr |= clSetKernelArg(this->m_PixelCalculation, 3, sizeof(cl_mem), &(this->m_ElementPositionsBuffer)); + 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_int), &(m_inputSlices)); + clErr |= clSetKernelArg(this->m_PixelCalculation, 9, sizeof(cl_int), &(reconstructionLines)); + clErr |= clSetKernelArg(this->m_PixelCalculation, 10, sizeof(cl_int), &(samplesPerLine)); + clErr |= clSetKernelArg(this->m_PixelCalculation, 11, sizeof(cl_float), &(angle)); + clErr |= clSetKernelArg(this->m_PixelCalculation, 12, sizeof(cl_float), &(probeRadius)); + clErr |= clSetKernelArg(this->m_PixelCalculation, 13, sizeof(cl_float), &(totalSamples_i)); + clErr |= clSetKernelArg(this->m_PixelCalculation, 14, sizeof(cl_float), &(horizontalExtent)); + clErr |= clSetKernelArg(this->m_PixelCalculation, 15, sizeof(cl_float), &(mult)); + clErr |= clSetKernelArg(this->m_PixelCalculation, 16, sizeof(cl_char), &(isPAImage)); + clErr |= clSetKernelArg(this->m_PixelCalculation, 17, sizeof(cl_mem), &(this->m_UsedLinesBuffer)); + } // execute the filter on a 2D/3D NDRange if (m_OutputDim[2] == 1 || m_ChunkSize[2] == 1) { if (!this->ExecuteKernelChunksInBatches(m_PixelCalculation, 2, m_ChunkSize, m_inputSlices, 50)) mitkThrow() << "openCL Error when executing Kernel"; } else { if (!this->ExecuteKernelChunksInBatches(m_PixelCalculation, 3, m_ChunkSize, m_inputSlices, 50)) mitkThrow() << "openCL Error when executing Kernel"; } // signalize the GPU-side data changed m_Output->Modified(GPU_DATA); } us::Module *mitk::PhotoacousticOCLBeamformingFilter::GetModule() { return us::GetModuleContext()->GetModule(); } bool mitk::PhotoacousticOCLBeamformingFilter::Initialize() { bool buildErr = true; cl_int clErr = 0; if (OclFilter::Initialize()) { - switch (m_Conf->GetAlgorithm()) + if (m_Conf->GetGeometry() == mitk::BeamformingSettings::ProbeGeometry::Linear) { - 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; + switch (m_Conf->GetAlgorithm()) + { + case BeamformingSettings::BeamformingAlgorithm::DAS: + { + MITK_INFO << "DAS bf"; + this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "ckDAS", &clErr); + break; + } + case BeamformingSettings::BeamformingAlgorithm::DMAS: + { + MITK_INFO << "DMAS bf"; + this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "ckDMAS", &clErr); + break; + } + case BeamformingSettings::BeamformingAlgorithm::sDMAS: + { + MITK_INFO << "sDMAS bf"; + 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; + } + } } - case BeamformingSettings::BeamformingAlgorithm::sDMAS: + else { - 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; - } + switch (m_Conf->GetAlgorithm()) + { + case BeamformingSettings::BeamformingAlgorithm::DAS: + { + MITK_INFO << "DAS bf"; + this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "ckDAS_g", &clErr); + break; + } + case BeamformingSettings::BeamformingAlgorithm::DMAS: + { + MITK_INFO << "DMAS bf"; + this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "ckDMAS_g", &clErr); + break; + } + case BeamformingSettings::BeamformingAlgorithm::sDMAS: + { + MITK_INFO << "sDMAS bf"; + this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "cksDMAS_g", &clErr); + break; + } + default: + { + MITK_INFO << "No beamforming algorithm specified, setting to DAS"; + this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "ckDAS_g", &clErr); + break; + } + } } buildErr |= CHECK_OCL_ERR(clErr); } CHECK_OCL_ERR(clErr); return (OclFilter::IsInitialized() && buildErr); } void mitk::PhotoacousticOCLBeamformingFilter::SetInput(mitk::Image::Pointer image) { OclDataSetToDataSetFilter::SetInput(image); m_InputImage = image; m_inputSlices = image->GetDimension(2); } void mitk::PhotoacousticOCLBeamformingFilter::SetInput(void* data, unsigned int* dimensions, unsigned int BpE) { OclDataSetToDataSetFilter::SetInput(data, dimensions[0] * dimensions[1] * dimensions[2], BpE); } mitk::Image::Pointer mitk::PhotoacousticOCLBeamformingFilter::GetOutputAsImage() { mitk::Image::Pointer outputImage = mitk::Image::New(); if (m_Output->IsModified(GPU_DATA)) { void* pData = m_Output->TransferDataToCPU(m_CommandQue); const unsigned int dimension = 3; unsigned int dimensions[3] = { (unsigned int)m_OutputDim[0], (unsigned int)m_OutputDim[1], (unsigned int)m_OutputDim[2] }; const mitk::SlicedGeometry3D::Pointer p_slg = m_InputImage->GetSlicedGeometry(); MITK_DEBUG << "Creating new MITK Image."; outputImage->Initialize(this->GetOutputType(), dimension, dimensions); outputImage->SetSpacing(p_slg->GetSpacing()); outputImage->SetImportVolume(pData, 0, 0, mitk::Image::ImportMemoryManagementType::CopyMemory); free(pData); } MITK_DEBUG << "Image Initialized."; return outputImage; } void* mitk::PhotoacousticOCLBeamformingFilter::GetOutput() { return OclDataSetToDataSetFilter::GetOutput(); } #endif diff --git a/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLDelayCalculation.cpp b/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLDelayCalculation.cpp index cc236090bc..e0131b7100 100644 --- a/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLDelayCalculation.cpp +++ b/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLDelayCalculation.cpp @@ -1,136 +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. ===================================================================*/ #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_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)); 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); + this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "ckDelayCalculationSphe", &clErr); buildErr |= CHECK_OCL_ERR(clErr); } return (OclFilter::IsInitialized() && buildErr); } diff --git a/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLUsedLinesCalculation.cpp b/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLUsedLinesCalculation.cpp index dda6a3eef9..8da3c4c197 100644 --- a/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLUsedLinesCalculation.cpp +++ b/Modules/PhotoacousticsAlgorithms/source/OpenCLFilter/mitkPhotoacousticOCLUsedLinesCalculation.cpp @@ -1,127 +1,167 @@ /*=================================================================== 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; } - // 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(); - char isPAImage = this->m_Conf->GetIsPhotoacousticImage(); - - 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; - - 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_char), &(isPAImage)); - clErr |= clSetKernelArg(this->m_PixelCalculation, 7, sizeof(cl_float), &(percentOfImageReconstructed)); + 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 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, 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()) { - this->m_PixelCalculation = clCreateKernel(this->m_ClProgram, "ckUsedLines", &clErr); - buildErr |= CHECK_OCL_ERR(clErr); + 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/filters/mitkBeamformingFilter.cpp b/Modules/PhotoacousticsAlgorithms/source/filters/mitkBeamformingFilter.cpp index 24a92884df..781aa4e5e2 100644 --- a/Modules/PhotoacousticsAlgorithms/source/filters/mitkBeamformingFilter.cpp +++ b/Modules/PhotoacousticsAlgorithms/source/filters/mitkBeamformingFilter.cpp @@ -1,298 +1,265 @@ /*=================================================================== mitkBeamformingFilter The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkProperties.h" #include "mitkImageReadAccessor.h" #include #include #include #include #include #include "mitkImageCast.h" #include "mitkBeamformingFilter.h" #include "mitkBeamformingUtils.h" mitk::BeamformingFilter::BeamformingFilter(mitk::BeamformingSettings::Pointer settings) : m_OutputData(nullptr), m_InputData(nullptr), m_Conf(settings) { MITK_INFO << "Instantiating BeamformingFilter..."; this->SetNumberOfIndexedInputs(1); this->SetNumberOfRequiredInputs(1); m_ProgressHandle = [](int, std::string) {}; #if defined(PHOTOACOUSTICS_USE_GPU) m_BeamformingOclFilter = mitk::PhotoacousticOCLBeamformingFilter::New(m_Conf); #else m_BeamformingOclFilter = mitk::PhotoacousticOCLBeamformingFilter::New(); #endif MITK_INFO << "Instantiating BeamformingFilter...[Done]"; } void mitk::BeamformingFilter::SetProgressHandle(std::function progressHandle) { m_ProgressHandle = progressHandle; } mitk::BeamformingFilter::~BeamformingFilter() { MITK_INFO << "Destructed BeamformingFilter"; } void mitk::BeamformingFilter::GenerateInputRequestedRegion() { Superclass::GenerateInputRequestedRegion(); mitk::Image* output = this->GetOutput(); mitk::Image* input = const_cast (this->GetInput()); if (!output->IsInitialized()) { return; } input->SetRequestedRegionToLargestPossibleRegion(); } void mitk::BeamformingFilter::GenerateOutputInformation() { mitk::Image::ConstPointer input = this->GetInput(); mitk::Image::Pointer output = this->GetOutput(); if ((output->IsInitialized()) && (this->GetMTime() <= m_TimeOfHeaderInitialization.GetMTime())) return; mitk::Vector3D spacing; - spacing[0] = m_Conf->GetPitchInMeters() * m_Conf->GetTransducerElements() * 1000 / m_Conf->GetReconstructionLines(); + spacing[0] = m_Conf->GetHorizontalExtent() / m_Conf->GetReconstructionLines() * 1000; float desiredYSpacing = m_Conf->GetReconstructionDepth() * 1000 / m_Conf->GetSamplesPerLine(); float maxYSpacing = m_Conf->GetSpeedOfSound() * m_Conf->GetTimeSpacing() * input->GetDimension(1) / m_Conf->GetSamplesPerLine() * 1000; spacing[1] = desiredYSpacing < maxYSpacing ? desiredYSpacing : maxYSpacing; spacing[2] = 1; unsigned int dim[] = { m_Conf->GetReconstructionLines(), m_Conf->GetSamplesPerLine(), input->GetDimension(2)}; output->Initialize(mitk::MakeScalarPixelType(), 3, dim); output->GetGeometry()->SetSpacing(spacing); output->GetGeometry()->Modified(); output->SetPropertyList(input->GetPropertyList()->Clone()); m_TimeOfHeaderInitialization.Modified(); } void mitk::BeamformingFilter::GenerateData() { mitk::Image::Pointer input = this->GetInput(); if (!(input->GetPixelType().GetTypeAsString() == "scalar (float)" || input->GetPixelType().GetTypeAsString() == " (float)")) { MITK_ERROR << "Pixel type of input needs to be float for this filter to work."; mitkThrow() << "Pixel type of input needs to be float for this filter to work."; } GenerateOutputInformation(); mitk::Image::Pointer output = this->GetOutput(); if (!output->IsInitialized()) return; auto begin = std::chrono::high_resolution_clock::now(); // debbuging the performance... if (!m_Conf->GetUseGPU()) { int progInterval = output->GetDimension(2) / 20 > 1 ? output->GetDimension(2) / 20 : 1; // the interval at which we update the gui progress bar float inputDim[2] = { (float)input->GetDimension(0), (float)input->GetDimension(1) }; float outputDim[2] = { (float)output->GetDimension(0), (float)output->GetDimension(1) }; for (unsigned int i = 0; i < output->GetDimension(2); ++i) // seperate Slices should get Beamforming seperately applied { mitk::ImageReadAccessor inputReadAccessor(input, input->GetSliceData(i)); m_InputData = (float*)inputReadAccessor.GetData(); m_OutputData = new float[m_Conf->GetReconstructionLines()*m_Conf->GetSamplesPerLine()]; // fill the image with zeros for (int l = 0; l < outputDim[0]; ++l) { for (int s = 0; s < outputDim[1]; ++s) { m_OutputData[l*(short)outputDim[1] + s] = 0; } } std::thread *threads = new std::thread[(short)outputDim[0]]; // every line will be beamformed in a seperate thread if (m_Conf->GetAlgorithm() == BeamformingSettings::BeamformingAlgorithm::DAS) { - if (m_Conf->GetDelayCalculationMethod() == BeamformingSettings::DelayCalc::QuadApprox) + for (short line = 0; line < outputDim[0]; ++line) { - 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); - } + 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) { - for (short line = 0; line < outputDim[0]; ++line) - { - threads[line] = std::thread(&BeamformingUtils::DMASSphericalLine, m_InputData, - m_OutputData, inputDim, outputDim, line, m_Conf); - } + 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) { - for (short line = 0; line < outputDim[0]; ++line) - { - threads[line] = std::thread(&BeamformingUtils::sDMASSphericalLine, m_InputData, - m_OutputData, inputDim, outputDim, line, m_Conf); - } + threads[line] = std::thread(&BeamformingUtils::sDMASSphericalLine, m_InputData, + m_OutputData, inputDim, outputDim, line, m_Conf); } } // wait for all lines to finish for (short line = 0; line < outputDim[0]; ++line) { threads[line].join(); } output->SetSlice(m_OutputData, i); if (i % progInterval == 0) m_ProgressHandle((int)((i + 1) / (float)output->GetDimension(2) * 100), "performing reconstruction"); delete[] m_OutputData; m_OutputData = nullptr; m_InputData = nullptr; } } #if defined(PHOTOACOUSTICS_USE_GPU) || DOXYGEN else { try { // first, we check whether the data is float, other formats are unsupported if (!(input->GetPixelType().GetTypeAsString() == "scalar (float)" || input->GetPixelType().GetTypeAsString() == " (float)")) { MITK_ERROR << "Pixel type is not float, abort"; mitkThrow() << "Pixel type is not float, abort"; } unsigned long availableMemory = m_BeamformingOclFilter->GetDeviceMemory(); unsigned int batchSize = m_Conf->GetGPUBatchSize(); unsigned int batches = (unsigned int)((float)input->GetDimension(2) / batchSize) + (input->GetDimension(2) % batchSize > 0); unsigned int batchDim[] = { input->GetDimension(0), input->GetDimension(1), batchSize }; unsigned int batchDimLast[] = { input->GetDimension(0), input->GetDimension(1), input->GetDimension(2) % batchSize }; // the following safeguard is probably only needed for absurdly small GPU memory if((unsigned long)batchSize * ((unsigned long)(batchDim[0] * batchDim[1]) * 4 + // single input image (float) (unsigned long)(m_Conf->GetReconstructionLines() * m_Conf->GetSamplesPerLine()) * 4) // single output image (float) > availableMemory - (unsigned long)(m_Conf->GetReconstructionLines() / 2 * m_Conf->GetSamplesPerLine()) * 2 - // Delays buffer (unsigned short) (unsigned long)(m_Conf->GetReconstructionLines() * m_Conf->GetSamplesPerLine()) * 3 * 2 - // UsedLines buffer (unsigned short) 50 * 1024 * 1024)// 50 MB buffer for local data, system purposes etc { MITK_ERROR << "device memory too small for GPU beamforming; try decreasing the batch size"; return; } mitk::ImageReadAccessor copy(input); for (unsigned int i = 0; i < batches; ++i) { m_ProgressHandle(100.f * (float)i / (float)batches, "performing reconstruction"); mitk::Image::Pointer inputBatch = mitk::Image::New(); unsigned int num_Slices = 1; if (i == batches - 1 && (input->GetDimension(2) % batchSize > 0)) { inputBatch->Initialize(mitk::MakeScalarPixelType(), 3, batchDimLast); num_Slices = batchDimLast[2]; } else { inputBatch->Initialize(mitk::MakeScalarPixelType(), 3, batchDim); num_Slices = batchDim[2]; } inputBatch->SetSpacing(input->GetGeometry()->GetSpacing()); inputBatch->SetImportVolume(&(((float*)copy.GetData())[input->GetDimension(0) * input->GetDimension(1) * batchSize * i])); m_BeamformingOclFilter->SetApodisation(m_Conf->GetApodizationFunction(), m_Conf->GetApodizationArraySize()); m_BeamformingOclFilter->SetInput(inputBatch); m_BeamformingOclFilter->Update(); void* out = m_BeamformingOclFilter->GetOutput(); for (unsigned int slice = 0; slice < num_Slices; ++slice) { output->SetImportSlice( &(((float*)out)[m_Conf->GetReconstructionLines() * m_Conf->GetSamplesPerLine() * slice]), batchSize * i + slice, 0, 0, mitk::Image::ImportMemoryManagementType::CopyMemory); } free(out); } } catch (mitk::Exception &e) { std::string errorMessage = "Caught unexpected exception "; errorMessage.append(e.what()); MITK_ERROR << errorMessage; float* dummyData = new float[m_Conf->GetReconstructionLines() * m_Conf->GetSamplesPerLine() * m_Conf->GetInputDim()[2]]; output->SetImportVolume(dummyData, 0, 0, mitk::Image::ImportMemoryManagementType::CopyMemory); } } #endif m_TimeOfHeaderInitialization.Modified(); auto end = std::chrono::high_resolution_clock::now(); MITK_INFO << "Beamforming of " << output->GetDimension(2) << " Images completed in " << ((float)std::chrono::duration_cast(end - begin).count()) / 1000000 << "ms" << std::endl; } diff --git a/Modules/PhotoacousticsAlgorithms/source/filters/mitkBeamformingSettings.cpp b/Modules/PhotoacousticsAlgorithms/source/filters/mitkBeamformingSettings.cpp index 7508605611..894802fe76 100644 --- a/Modules/PhotoacousticsAlgorithms/source/filters/mitkBeamformingSettings.cpp +++ b/Modules/PhotoacousticsAlgorithms/source/filters/mitkBeamformingSettings.cpp @@ -1,96 +1,141 @@ /*=================================================================== mitkBeamformingSettings 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 "mitkBeamformingSettings.h" #include "mitkBeamformingUtils.h" #include "itkMutexLock.h" mitk::BeamformingSettings::BeamformingSettings(float pitchInMeters, float speedOfSound, float timeSpacing, float angle, bool isPhotoacousticImage, unsigned int samplesPerLine, unsigned int reconstructionLines, unsigned int* inputDim, float reconstructionDepth, bool useGPU, unsigned int GPUBatchSize, - DelayCalc delayCalculationMethod, Apodization apod, unsigned int apodizationArraySize, - BeamformingAlgorithm algorithm + BeamformingAlgorithm algorithm, + ProbeGeometry geometry, + float probeRadius ) : m_PitchInMeters(pitchInMeters), m_SpeedOfSound(speedOfSound), m_TimeSpacing(timeSpacing), m_Angle(angle), m_IsPhotoacousticImage(isPhotoacousticImage), m_SamplesPerLine(samplesPerLine), m_ReconstructionLines(reconstructionLines), m_ReconstructionDepth(reconstructionDepth), m_UseGPU(useGPU), m_GPUBatchSize(GPUBatchSize), - m_DelayCalculationMethod(delayCalculationMethod), m_Apod(apod), m_ApodizationArraySize(apodizationArraySize), - m_Algorithm(algorithm) + m_Algorithm(algorithm), + m_Geometry(geometry), + m_ProbeRadius(probeRadius), + m_MinMaxLines(nullptr) { if (inputDim == nullptr) { MITK_ERROR << "No input dimension given."; mitkThrow() << "No input dimension given."; } switch (GetApod()) { case BeamformingSettings::Apodization::Hann: m_ApodizationFunction = mitk::BeamformingUtils::VonHannFunction(GetApodizationArraySize()); break; case BeamformingSettings::Apodization::Hamm: m_ApodizationFunction = mitk::BeamformingUtils::HammFunction(GetApodizationArraySize()); break; case BeamformingSettings::Apodization::Box: default: m_ApodizationFunction = mitk::BeamformingUtils::BoxFunction(GetApodizationArraySize()); break; } - m_InputDim = new unsigned int[3]{ inputDim[0], inputDim[1], inputDim[2] }; - m_TransducerElements = m_InputDim[0]; + + m_ElementHeights = new float[m_TransducerElements]; + m_ElementPositions = new float[m_TransducerElements]; + + if (m_Geometry == ProbeGeometry::Concave) + { + float openingAngle = (m_TransducerElements * m_PitchInMeters) / (probeRadius * 2 * itk::Math::pi) * 2 * itk::Math::pi; + m_HorizontalExtent = std::sin(openingAngle / 2.f) * probeRadius * 2.f; + + float elementAngle = 0; + + for (unsigned int i = 0; i < m_TransducerElements; ++i) + { + elementAngle = ((i- m_TransducerElements /2.f) * m_PitchInMeters) / (probeRadius * 2 * itk::Math::pi) * 2 * itk::Math::pi; + m_ElementHeights[i] = probeRadius - std::cos(elementAngle) * probeRadius; + m_ElementPositions[i] = m_HorizontalExtent/2.f + std::sin(elementAngle) * probeRadius; + } + } + else + { + m_HorizontalExtent = m_PitchInMeters * m_TransducerElements; + for (unsigned int i = 0; i < m_TransducerElements; ++i) + { + m_ElementHeights[i] = 0; + m_ElementPositions[i] = i * m_PitchInMeters; + } + } } mitk::BeamformingSettings::~BeamformingSettings() { MITK_INFO << "Destructing beamforming settings..."; //Free memory if (m_ApodizationFunction != nullptr) { MITK_INFO << "Deleting apodization function..."; delete[] m_ApodizationFunction; MITK_INFO << "Deleting apodization function...[Done]"; } if (m_InputDim != nullptr) { MITK_INFO << "Deleting input dim..."; delete[] m_InputDim; MITK_INFO << "Deleting input dim...[Done]"; } + if (m_ElementHeights != nullptr || m_ElementPositions != nullptr) + { + MITK_INFO << "Deleting element geometry..."; + if (m_ElementHeights != nullptr) + delete[] m_ElementHeights; - MITK_INFO << "Destructing beamforming settings...[Done]"; + if (m_ElementPositions != nullptr) + delete[] m_ElementPositions; + MITK_INFO << "Destructing beamforming settings...[Done]"; + } + if (m_MinMaxLines) + delete[] m_MinMaxLines; } + +unsigned short* mitk::BeamformingSettings::GetMinMaxLines() +{ + if (!m_MinMaxLines) + m_MinMaxLines = mitk::BeamformingUtils::MinMaxLines(this); + return m_MinMaxLines; +} \ No newline at end of file diff --git a/Modules/PhotoacousticsAlgorithms/source/utils/mitkBeamformingUtils.cpp b/Modules/PhotoacousticsAlgorithms/source/utils/mitkBeamformingUtils.cpp index b000906692..3162ce7f32 100644 --- a/Modules/PhotoacousticsAlgorithms/source/utils/mitkBeamformingUtils.cpp +++ b/Modules/PhotoacousticsAlgorithms/source/utils/mitkBeamformingUtils.cpp @@ -1,563 +1,435 @@ /*=================================================================== 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) +unsigned short* mitk::BeamformingUtils::MinMaxLines(const mitk::BeamformingSettings::Pointer config) { - const float* apodisation = config->GetApodizationFunction(); - const short apodArraySize = config->GetApodizationArraySize(); - float& inputS = inputDim[1]; - float& inputL = inputDim[0]; + int outputL = (int)config->GetReconstructionLines(); + int outputS = (int)config->GetSamplesPerLine(); - float& outputS = outputDim[1]; - float& outputL = outputDim[0]; + unsigned short* dDest = new unsigned short[outputL * outputS * 2]; - short AddSample = 0; - short maxLine = 0; - short minLine = 0; - float delayMultiplicator = 0; - float l_i = 0; - float s_i = 0; + int inputL = (int)config->GetInputDim()[0]; + int inputS = (int)config->GetInputDim()[1]; + + float horizontalExtent = config->GetHorizontalExtent(); + float verticalExtent = config->GetReconstructionDepth(); - 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(); + float* elementHeights = config->GetElementHeights(); + float* elementPositions = config->GetElementPositions(); - short usedLines = (maxLine - minLine); - - float percentOfImageReconstructed = (float)(config->GetReconstructionDepth()) / - (float)(inputS * config->GetSpeedOfSound() * config->GetTimeSpacing() / (float)(2 - (int)config->GetIsPhotoacousticImage())); - percentOfImageReconstructed = percentOfImageReconstructed <= 1 ? percentOfImageReconstructed : 1; + float cos_deg = std::cos(config->GetAngle() / 2.f / 360 * 2 * itk::Math::pi); - l_i = (float)line / outputL * inputL; + float cos = 0; + float a = 0; + float d = 0; - for (short sample = 0; sample < outputS; ++sample) + for (int x = 0; x < outputL; ++x) { - s_i = (float)sample / outputS * inputS / (float)(2 - (int)config->GetIsPhotoacousticImage()) * percentOfImageReconstructed; + for (int y = 0; y < outputS; ++y) + { + float l_p = (float)x / outputL * horizontalExtent; + float s_p = (float)y / (float)outputS * verticalExtent; - part = part_multiplicator*s_i; + int maxLine = inputL; + int minLine = 0; - if (part < 1) - part = 1; + 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); - maxLine = (short)std::min((l_i + part) + 1, inputL); - minLine = (short)std::max((l_i - part), 0.0f); - usedLines = (maxLine - minLine); + 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); - apod_mult = (float)apodArraySize / (float)usedLines; + 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) + { + 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); - delayMultiplicator = pow((1 / (config->GetTimeSpacing()*config->GetSpeedOfSound()) * - (config->GetPitchInMeters()*config->GetTransducerElements()) / inputL), 2) / s_i / 2; + if (cos > cos_deg) + { + minLine = l_s; + break; + } + } - 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; + 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) + { + 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 } - output[sample*(short)outputL + line] = output[sample*(short)outputL + line] / usedLines; } + 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_i = 0; + float l_p = 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(); + short usedLines = (maxLine - minLine); - float percentOfImageReconstructed = (float)(config->GetReconstructionDepth()) / - (float)(inputS * config->GetSpeedOfSound() * config->GetTimeSpacing() / (float)(2 - (int)config->GetIsPhotoacousticImage())); - percentOfImageReconstructed = percentOfImageReconstructed <= 1 ? percentOfImageReconstructed : 1; + 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; + l_p = (float)line / outputL * config->GetHorizontalExtent(); for (short sample = 0; sample < outputS; ++sample) { - s_i = (float)sample / outputS * inputS / (float)(2 - (int)config->GetIsPhotoacousticImage()) * percentOfImageReconstructed; - - part = part_multiplicator*s_i; + s_i = (float)sample / outputS * totalSamples_i; - if (part < 1) - part = 1; - - maxLine = (short)std::min((l_i + part) + 1, inputL); - minLine = (short)std::max((l_i - part), 0.0f); + 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, 2) + pow(s_i-elementHeights[l_s]/(config->GetSpeedOfSound()*config->GetTimeSpacing()), 2) + - pow((1 / (config->GetTimeSpacing()*config->GetSpeedOfSound()) * - (((float)l_s - l_i)*config->GetPitchInMeters()*(float)config->GetTransducerElements()) / inputL), 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::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); - - float percentOfImageReconstructed = (float)(config->GetReconstructionDepth()) / - (float)(inputS * config->GetSpeedOfSound() * config->GetTimeSpacing() / (float)(2 - (int)config->GetIsPhotoacousticImage())); - percentOfImageReconstructed = percentOfImageReconstructed <= 1 ? percentOfImageReconstructed : 1; - - l_i = (float)line / outputL * inputL; - - for (short sample = 0; sample < outputS; ++sample) - { - s_i = (float)sample / outputS * inputS / (float)(2 - (int)config->GetIsPhotoacousticImage()) * percentOfImageReconstructed; - - 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(); + 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_i = 0; + float l_p = 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(); + float mult = 0; short usedLines = (maxLine - minLine); - float percentOfImageReconstructed = (float)(config->GetReconstructionDepth()) / - (float)(inputS * config->GetSpeedOfSound() * config->GetTimeSpacing() / (float)(2 - (int)config->GetIsPhotoacousticImage())); - percentOfImageReconstructed = percentOfImageReconstructed <= 1 ? percentOfImageReconstructed : 1; + 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; + l_p = (float)line / outputL * config->GetHorizontalExtent(); for (short sample = 0; sample < outputS; ++sample) { - s_i = (float)sample / outputS * inputS / (float)(2 - (int)config->GetIsPhotoacousticImage()) * percentOfImageReconstructed; + s_i = (float)sample / outputS * totalSamples_i; - 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); + 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] = (short)sqrt( - pow(s_i, 2) + AddSample[l_s] = (int)sqrt( + pow(s_i - elementHeights[l_s + minLine] / (config->GetSpeedOfSound()*config->GetTimeSpacing()), 2) + - pow((1 / (config->GetTimeSpacing()*config->GetSpeedOfSound()) * - (((float)minLine + (float)l_s - l_i)*config->GetPitchInMeters()*(float)config->GetTransducerElements()) / inputL), 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::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); - - float percentOfImageReconstructed = (float)(config->GetReconstructionDepth()) / - (float)(inputS * config->GetSpeedOfSound() * config->GetTimeSpacing() / (float)(2 - (int)config->GetIsPhotoacousticImage())); - percentOfImageReconstructed = percentOfImageReconstructed <= 1 ? percentOfImageReconstructed : 1; - - l_i = (float)line / outputL * inputL; - - for (short sample = 0; sample < outputS; ++sample) - { - s_i = (float)sample / outputS * inputS / (float)(2 - (int)config->GetIsPhotoacousticImage()) * percentOfImageReconstructed; - - 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(); + 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_i = 0; + float l_p = 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 percentOfImageReconstructed = (float)(config->GetReconstructionDepth()) / - (float)(inputS * config->GetSpeedOfSound() * config->GetTimeSpacing() / (float)(2 - (int)config->GetIsPhotoacousticImage())); - percentOfImageReconstructed = percentOfImageReconstructed <= 1 ? percentOfImageReconstructed : 1; + 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; + l_p = (float)line / outputL * config->GetHorizontalExtent(); for (short sample = 0; sample < outputS; ++sample) { - s_i = (float)sample / outputS * inputS / (float)(2 - (int)config->GetIsPhotoacousticImage()) * percentOfImageReconstructed; - - part = part_multiplicator*s_i; - - if (part < 1) - part = 1; + s_i = (float)sample / outputS * totalSamples_i; - maxLine = (short)std::min((l_i + part) + 1, inputL); - minLine = (short)std::max((l_i - part), 0.0f); + 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] = (short)sqrt( - pow(s_i, 2) + AddSample[l_s] = (int)sqrt( + pow(s_i - elementHeights[l_s + minLine] / (config->GetSpeedOfSound()*config->GetTimeSpacing()), 2) + - pow((1 / (config->GetTimeSpacing()*config->GetSpeedOfSound()) * - (((float)minLine + (float)l_s - l_i)*config->GetPitchInMeters()*(float)config->GetTransducerElements()) / inputL), 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; } } diff --git a/Modules/PhotoacousticsAlgorithms/source/utils/mitkPhotoacousticFilterService.cpp b/Modules/PhotoacousticsAlgorithms/source/utils/mitkPhotoacousticFilterService.cpp index eb3940e9cc..ddcfed7aa4 100644 --- a/Modules/PhotoacousticsAlgorithms/source/utils/mitkPhotoacousticFilterService.cpp +++ b/Modules/PhotoacousticsAlgorithms/source/utils/mitkPhotoacousticFilterService.cpp @@ -1,307 +1,306 @@ /*=================================================================== 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 "mitkPhotoacousticFilterService.h" #include "mitkITKImageImport.h" #include #include #include "./OpenCLFilter/mitkPhotoacousticBModeFilter.h" #include "mitkConvert2Dto3DImageFilter.h" #include #include "../ITKFilter/ITKUltrasound/itkBModeImageFilter.h" #include "../ITKFilter/itkPhotoacousticBModeImageFilter.h" #include // itk dependencies #include "itkImage.h" #include "itkResampleImageFilter.h" #include "itkCastImageFilter.h" #include "itkCropImageFilter.h" #include "itkRescaleIntensityImageFilter.h" #include "itkIntensityWindowingImageFilter.h" #include #include "itkBSplineInterpolateImageFunction.h" #include // needed itk image filters #include "mitkImageCast.h" mitk::PhotoacousticFilterService::PhotoacousticFilterService() { MITK_INFO << "[PhotoacousticFilterService] created filter service"; } mitk::PhotoacousticFilterService::~PhotoacousticFilterService() { MITK_INFO << "[PhotoacousticFilterService] destructed filter service"; } mitk::Image::Pointer mitk::PhotoacousticFilterService::ApplyBmodeFilter( mitk::Image::Pointer inputImage, BModeMethod method, bool UseLogFilter) { // the image needs to be of floating point type for the envelope filter to work; the casting is done automatically by the CastToItkImage typedef itk::Image< float, 3 > itkFloatImageType; auto floatImage = ConvertToFloat(inputImage); if (method == BModeMethod::Abs) { PhotoacousticBModeFilter::Pointer filter = PhotoacousticBModeFilter::New(); filter->UseLogFilter(UseLogFilter); filter->SetInput(floatImage); filter->Update(); return filter->GetOutput(); } typedef itk::BModeImageFilter < itkFloatImageType, itkFloatImageType > BModeFilterType; BModeFilterType::Pointer bModeFilter = BModeFilterType::New(); // LogFilter typedef itk::PhotoacousticBModeImageFilter < itkFloatImageType, itkFloatImageType > PhotoacousticBModeImageFilter; PhotoacousticBModeImageFilter::Pointer photoacousticBModeFilter = PhotoacousticBModeImageFilter::New(); // No LogFilter typedef itk::ResampleImageFilter < itkFloatImageType, itkFloatImageType > ResampleImageFilter; ResampleImageFilter::Pointer resampleImageFilter = ResampleImageFilter::New(); itkFloatImageType::Pointer itkImage; mitk::CastToItkImage(floatImage, itkImage); if (UseLogFilter) { bModeFilter->SetInput(itkImage); bModeFilter->SetDirection(1); itkImage = bModeFilter->GetOutput(); } else { photoacousticBModeFilter->SetInput(itkImage); photoacousticBModeFilter->SetDirection(1); itkImage = photoacousticBModeFilter->GetOutput(); } return mitk::GrabItkImageMemory(itkImage); } mitk::Image::Pointer mitk::PhotoacousticFilterService::ApplyResampling( mitk::Image::Pointer inputImage, double *outputSpacing) { typedef itk::Image< float, 3 > itkFloatImageType; auto floatImage = ConvertToFloat(inputImage); typedef itk::ResampleImageFilter < itkFloatImageType, itkFloatImageType > ResampleImageFilter; ResampleImageFilter::Pointer resampleImageFilter = ResampleImageFilter::New(); itkFloatImageType::Pointer itkImage; mitk::CastToItkImage(floatImage, itkImage); itkFloatImageType::SpacingType outputSpacingItk; itkFloatImageType::SizeType inputSizeItk = itkImage->GetLargestPossibleRegion().GetSize(); itkFloatImageType::SizeType outputSizeItk = inputSizeItk; outputSpacingItk[0] = outputSpacing[0]; outputSpacingItk[1] = outputSpacing[1]; outputSpacingItk[2] = itkImage->GetSpacing()[2]; outputSizeItk[0] = outputSizeItk[0] * (floatImage->GetGeometry()->GetSpacing()[0] / outputSpacing[0]); outputSizeItk[1] = outputSizeItk[1] * (floatImage->GetGeometry()->GetSpacing()[1] / outputSpacing[1]); resampleImageFilter->SetInput(itkImage); resampleImageFilter->SetSize(outputSizeItk); resampleImageFilter->SetOutputSpacing(outputSpacingItk); resampleImageFilter->UpdateLargestPossibleRegion(); return mitk::GrabItkImageMemory(resampleImageFilter->GetOutput()); } mitk::Image::Pointer mitk::PhotoacousticFilterService::ApplyResamplingToDim( mitk::Image::Pointer inputImage, double *outputDimension) { typedef itk::Image< float, 3 > itkFloatImageType; auto floatImage = ConvertToFloat(inputImage); typedef itk::ResampleImageFilter < itkFloatImageType, itkFloatImageType > ResampleImageFilter; ResampleImageFilter::Pointer resampleImageFilter = ResampleImageFilter::New(); itkFloatImageType::Pointer itkImage; mitk::CastToItkImage(floatImage, itkImage); itkFloatImageType::SpacingType outputSpacingItk; itkFloatImageType::SizeType inputSizeItk = itkImage->GetLargestPossibleRegion().GetSize(); itkFloatImageType::SizeType outputSizeItk = inputSizeItk; outputSizeItk[0] = outputDimension[0]; outputSizeItk[1] = outputDimension[1]; MITK_INFO << outputSizeItk[0] << " " << outputSizeItk[1]; outputSpacingItk[0] = (double)inputSizeItk[0] / (double)outputSizeItk[0] * floatImage->GetGeometry()->GetSpacing()[0]; outputSpacingItk[1] = (double)inputSizeItk[1] / (double)outputSizeItk[1] * floatImage->GetGeometry()->GetSpacing()[1]; outputSpacingItk[2] = itkImage->GetSpacing()[2]; MITK_INFO << outputSpacingItk[0] << " " << outputSpacingItk[1]; resampleImageFilter->SetInput(itkImage); resampleImageFilter->SetSize(outputSizeItk); resampleImageFilter->SetOutputSpacing(outputSpacingItk); resampleImageFilter->UpdateLargestPossibleRegion(); return mitk::GrabItkImageMemory(resampleImageFilter->GetOutput()); } mitk::Image::Pointer mitk::PhotoacousticFilterService::ApplyCropping( mitk::Image::Pointer inputImage, int above, int below, int right, int left, int zStart, int zEnd, int* errCode) { *errCode = 0; try { auto floatImage = ConvertToFloat(inputImage); mitk::CropImageFilter::Pointer cropImageFilter = mitk::CropImageFilter::New(); cropImageFilter->SetInput(floatImage); cropImageFilter->SetXPixelsCropStart(left); cropImageFilter->SetXPixelsCropEnd(right); cropImageFilter->SetYPixelsCropStart(above); cropImageFilter->SetYPixelsCropEnd(below); cropImageFilter->SetZPixelsCropStart(zStart); cropImageFilter->SetZPixelsCropEnd(zEnd); cropImageFilter->Update(); return cropImageFilter->GetOutput(); } catch (mitk::Exception &e) { std::string errorMessage = "Caught unexpected exception "; errorMessage.append(e.what()); MITK_ERROR << errorMessage; *errCode = -1; mitk::Image::Pointer ret = mitk::Image::New(); unsigned int dim[3] = { 1,1,1 }; ret->Initialize(MakeScalarPixelType(), 3, dim); return ret; } } mitk::Image::Pointer mitk::PhotoacousticFilterService::ExtendImage(mitk::Image::Pointer inputImage, float pixelColor, unsigned int outputDimensionY) { mitk::Image::Pointer outputImage = mitk::Image::New(); unsigned int dim[] = {inputImage->GetDimension(0), outputDimensionY, inputImage->GetDimension(2)}; outputImage->Initialize(inputImage->GetPixelType(), 3, dim); float *sliceData = new float[dim[0] * dim[1]]; for (size_t i = inputImage->GetDimension(1) * dim[0]; i < dim[0] * dim[1]; ++i) { sliceData[i] = pixelColor; } for (unsigned int slice = 0; slice < dim[2]; ++slice) { mitk::ImageReadAccessor cpy(inputImage, inputImage->GetSliceData(slice)); cpy.GetData(); std::memcpy((void*)sliceData, cpy.GetData(), sizeof(float) * inputImage->GetDimension(1) * dim[0]); outputImage->SetSlice(sliceData, slice); } delete[] sliceData; return outputImage; } mitk::Image::Pointer mitk::PhotoacousticFilterService::ApplyBeamforming( mitk::Image::Pointer inputImage, BeamformingSettings::Pointer config, std::function progressHandle) { Image::Pointer processedImage = mitk::Image::New(); if (inputImage->GetDimension() != 3) { mitk::Convert2Dto3DImageFilter::Pointer dimensionImageFilter = mitk::Convert2Dto3DImageFilter::New(); dimensionImageFilter->SetInput(inputImage); dimensionImageFilter->Update(); processedImage = dimensionImageFilter->GetOutput(); } else { processedImage = inputImage; } m_BeamformingFilter = mitk::BeamformingFilter::New(config); m_BeamformingFilter->SetInput(ConvertToFloat(processedImage)); m_BeamformingFilter->SetProgressHandle(progressHandle); m_BeamformingFilter->UpdateLargestPossibleRegion(); processedImage = m_BeamformingFilter->GetOutput(); return processedImage; } mitk::Image::Pointer mitk::PhotoacousticFilterService::ApplyBandpassFilter( mitk::Image::Pointer data, float BPHighPass, float BPLowPass, float alphaHighPass, float alphaLowPass, float TimeSpacing, float SpeedOfSound, bool IsBFImage) { try { auto floatData = ConvertToFloat(data); mitk::BandpassFilter::Pointer bandpassFilter = mitk::BandpassFilter::New(); bandpassFilter->SetInput(floatData); bandpassFilter->SetHighPass(BPHighPass); bandpassFilter->SetLowPass(BPLowPass); bandpassFilter->SetHighPassAlpha(alphaHighPass); bandpassFilter->SetLowPassAlpha(alphaLowPass); bandpassFilter->SetSpeedOfSound(SpeedOfSound); bandpassFilter->SetTimeSpacing(TimeSpacing); bandpassFilter->SetIsBFImage(IsBFImage); bandpassFilter->Update(); return bandpassFilter->GetOutput(); } - catch (mitk::Exception &e) { std::string errorMessage = "Caught unexpected exception "; errorMessage.append(e.what()); MITK_ERROR << errorMessage; return data; } } mitk::Image::Pointer mitk::PhotoacousticFilterService::ConvertToFloat(mitk::Image::Pointer inputImage) { if ((inputImage->GetPixelType().GetTypeAsString() == "scalar (float)" || inputImage->GetPixelType().GetTypeAsString() == " (float)")) { return inputImage; } mitk::CastToFloatImageFilter::Pointer castToFloatImageFilter = mitk::CastToFloatImageFilter::New(); castToFloatImageFilter->SetInput(inputImage); castToFloatImageFilter->Update(); return castToFloatImageFilter->GetOutput(); } diff --git a/Plugins/org.mitk.gui.qt.photoacoustics.imageprocessing/src/internal/PAImageProcessing.cpp b/Plugins/org.mitk.gui.qt.photoacoustics.imageprocessing/src/internal/PAImageProcessing.cpp index 1a48655e57..f4834082f3 100644 --- a/Plugins/org.mitk.gui.qt.photoacoustics.imageprocessing/src/internal/PAImageProcessing.cpp +++ b/Plugins/org.mitk.gui.qt.photoacoustics.imageprocessing/src/internal/PAImageProcessing.cpp @@ -1,1122 +1,1147 @@ /*=================================================================== 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. ===================================================================*/ // Blueberry #include #include // Qmitk #include "PAImageProcessing.h" // Qt #include #include #include #include //mitk image #include #include "mitkPhotoacousticFilterService.h" #include "mitkCastToFloatImageFilter.h" #include "mitkBeamformingFilter.h" //other #include #include #include #define GPU_BATCH_SIZE 32 const std::string PAImageProcessing::VIEW_ID = "org.mitk.views.paimageprocessing"; PAImageProcessing::PAImageProcessing() : m_ResampleSpacing(0), m_UseLogfilter(false), m_FilterBank(mitk::PhotoacousticFilterService::New()) { qRegisterMetaType(); qRegisterMetaType(); } void PAImageProcessing::SetFocus() { m_Controls.buttonApplyBModeFilter->setFocus(); } void PAImageProcessing::CreateQtPartControl(QWidget *parent) { // create GUI widgets from the Qt Designer's .ui file m_Controls.setupUi(parent); connect(m_Controls.buttonApplyBModeFilter, SIGNAL(clicked()), this, SLOT(StartBmodeThread())); + connect(m_Controls.Geometry, SIGNAL(currentIndexChanged(int)), this, SLOT(ChangedProbe())); connect(m_Controls.DoResampling, SIGNAL(clicked()), this, SLOT(UseResampling())); connect(m_Controls.Logfilter, SIGNAL(clicked()), this, SLOT(UseLogfilter())); connect(m_Controls.ResamplingValue, SIGNAL(valueChanged(double)), this, SLOT(SetResampling())); connect(m_Controls.buttonApplyBeamforming, SIGNAL(clicked()), this, SLOT(StartBeamformingThread())); connect(m_Controls.buttonApplyCropFilter, SIGNAL(clicked()), this, SLOT(StartCropThread())); connect(m_Controls.buttonApplyBandpass, SIGNAL(clicked()), this, SLOT(StartBandpassThread())); connect(m_Controls.UseImageSpacing, SIGNAL(clicked()), this, SLOT(UseImageSpacing())); connect(m_Controls.ScanDepth, SIGNAL(valueChanged(double)), this, SLOT(UpdateImageInfo())); connect(m_Controls.SpeedOfSound, SIGNAL(valueChanged(double)), this, SLOT(UpdateImageInfo())); connect(m_Controls.SpeedOfSound, SIGNAL(valueChanged(double)), this, SLOT(ChangedSOSBeamforming())); connect(m_Controls.BPSpeedOfSound, SIGNAL(valueChanged(double)), this, SLOT(ChangedSOSBandpass())); connect(m_Controls.Samples, SIGNAL(valueChanged(int)), this, SLOT(UpdateImageInfo())); connect(m_Controls.UseImageSpacing, SIGNAL(clicked()), this, SLOT(UpdateImageInfo())); connect(m_Controls.boundLow, SIGNAL(valueChanged(int)), this, SLOT(LowerSliceBoundChanged())); connect(m_Controls.boundHigh, SIGNAL(valueChanged(int)), this, SLOT(UpperSliceBoundChanged())); connect(m_Controls.Partial, SIGNAL(clicked()), this, SLOT(SliceBoundsEnabled())); connect(m_Controls.BatchProcessing, SIGNAL(clicked()), this, SLOT(BatchProcessing())); connect(m_Controls.StepBeamforming, SIGNAL(clicked()), this, SLOT(UpdateSaveBoxes())); connect(m_Controls.StepCropping, SIGNAL(clicked()), this, SLOT(UpdateSaveBoxes())); connect(m_Controls.StepBandpass, SIGNAL(clicked()), this, SLOT(UpdateSaveBoxes())); connect(m_Controls.StepBMode, SIGNAL(clicked()), this, SLOT(UpdateSaveBoxes())); connect(m_Controls.UseSignalDelay, SIGNAL(clicked()), this, SLOT(UseSignalDelay())); connect(m_Controls.IsBFImage, SIGNAL(clicked()), this, SLOT(UpdateImageInfo())); UpdateSaveBoxes(); UseSignalDelay(); m_Controls.DoResampling->setChecked(false); m_Controls.ResamplingValue->setEnabled(false); m_Controls.progressBar->setMinimum(0); m_Controls.progressBar->setMaximum(100); m_Controls.progressBar->setVisible(false); m_Controls.UseImageSpacing->setToolTip("Image spacing of y-Axis must be in us, x-Axis in mm."); m_Controls.UseImageSpacing->setToolTipDuration(5000); m_Controls.ProgressInfo->setVisible(false); m_Controls.UseGPUBmode->hide(); #ifndef PHOTOACOUSTICS_USE_GPU m_Controls.UseGPUBf->setEnabled(false); m_Controls.UseGPUBf->setChecked(false); m_Controls.UseGPUBmode->setEnabled(false); m_Controls.UseGPUBmode->setChecked(false); #endif UseImageSpacing(); + ChangedProbe(); +} + +void PAImageProcessing::ChangedProbe() +{ + if (m_Controls.Geometry->currentText() == "Concave") + { + m_Controls.ProbeRadius->setEnabled(true); + } + else + { + m_Controls.ProbeRadius->setEnabled(false); + } } void PAImageProcessing::UseSignalDelay() { if (m_Controls.UseSignalDelay->isChecked()) { m_Controls.SignalDelay->setEnabled(true); } else { m_Controls.SignalDelay->setEnabled(false); } } void PAImageProcessing::ChangedSOSBandpass() { m_Controls.SpeedOfSound->setValue(m_Controls.BPSpeedOfSound->value()); } void PAImageProcessing::ChangedSOSBeamforming() { m_Controls.BPSpeedOfSound->setValue(m_Controls.SpeedOfSound->value()); } std::vector splitpath( const std::string& str , const std::set delimiters) { std::vector result; char const* pch = str.c_str(); char const* start = pch; for (; *pch; ++pch) { if (delimiters.find(*pch) != delimiters.end()) { if (start != pch) { std::string str(start, pch); result.push_back(str); } else { result.push_back(""); } start = pch + 1; } } result.push_back(start); return result; } void PAImageProcessing::UpdateSaveBoxes() { if (m_Controls.StepBeamforming->isChecked()) m_Controls.SaveBeamforming->setEnabled(true); else m_Controls.SaveBeamforming->setEnabled(false); if (m_Controls.StepCropping->isChecked()) m_Controls.SaveCropping->setEnabled(true); else m_Controls.SaveCropping->setEnabled(false); if (m_Controls.StepBandpass->isChecked()) m_Controls.SaveBandpass->setEnabled(true); else m_Controls.SaveBandpass->setEnabled(false); if (m_Controls.StepBMode->isChecked()) m_Controls.SaveBMode->setEnabled(true); else m_Controls.SaveBMode->setEnabled(false); } void PAImageProcessing::BatchProcessing() { QFileDialog LoadDialog(nullptr, "Select Files to be processed"); LoadDialog.setFileMode(QFileDialog::FileMode::ExistingFiles); LoadDialog.setNameFilter(tr("Images (*.nrrd)")); LoadDialog.setViewMode(QFileDialog::Detail); QStringList fileNames; if (LoadDialog.exec()) fileNames = LoadDialog.selectedFiles(); QString saveDir = QFileDialog::getExistingDirectory(nullptr, tr("Select Directory To Save To"), "", QFileDialog::ShowDirsOnly | QFileDialog::DontResolveSymlinks); DisableControls(); std::set delims{ '/' }; bool doSteps[] = { m_Controls.StepBeamforming->isChecked(), m_Controls.StepCropping->isChecked() , m_Controls.StepBandpass->isChecked(), m_Controls.StepBMode->isChecked() }; bool saveSteps[] = { m_Controls.SaveBeamforming->isChecked(), m_Controls.SaveCropping->isChecked() , m_Controls.SaveBandpass->isChecked(), m_Controls.SaveBMode->isChecked() }; for (int fileNumber = 0; fileNumber < fileNames.size(); ++fileNumber) { m_Controls.progressBar->setValue(0); m_Controls.progressBar->setVisible(true); m_Controls.ProgressInfo->setVisible(true); m_Controls.ProgressInfo->setText("loading file"); QString filename = fileNames.at(fileNumber); auto split = splitpath(filename.toStdString(), delims); std::string imageName = split.at(split.size() - 1); // remove ".nrrd" imageName = imageName.substr(0, imageName.size() - 5); mitk::Image::Pointer image = mitk::IOUtil::Load(filename.toStdString().c_str()); auto BFconfig = CreateBeamformingSettings(image); // Beamforming if (doSteps[0]) { if (m_Controls.UseSignalDelay->isChecked()) { float signalDelay = m_Controls.SignalDelay->value(); if (signalDelay != 0) { int cropPixels = std::round(signalDelay / BFconfig->GetTimeSpacing() / 1000000); MITK_INFO << cropPixels; int errCode = 0; image = m_FilterBank->ApplyCropping(image, cropPixels, 0, 0, 0, 0, 0, &errCode); if (errCode == -1) { QMessageBox Msgbox; Msgbox.setText("It has been attempted to cut off more pixels than the image contains. Aborting batch processing."); Msgbox.exec(); m_Controls.progressBar->setVisible(false); EnableControls(); return; } BFconfig = mitk::BeamformingSettings::New(BFconfig->GetPitchInMeters(), BFconfig->GetSpeedOfSound(), BFconfig->GetTimeSpacing(), BFconfig->GetAngle(), BFconfig->GetIsPhotoacousticImage(), BFconfig->GetSamplesPerLine(), BFconfig->GetReconstructionLines(), image->GetDimensions(), BFconfig->GetReconstructionDepth(), - BFconfig->GetUseGPU(), BFconfig->GetGPUBatchSize(), BFconfig->GetDelayCalculationMethod(), BFconfig->GetApod(), - BFconfig->GetApodizationArraySize(), BFconfig->GetAlgorithm()); + BFconfig->GetUseGPU(), BFconfig->GetGPUBatchSize(), BFconfig->GetApod(), + BFconfig->GetApodizationArraySize(), BFconfig->GetAlgorithm(), BFconfig->GetGeometry(), BFconfig->GetProbeRadius()); } } std::function progressHandle = [this](int progress, std::string progressInfo) { this->UpdateProgress(progress, progressInfo); }; m_Controls.progressBar->setValue(100); image = m_FilterBank->ApplyBeamforming(image, BFconfig, progressHandle); if (saveSteps[0]) { std::string saveFileName = saveDir.toStdString() + "/" + imageName + " beamformed" + ".nrrd"; mitk::IOUtil::Save(image, saveFileName); } } // Cropping if (doSteps[1]) { m_Controls.ProgressInfo->setText("cropping image"); int errCode = 0; - image = m_FilterBank->ApplyCropping(image, m_Controls.CutoffAbove->value(), m_Controls.CutoffBelow->value(), 0, 0, 0, 0, &errCode); + image = m_FilterBank->ApplyCropping(image, m_Controls.CutoffAbove->value(), m_Controls.CutoffBelow->value(), m_Controls.CutoffRight->value(), m_Controls.CutoffLeft->value(), 0, 0, &errCode); if (errCode == -1) { QMessageBox Msgbox; Msgbox.setText("It has been attempted to cut off more pixels than the image contains. Aborting batch processing."); Msgbox.exec(); m_Controls.progressBar->setVisible(false); EnableControls(); return; } if (saveSteps[1]) { std::string saveFileName = saveDir.toStdString() + "/" + imageName + " cropped" + ".nrrd"; mitk::IOUtil::Save(image, saveFileName); } } // Bandpass if (doSteps[2]) { m_Controls.ProgressInfo->setText("applying bandpass"); float recordTime = image->GetDimension(1)*image->GetGeometry()->GetSpacing()[1] / 1000 / m_Controls.BPSpeedOfSound->value(); // add a safeguard so the program does not chrash when applying a Bandpass that reaches out of the bounds of the image float maxFrequency = 1 / (recordTime / image->GetDimension(1)) * image->GetDimension(1) / 2 / 2 / 1000; float BPHighPass = 1000000 * m_Controls.BPhigh->value(); // [Hz] float BPLowPass = maxFrequency - 1000000 * m_Controls.BPlow->value(); // [Hz] if (BPLowPass > maxFrequency) { QMessageBox Msgbox; Msgbox.setText("LowPass too low, disabled it."); Msgbox.exec(); BPLowPass = 0; } if (BPLowPass < 0) { QMessageBox Msgbox; Msgbox.setText("LowPass too high, disabled it."); Msgbox.exec(); BPLowPass = 0; } if (BPHighPass > maxFrequency) { QMessageBox Msgbox; Msgbox.setText("HighPass too high, disabled it."); Msgbox.exec(); BPHighPass = 0; } if (BPHighPass > maxFrequency - BPLowPass) { QMessageBox Msgbox; Msgbox.setText("HighPass higher than LowPass, disabled both."); Msgbox.exec(); BPHighPass = 0; BPLowPass = 0; } image = m_FilterBank->ApplyBandpassFilter(image, BPHighPass, BPLowPass, m_Controls.BPFalloffHigh->value(), m_Controls.BPFalloffLow->value(), BFconfig->GetTimeSpacing(), BFconfig->GetSpeedOfSound(), m_Controls.IsBFImage->isChecked()); if (saveSteps[2]) { std::string saveFileName = saveDir.toStdString() + "/" + imageName + " bandpassed" + ".nrrd"; mitk::IOUtil::Save(image, saveFileName); } } // Bmode if (doSteps[3]) { m_Controls.ProgressInfo->setText("applying bmode filter"); if (m_Controls.BModeMethod->currentText() == "Absolute Filter") image = m_FilterBank->ApplyBmodeFilter(image, mitk::PhotoacousticFilterService::BModeMethod::Abs, m_UseLogfilter); else if (m_Controls.BModeMethod->currentText() == "Envelope Detection") image = m_FilterBank->ApplyBmodeFilter(image, mitk::PhotoacousticFilterService::BModeMethod::EnvelopeDetection, m_UseLogfilter); if (m_ResampleSpacing != 0) { double desiredSpacing[2]{ image->GetGeometry()->GetSpacing()[0], m_ResampleSpacing }; image = m_FilterBank->ApplyResampling(image, desiredSpacing); } if (saveSteps[3]) { std::string saveFileName = saveDir.toStdString() + "/" + imageName + " bmode" + ".nrrd"; mitk::IOUtil::Save(image, saveFileName); } } m_Controls.progressBar->setVisible(false); } EnableControls(); } void PAImageProcessing::StartBeamformingThread() { QList nodes = this->GetDataManagerSelection(); if (nodes.empty()) return; mitk::DataStorage::Pointer storage = this->GetDataStorage(); mitk::DataNode::Pointer node = nodes.front(); if (!node) { // Nothing selected. Inform the user and return QMessageBox::information(NULL, "Template", "Please load and select an image before starting image processing."); return; } mitk::BaseData* data = node->GetData(); if (data) { // test if this data item is an image or not (could also be a surface or something totally different) mitk::Image* image = dynamic_cast(data); if (image) { auto BFconfig = CreateBeamformingSettings(image); std::stringstream message; std::string name; message << "Performing beamforming for image "; if (node->GetName(name)) { // a property called "name" was found for this DataNode message << "'" << name << "'"; m_OldNodeName = name; } else m_OldNodeName = " "; message << "."; MITK_INFO << message.str(); m_Controls.progressBar->setValue(0); m_Controls.progressBar->setVisible(true); m_Controls.ProgressInfo->setVisible(true); m_Controls.ProgressInfo->setText("started"); DisableControls(); BeamformingThread *thread = new BeamformingThread(); connect(thread, &BeamformingThread::result, this, &PAImageProcessing::HandleResults); connect(thread, &BeamformingThread::updateProgress, this, &PAImageProcessing::UpdateProgress); connect(thread, &BeamformingThread::finished, thread, &QObject::deleteLater); thread->setConfig(BFconfig); if (m_Controls.UseSignalDelay->isChecked()) thread->setSignalDelay(m_Controls.SignalDelay->value()); thread->setInputImage(image); thread->setFilterBank(m_FilterBank); MITK_INFO << "Started new thread for Beamforming"; thread->start(); } } } void PAImageProcessing::HandleResults(mitk::Image::Pointer image, std::string nameExtension) { if (image == nullptr) { QMessageBox Msgbox; Msgbox.setText("An error has occurred during processing; please see the console output."); Msgbox.exec(); // disable progress bar m_Controls.progressBar->setVisible(false); m_Controls.ProgressInfo->setVisible(false); EnableControls(); return; } MITK_INFO << "Handling results..."; auto newNode = mitk::DataNode::New(); newNode->SetData(image); newNode->SetName(m_OldNodeName + nameExtension); // update level window for the current dynamic range mitk::LevelWindow levelWindow; newNode->GetLevelWindow(levelWindow); levelWindow.SetAuto(image, true, true); newNode->SetLevelWindow(levelWindow); // add new node to data storage this->GetDataStorage()->Add(newNode); // disable progress bar m_Controls.progressBar->setVisible(false); m_Controls.ProgressInfo->setVisible(false); EnableControls(); // update rendering mitk::RenderingManager::GetInstance()->InitializeViews(image->GetGeometry(), mitk::RenderingManager::REQUEST_UPDATE_ALL, true); mitk::RenderingManager::GetInstance()->RequestUpdateAll(); MITK_INFO << "Handling results...[Done]"; } void PAImageProcessing::StartBmodeThread() { QList nodes = this->GetDataManagerSelection(); if (nodes.empty()) return; mitk::DataStorage::Pointer storage = this->GetDataStorage(); mitk::DataNode::Pointer node = nodes.front(); if (!node) { // Nothing selected. Inform the user and return QMessageBox::information(NULL, "Template", "Please load and select an image before starting image processing."); return; } mitk::BaseData* data = node->GetData(); if (data) { // test if this data item is an image or not (could also be a surface or something totally different) mitk::Image* image = dynamic_cast(data); if (image) { std::stringstream message; std::string name; message << "Performing image processing for image "; if (node->GetName(name)) { // a property called "name" was found for this DataNode message << "'" << name << "'"; m_OldNodeName = name; } else m_OldNodeName = " "; message << "."; MITK_INFO << message.str(); DisableControls(); BmodeThread *thread = new BmodeThread(); connect(thread, &BmodeThread::result, this, &PAImageProcessing::HandleResults); connect(thread, &BmodeThread::finished, thread, &QObject::deleteLater); bool useGPU = m_Controls.UseGPUBmode->isChecked(); if (m_Controls.BModeMethod->currentText() == "Absolute Filter") thread->setConfig(m_UseLogfilter, m_ResampleSpacing, mitk::PhotoacousticFilterService::BModeMethod::Abs, useGPU); else if (m_Controls.BModeMethod->currentText() == "Envelope Detection") thread->setConfig(m_UseLogfilter, m_ResampleSpacing, mitk::PhotoacousticFilterService::BModeMethod::EnvelopeDetection, useGPU); thread->setInputImage(image); thread->setFilterBank(m_FilterBank); MITK_INFO << "Started new thread for Image Processing"; thread->start(); } } } void PAImageProcessing::StartCropThread() { QList nodes = this->GetDataManagerSelection(); if (nodes.empty()) return; mitk::DataStorage::Pointer storage = this->GetDataStorage(); mitk::DataNode::Pointer node = nodes.front(); if (!node) { // Nothing selected. Inform the user and return QMessageBox::information(NULL, "Template", "Please load and select an image before starting image cropping."); return; } mitk::BaseData* data = node->GetData(); if (data) { // test if this data item is an image or not (could also be a surface or something totally different) mitk::Image* image = dynamic_cast(data); if (image) { std::stringstream message; std::string name; message << "Performing image cropping for image "; if (node->GetName(name)) { // a property called "name" was found for this DataNode message << "'" << name << "'"; m_OldNodeName = name; } else m_OldNodeName = " "; message << "."; MITK_INFO << message.str(); DisableControls(); CropThread *thread = new CropThread(); connect(thread, &CropThread::result, this, &PAImageProcessing::HandleResults); connect(thread, &CropThread::finished, thread, &QObject::deleteLater); if(m_Controls.Partial->isChecked()) - thread->setConfig(m_Controls.CutoffAbove->value(), m_Controls.CutoffBelow->value(), m_Controls.boundLow->value(), m_Controls.boundHigh->value()); + thread->setConfig(m_Controls.CutoffAbove->value(), m_Controls.CutoffBelow->value(), m_Controls.CutoffRight->value(), m_Controls.CutoffLeft->value(), m_Controls.boundLow->value(), m_Controls.boundHigh->value()); else - thread->setConfig(m_Controls.CutoffAbove->value(), m_Controls.CutoffBelow->value(), 0, image->GetDimension(2) - 1); + thread->setConfig(m_Controls.CutoffAbove->value(), m_Controls.CutoffBelow->value(), m_Controls.CutoffRight->value(), m_Controls.CutoffLeft->value(), 0, image->GetDimension(2) - 1); thread->setInputImage(image); thread->setFilterBank(m_FilterBank); MITK_INFO << "Started new thread for Image Cropping"; thread->start(); } } } void PAImageProcessing::StartBandpassThread() { QList nodes = this->GetDataManagerSelection(); if (nodes.empty()) return; mitk::DataStorage::Pointer storage = this->GetDataStorage(); mitk::DataNode::Pointer node = nodes.front(); if (!node) { // Nothing selected. Inform the user and return QMessageBox::information(NULL, "Template", "Please load and select an image before applying a bandpass filter."); return; } mitk::BaseData* data = node->GetData(); if (data) { // test if this data item is an image or not (could also be a surface or something totally different) mitk::Image* image = dynamic_cast(data); if (image) { auto config = CreateBeamformingSettings(image); std::stringstream message; std::string name; message << "Performing Bandpass filter on image "; if (node->GetName(name)) { // a property called "name" was found for this DataNode message << "'" << name << "'"; m_OldNodeName = name; } else m_OldNodeName = " "; message << "."; MITK_INFO << message.str(); DisableControls(); BandpassThread *thread = new BandpassThread(); connect(thread, &BandpassThread::result, this, &PAImageProcessing::HandleResults); connect(thread, &BandpassThread::finished, thread, &QObject::deleteLater); float BPHighPass = 1000000.0f * m_Controls.BPhigh->value(); // [Now in Hz] float BPLowPass = 1000000.0f * m_Controls.BPlow->value(); // [Now in Hz] thread->setConfig(BPHighPass, BPLowPass, m_Controls.BPFalloffLow->value(), m_Controls.BPFalloffHigh->value(), config->GetTimeSpacing(), config->GetSpeedOfSound(), m_Controls.IsBFImage->isChecked()); thread->setInputImage(image); thread->setFilterBank(m_FilterBank); MITK_INFO << "Started new thread for Bandpass filter"; thread->start(); } } } void PAImageProcessing::SliceBoundsEnabled() { if (!m_Controls.Partial->isChecked()) { m_Controls.boundLow->setEnabled(false); m_Controls.boundHigh->setEnabled(false); return; } else { m_Controls.boundLow->setEnabled(true); m_Controls.boundHigh->setEnabled(true); } } void PAImageProcessing::UpperSliceBoundChanged() { if (m_Controls.boundLow->value() > m_Controls.boundHigh->value()) { m_Controls.boundLow->setValue(m_Controls.boundHigh->value()); } } void PAImageProcessing::LowerSliceBoundChanged() { if (m_Controls.boundLow->value() > m_Controls.boundHigh->value()) { m_Controls.boundHigh->setValue(m_Controls.boundLow->value()); } } void PAImageProcessing::UpdateProgress(int progress, std::string progressInfo) { if (progress < 100) m_Controls.progressBar->setValue(progress); else m_Controls.progressBar->setValue(100); m_Controls.ProgressInfo->setText(progressInfo.c_str()); qApp->processEvents(); } void PAImageProcessing::PAMessageBox(std::string message) { if (0 != message.compare("noMessage")) { QMessageBox msgBox; msgBox.setText(message.c_str()); msgBox.exec(); } } void PAImageProcessing::UpdateImageInfo() { QList nodes = this->GetDataManagerSelection(); if (nodes.empty()) return; mitk::DataNode::Pointer node = nodes.front(); if (!node) { // Nothing selected return; } mitk::BaseData* data = node->GetData(); if (data) { // test if this data item is an image or not (could also be a surface or something totally different) mitk::Image* image = dynamic_cast(data); if (image) { // beamforming configs if (m_Controls.UseImageSpacing->isChecked()) { m_Controls.ElementCount->setValue(image->GetDimension(0)); m_Controls.Pitch->setValue(image->GetGeometry()->GetSpacing()[0]); } m_Controls.boundLow->setMaximum(image->GetDimension(2) - 1); m_Controls.boundHigh->setMaximum(image->GetDimension(2) - 1); float speedOfSound = m_Controls.SpeedOfSound->value(); // [m/s] std::stringstream frequency; float timeSpacing; if (m_Controls.UseImageSpacing->isChecked()) { timeSpacing = image->GetGeometry()->GetSpacing()[1] / 1000000.0f; MITK_INFO << "Calculated Scan Depth of " << (image->GetDimension(1)*image->GetGeometry()->GetSpacing()[1] / 1000000) * speedOfSound * 100 / 2 << "cm"; } else { timeSpacing = (2 * m_Controls.ScanDepth->value() / 1000 / speedOfSound) / image->GetDimension(1); } float maxFrequency = (1 / timeSpacing) / 2; if(m_Controls.IsBFImage->isChecked()) maxFrequency = ( 1 / (image->GetGeometry()->GetSpacing()[1] / 1e3 / speedOfSound)) / 2; frequency << maxFrequency / 1e6; //[MHz] frequency << "MHz"; m_Controls.BPhigh->setMaximum(maxFrequency / 1e6); m_Controls.BPlow->setMaximum(maxFrequency / 1e6); frequency << " is the maximal allowed frequency for the selected image."; m_Controls.BPhigh->setToolTip(frequency.str().c_str()); m_Controls.BPlow->setToolTip(frequency.str().c_str()); m_Controls.BPhigh->setToolTipDuration(5000); m_Controls.BPlow->setToolTipDuration(5000); } } } void PAImageProcessing::OnSelectionChanged(berry::IWorkbenchPart::Pointer /*source*/, const QList& nodes) { // iterate all selected objects, adjust warning visibility foreach(mitk::DataNode::Pointer node, nodes) { if (node.IsNotNull() && dynamic_cast(node->GetData())) { m_Controls.labelWarning->setVisible(false); m_Controls.buttonApplyBModeFilter->setEnabled(true); m_Controls.labelWarning2->setVisible(false); m_Controls.buttonApplyCropFilter->setEnabled(true); m_Controls.labelWarning3->setVisible(false); m_Controls.buttonApplyBandpass->setEnabled(true); m_Controls.labelWarning4->setVisible(false); m_Controls.buttonApplyBeamforming->setEnabled(true); UpdateImageInfo(); return; } } m_Controls.labelWarning->setVisible(true); m_Controls.buttonApplyBModeFilter->setEnabled(false); m_Controls.labelWarning2->setVisible(true); m_Controls.buttonApplyCropFilter->setEnabled(false); m_Controls.labelWarning3->setVisible(true); m_Controls.buttonApplyBandpass->setEnabled(false); m_Controls.labelWarning4->setVisible(true); m_Controls.buttonApplyBeamforming->setEnabled(false); } void PAImageProcessing::UseResampling() { if (m_Controls.DoResampling->isChecked()) { m_Controls.ResamplingValue->setEnabled(true); m_ResampleSpacing = m_Controls.ResamplingValue->value(); } else { m_Controls.ResamplingValue->setEnabled(false); m_ResampleSpacing = 0; } } void PAImageProcessing::UseLogfilter() { m_UseLogfilter = m_Controls.Logfilter->isChecked(); } void PAImageProcessing::SetResampling() { m_ResampleSpacing = m_Controls.ResamplingValue->value(); } mitk::BeamformingSettings::Pointer PAImageProcessing::CreateBeamformingSettings(mitk::Image::Pointer image) { mitk::BeamformingSettings::BeamformingAlgorithm algorithm = mitk::BeamformingSettings::BeamformingAlgorithm::DAS; if ("DAS" == m_Controls.BFAlgorithm->currentText()) algorithm = mitk::BeamformingSettings::BeamformingAlgorithm::DAS; else if ("DMAS" == m_Controls.BFAlgorithm->currentText()) algorithm = mitk::BeamformingSettings::BeamformingAlgorithm::DMAS; else if ("sDMAS" == m_Controls.BFAlgorithm->currentText()) algorithm = mitk::BeamformingSettings::BeamformingAlgorithm::sDMAS; - mitk::BeamformingSettings::DelayCalc delay = mitk::BeamformingSettings::DelayCalc::Spherical; - mitk::BeamformingSettings::Apodization apod = mitk::BeamformingSettings::Apodization::Box; if ("Von Hann" == m_Controls.Apodization->currentText()) { apod = mitk::BeamformingSettings::Apodization::Hann; } else if ("Hamming" == m_Controls.Apodization->currentText()) { apod = mitk::BeamformingSettings::Apodization::Hamm; } else if ("Box" == m_Controls.Apodization->currentText()) { apod = mitk::BeamformingSettings::Apodization::Box; } float pitchInMeters = m_Controls.Pitch->value() / 1000; // [m] float speedOfSound = m_Controls.SpeedOfSound->value(); // [m/s] unsigned int samplesPerLine = m_Controls.Samples->value(); unsigned int reconstructionLines = m_Controls.Lines->value(); unsigned int apodizatonArraySize = m_Controls.Lines->value(); float angle = m_Controls.Angle->value(); // [deg] bool useGPU = m_Controls.UseGPUBf->isChecked(); float timeSpacing; if (m_Controls.UseImageSpacing->isChecked()) { timeSpacing = image->GetGeometry()->GetSpacing()[1] / 1000000.0f; MITK_INFO << "Calculated Scan Depth of " << (image->GetDimension(1)*image->GetGeometry()->GetSpacing()[1] / 1000000) * speedOfSound * 100 / 2 << "cm"; } else { timeSpacing = (2 * m_Controls.ScanDepth->value() / 1000 / speedOfSound) / image->GetDimension(1); } bool isPAImage = true; if ("US Image" == m_Controls.ImageType->currentText()) { isPAImage = false; } else if ("PA Image" == m_Controls.ImageType->currentText()) { isPAImage = true; } float reconstructionDepth = m_Controls.ReconstructionDepth->value() / 1000.f; // [m] + mitk::BeamformingSettings::ProbeGeometry geometry = mitk::BeamformingSettings::ProbeGeometry::Linear; + if ("Linear" == m_Controls.Geometry->currentText()) + { + geometry = mitk::BeamformingSettings::ProbeGeometry::Linear; + } + else if ("Concave" == m_Controls.Geometry->currentText()) + { + geometry = mitk::BeamformingSettings::ProbeGeometry::Concave; + } + float probeRadius = m_Controls.ProbeRadius->value()/1000.f; // [m] + return mitk::BeamformingSettings::New(pitchInMeters, speedOfSound, timeSpacing, angle, isPAImage, samplesPerLine, reconstructionLines, - image->GetDimensions(), reconstructionDepth, useGPU, GPU_BATCH_SIZE, delay, apod, - apodizatonArraySize, algorithm); + image->GetDimensions(), reconstructionDepth, useGPU, GPU_BATCH_SIZE, apod, + apodizatonArraySize, algorithm, geometry, probeRadius); } void PAImageProcessing::EnableControls() { m_Controls.BatchProcessing->setEnabled(true); m_Controls.StepBeamforming->setEnabled(true); m_Controls.StepBandpass->setEnabled(true); m_Controls.StepCropping->setEnabled(true); m_Controls.StepBMode->setEnabled(true); UpdateSaveBoxes(); m_Controls.DoResampling->setEnabled(true); UseResampling(); m_Controls.Logfilter->setEnabled(true); m_Controls.BModeMethod->setEnabled(true); m_Controls.buttonApplyBModeFilter->setEnabled(true); m_Controls.CutoffAbove->setEnabled(true); m_Controls.CutoffBelow->setEnabled(true); m_Controls.buttonApplyCropFilter->setEnabled(true); m_Controls.BPSpeedOfSound->setEnabled(true); m_Controls.buttonApplyBandpass->setEnabled(true); m_Controls.Partial->setEnabled(true); m_Controls.boundHigh->setEnabled(true); m_Controls.boundLow->setEnabled(true); m_Controls.BFAlgorithm->setEnabled(true); m_Controls.ReconstructionDepth->setEnabled(true); m_Controls.ImageType->setEnabled(true); m_Controls.Apodization->setEnabled(true); #ifdef PHOTOACOUSTICS_USE_GPU m_Controls.UseGPUBf->setEnabled(true); m_Controls.UseGPUBmode->setEnabled(true); #endif m_Controls.BPhigh->setEnabled(true); m_Controls.BPlow->setEnabled(true); m_Controls.BPFalloffLow->setEnabled(true); m_Controls.BPFalloffHigh->setEnabled(true); m_Controls.UseImageSpacing->setEnabled(true); UseImageSpacing(); m_Controls.Pitch->setEnabled(true); m_Controls.ElementCount->setEnabled(true); m_Controls.SpeedOfSound->setEnabled(true); m_Controls.Samples->setEnabled(true); m_Controls.Lines->setEnabled(true); m_Controls.Angle->setEnabled(true); m_Controls.buttonApplyBeamforming->setEnabled(true); m_Controls.UseSignalDelay->setEnabled(true); m_Controls.SignalDelay->setEnabled(true); } void PAImageProcessing::DisableControls() { m_Controls.BatchProcessing->setEnabled(false); m_Controls.StepBeamforming->setEnabled(false); m_Controls.StepBandpass->setEnabled(false); m_Controls.StepCropping->setEnabled(false); m_Controls.StepBMode->setEnabled(false); m_Controls.SaveBeamforming->setEnabled(false); m_Controls.SaveBandpass->setEnabled(false); m_Controls.SaveCropping->setEnabled(false); m_Controls.SaveBMode->setEnabled(false); m_Controls.DoResampling->setEnabled(false); m_Controls.ResamplingValue->setEnabled(false); m_Controls.Logfilter->setEnabled(false); m_Controls.BModeMethod->setEnabled(false); m_Controls.buttonApplyBModeFilter->setEnabled(false); m_Controls.CutoffAbove->setEnabled(false); m_Controls.CutoffBelow->setEnabled(false); m_Controls.buttonApplyCropFilter->setEnabled(false); m_Controls.BPSpeedOfSound->setEnabled(false); m_Controls.buttonApplyBandpass->setEnabled(false); m_Controls.Partial->setEnabled(false); m_Controls.boundHigh->setEnabled(false); m_Controls.boundLow->setEnabled(false); m_Controls.BFAlgorithm->setEnabled(false); m_Controls.ReconstructionDepth->setEnabled(false); m_Controls.ImageType->setEnabled(false); m_Controls.Apodization->setEnabled(false); #ifdef PHOTOACOUSTICS_USE_GPU m_Controls.UseGPUBf->setEnabled(false); m_Controls.UseGPUBmode->setEnabled(false); #endif m_Controls.BPhigh->setEnabled(false); m_Controls.BPlow->setEnabled(false); m_Controls.BPFalloffLow->setEnabled(false); m_Controls.BPFalloffHigh->setEnabled(false); m_Controls.UseImageSpacing->setEnabled(false); m_Controls.ScanDepth->setEnabled(false); m_Controls.Pitch->setEnabled(false); m_Controls.ElementCount->setEnabled(false); m_Controls.SpeedOfSound->setEnabled(false); m_Controls.Samples->setEnabled(false); m_Controls.Lines->setEnabled(false); m_Controls.Angle->setEnabled(false); m_Controls.buttonApplyBeamforming->setEnabled(false); m_Controls.UseSignalDelay->setEnabled(false); m_Controls.SignalDelay->setEnabled(false); } void PAImageProcessing::UseImageSpacing() { if (m_Controls.UseImageSpacing->isChecked()) { m_Controls.ScanDepth->setDisabled(true); } else { m_Controls.ScanDepth->setEnabled(true); } } #include void BeamformingThread::run() { if (m_SignalDelay != 0) { int cropPixels = std::round(m_SignalDelay / m_BFconfig->GetTimeSpacing() / 1000000); MITK_INFO << cropPixels; int errCode = 0; m_InputImage = m_FilterBank->ApplyCropping(m_InputImage, cropPixels, 0, 0, 0, 0, 0, &errCode); m_BFconfig = mitk::BeamformingSettings::New(m_BFconfig->GetPitchInMeters(), m_BFconfig->GetSpeedOfSound(), m_BFconfig->GetTimeSpacing(), m_BFconfig->GetAngle(), m_BFconfig->GetIsPhotoacousticImage(), m_BFconfig->GetSamplesPerLine(), m_BFconfig->GetReconstructionLines(), m_InputImage->GetDimensions(), m_BFconfig->GetReconstructionDepth(), - m_BFconfig->GetUseGPU(), m_BFconfig->GetGPUBatchSize(), m_BFconfig->GetDelayCalculationMethod(), m_BFconfig->GetApod(), - m_BFconfig->GetApodizationArraySize(), m_BFconfig->GetAlgorithm()); + m_BFconfig->GetUseGPU(), m_BFconfig->GetGPUBatchSize(), m_BFconfig->GetApod(), + m_BFconfig->GetApodizationArraySize(), m_BFconfig->GetAlgorithm(), m_BFconfig->GetGeometry(), m_BFconfig->GetProbeRadius()); } mitk::Image::Pointer resultImage; std::function progressHandle = [this](int progress, std::string progressInfo) { emit updateProgress(progress, progressInfo); }; resultImage = m_FilterBank->ApplyBeamforming(m_InputImage, m_BFconfig, progressHandle); emit result(resultImage, "_bf"); } void BeamformingThread::setConfig(mitk::BeamformingSettings::Pointer BFconfig) { m_BFconfig = BFconfig; } void BeamformingThread::setSignalDelay(float delay) { m_SignalDelay = delay; } void BeamformingThread::setInputImage(mitk::Image::Pointer image) { m_InputImage = image; } void BmodeThread::run() { mitk::Image::Pointer resultImage = m_FilterBank->ApplyBmodeFilter(m_InputImage, m_Method, m_UseLogfilter); if (m_ResampleSpacing != 0) { double desiredSpacing[2]{ m_InputImage->GetGeometry()->GetSpacing()[0], m_ResampleSpacing }; resultImage = m_FilterBank->ApplyResampling(resultImage, desiredSpacing); } emit result(resultImage, "_bmode"); } void BmodeThread::setConfig(bool useLogfilter, double resampleSpacing, mitk::PhotoacousticFilterService::BModeMethod method, bool useGPU) { m_UseLogfilter = useLogfilter; m_ResampleSpacing = resampleSpacing; m_Method = method; m_UseGPU = useGPU; } void BmodeThread::setInputImage(mitk::Image::Pointer image) { m_InputImage = image; } void CropThread::run() { mitk::Image::Pointer resultImage; - int errCode = 0; - resultImage = m_FilterBank->ApplyCropping(m_InputImage, m_CutAbove, m_CutBelow, 0, 0, m_CutSliceFirst, (m_InputImage->GetDimension(2) - 1) - m_CutSliceLast, &errCode); + resultImage = m_FilterBank->ApplyCropping(m_InputImage, m_CutAbove, m_CutBelow, m_CutRight, m_CutLeft, m_CutSliceFirst, (m_InputImage->GetDimension(2) - 1) - m_CutSliceLast, &errCode); if (errCode == -1) { emit result(nullptr, "_cropped"); return; } emit result(resultImage, "_cropped"); } -void CropThread::setConfig(unsigned int CutAbove, unsigned int CutBelow, unsigned int CutSliceFirst, unsigned int CutSliceLast) +void CropThread::setConfig(unsigned int CutAbove, unsigned int CutBelow, unsigned int CutRight, unsigned int CutLeft, unsigned int CutSliceFirst, unsigned int CutSliceLast) { m_CutAbove = CutAbove; m_CutBelow = CutBelow; + m_CutRight = CutRight; + m_CutLeft = CutLeft; + m_CutSliceLast = CutSliceLast; m_CutSliceFirst = CutSliceFirst; } void CropThread::setInputImage(mitk::Image::Pointer image) { m_InputImage = image; } void BandpassThread::run() { mitk::Image::Pointer resultImage = m_FilterBank->ApplyBandpassFilter(m_InputImage, m_BPHighPass, m_BPLowPass, m_TukeyAlphaHighPass, m_TukeyAlphaLowPass, m_TimeSpacing, m_SpeedOfSound, m_IsBFImage); emit result(resultImage, "_bandpassed"); } void BandpassThread::setConfig(float BPHighPass, float BPLowPass, float TukeyAlphaHighPass, float TukeyAlphaLowPass, float TimeSpacing, float SpeedOfSound, bool IsBFImage) { m_BPHighPass = BPHighPass; m_BPLowPass = BPLowPass; m_TukeyAlphaHighPass = TukeyAlphaHighPass; m_TukeyAlphaLowPass = TukeyAlphaLowPass; m_TimeSpacing = TimeSpacing; m_SpeedOfSound = SpeedOfSound; m_IsBFImage = IsBFImage; } void BandpassThread::setInputImage(mitk::Image::Pointer image) { m_InputImage = image; } diff --git a/Plugins/org.mitk.gui.qt.photoacoustics.imageprocessing/src/internal/PAImageProcessing.h b/Plugins/org.mitk.gui.qt.photoacoustics.imageprocessing/src/internal/PAImageProcessing.h index 159ae2f2be..95d9e878f3 100644 --- a/Plugins/org.mitk.gui.qt.photoacoustics.imageprocessing/src/internal/PAImageProcessing.h +++ b/Plugins/org.mitk.gui.qt.photoacoustics.imageprocessing/src/internal/PAImageProcessing.h @@ -1,244 +1,247 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef PAImageProcessing_h #define PAImageProcessing_h #include #include #include #include #include "ui_PAImageProcessingControls.h" #include "mitkBeamformingFilter.h" #include "mitkBeamformingSettings.h" Q_DECLARE_METATYPE(mitk::Image::Pointer) Q_DECLARE_METATYPE(std::string) /*! * \brief Plugin implementing an interface for the Photoacoustic Algorithms Module * * Beamforming, Image processing as B-Mode filtering, cropping, resampling, as well as batch processing can be performed using this plugin. */ class PAImageProcessing : public QmitkAbstractView { // this is needed for all Qt objects that should have a Qt meta-object // (everything that derives from QObject and wants to have signal/slots) Q_OBJECT public: static const std::string VIEW_ID; PAImageProcessing(); protected slots: void UpperSliceBoundChanged(); void LowerSliceBoundChanged(); void SliceBoundsEnabled(); + void ChangedProbe(); void UseResampling(); void UseLogfilter(); void SetResampling(); void UseImageSpacing(); void UpdateImageInfo(); void UseSignalDelay(); /** \brief Beamforming is being performed in a separate thread to keep the workbench from freezing. */ void StartBeamformingThread(); /** \brief B-mode filtering is being performed in a separate thread to keep the workbench from freezing. */ void StartBmodeThread(); /** \brief Cropping is being performed in a separate thread to keep the workbench from freezing. */ void StartCropThread(); /** \brief Method called when the bandpass thread finishes; * it adds the image to a new data node and registers it to the worbench's data storage */ void HandleResults(mitk::Image::Pointer image, std::string nameExtension); /** \brief Bandpassing is being performed in a separate thread to keep the workbench from freezing. */ void StartBandpassThread(); void UpdateProgress(int progress, std::string progressInfo); void PAMessageBox(std::string message); void BatchProcessing(); void UpdateSaveBoxes(); void ChangedSOSBandpass(); void ChangedSOSBeamforming(); protected: virtual void CreateQtPartControl(QWidget *parent) override; virtual void SetFocus() override; /** \brief called by QmitkFunctionality when DataManager's selection has changed. * On a change some parameters are internally updated to calculate bounds for GUI elements as the slice selector for beamforming or * the bandpass filter settings. */ virtual void OnSelectionChanged(berry::IWorkbenchPart::Pointer source, const QList& nodes) override; /** \brief Instance of the GUI controls */ Ui::PAImageProcessingControls m_Controls; float m_ResampleSpacing; bool m_UseLogfilter; std::string m_OldNodeName; /** \brief Method for updating the BFconfig by using a selected image and the GUI configuration. */ mitk::BeamformingSettings::Pointer CreateBeamformingSettings(mitk::Image::Pointer image); void EnableControls(); void DisableControls(); /** \brief Class through which the filters are called. */ mitk::PhotoacousticFilterService::Pointer m_FilterBank; }; class BeamformingThread : public QThread { Q_OBJECT void run() Q_DECL_OVERRIDE; signals: void result(mitk::Image::Pointer, std::string nameExtension); void updateProgress(int, std::string); void message(std::string); public: BeamformingThread() : m_SignalDelay(0) {} void setConfig(mitk::BeamformingSettings::Pointer BFconfig); void setSignalDelay(float delay); void setInputImage(mitk::Image::Pointer image); void setFilterBank(mitk::PhotoacousticFilterService::Pointer filterBank) { m_FilterBank = filterBank; } protected: mitk::BeamformingSettings::Pointer m_BFconfig; mitk::Image::Pointer m_InputImage; int m_Cutoff; float m_SignalDelay; // [us] mitk::PhotoacousticFilterService::Pointer m_FilterBank; }; class BmodeThread : public QThread { Q_OBJECT void run() Q_DECL_OVERRIDE; signals: void result(mitk::Image::Pointer, std::string nameExtension); public: enum BModeMethod { ShapeDetection, Abs }; void setConfig(bool useLogfilter, double resampleSpacing, mitk::PhotoacousticFilterService::BModeMethod method, bool useGPU); void setInputImage(mitk::Image::Pointer image); void setFilterBank(mitk::PhotoacousticFilterService::Pointer filterBank) { m_FilterBank = filterBank; } protected: mitk::Image::Pointer m_InputImage; mitk::PhotoacousticFilterService::BModeMethod m_Method; bool m_UseLogfilter; double m_ResampleSpacing; bool m_UseGPU; mitk::PhotoacousticFilterService::Pointer m_FilterBank; }; class CropThread : public QThread { Q_OBJECT void run() Q_DECL_OVERRIDE; signals: void result(mitk::Image::Pointer, std::string nameExtension); public: - void setConfig(unsigned int CutAbove, unsigned int CutBelow, unsigned int CutSliceFirst, unsigned int CutSliceLast); + void setConfig(unsigned int CutAbove, unsigned int CutBelow, unsigned int CutRight, unsigned int CutLeft, unsigned int CutSliceFirst, unsigned int CutSliceLast); void setInputImage(mitk::Image::Pointer image); void setFilterBank(mitk::PhotoacousticFilterService::Pointer filterBank) { m_FilterBank = filterBank; } protected: mitk::Image::Pointer m_InputImage; unsigned int m_CutAbove; unsigned int m_CutBelow; + unsigned int m_CutRight; + unsigned int m_CutLeft; unsigned int m_CutSliceLast; unsigned int m_CutSliceFirst; mitk::PhotoacousticFilterService::Pointer m_FilterBank; }; class BandpassThread : public QThread { Q_OBJECT void run() Q_DECL_OVERRIDE; signals: void result(mitk::Image::Pointer, std::string nameExtension); public: void setConfig(float BPHighPass, float BPLowPass, float TukeyAlphaHighPass, float TukeyAlphaLowPass, float TimeSpacing, float SpeedOfSound, bool IsBFImage); void setInputImage(mitk::Image::Pointer image); void setFilterBank(mitk::PhotoacousticFilterService::Pointer filterBank) { m_FilterBank = filterBank; } protected: mitk::Image::Pointer m_InputImage; float m_BPHighPass; float m_BPLowPass; float m_TukeyAlphaHighPass; float m_TukeyAlphaLowPass; float m_TimeSpacing; float m_SpeedOfSound; bool m_IsBFImage; mitk::PhotoacousticFilterService::Pointer m_FilterBank; }; #endif // PAImageProcessing_h diff --git a/Plugins/org.mitk.gui.qt.photoacoustics.imageprocessing/src/internal/PAImageProcessingControls.ui b/Plugins/org.mitk.gui.qt.photoacoustics.imageprocessing/src/internal/PAImageProcessingControls.ui index bba077088b..ab889470f6 100644 --- a/Plugins/org.mitk.gui.qt.photoacoustics.imageprocessing/src/internal/PAImageProcessingControls.ui +++ b/Plugins/org.mitk.gui.qt.photoacoustics.imageprocessing/src/internal/PAImageProcessingControls.ui @@ -1,1021 +1,1078 @@ PAImageProcessingControls 0 0 385 1278 0 0 QmitkTemplate <html><head/><body><p><span style=" font-weight:600;">Batch Processing</span></p></body></html> Start Batch Processing Bandpass true Crop true Save false Save false Save false Beamform true BMode true Save true <html><head/><body><p><span style=" font-weight:600;">B-mode Filter Settings</span></p></body></html> Envelope Detection Envelope Detection Absolute Filter Do Resampling true 0 0 13 0 11 3 0.010000000000000 1.000000000000000 0.010000000000000 0.150000000000000 [mm] Resampled Depth Spacing Logarithmic Compression Use GPU QLabel { color: rgb(255, 0, 0) } <html><head/><body><p align="center"><span style=" font-size:10pt; font-weight:600;">Please select an image!</span></p></body></html> 0 0 Do image processing Apply B-mode Filter Qt::Horizontal <html><head/><body><p><span style=" font-weight:600;">Bandpass Filter Settings</span></p></body></html> QLayout::SetDefaultConstraint 0 0 0 3 0.010000000000000 200.000000000000000 0.100000000000000 1.000000000000000 [MHz] f High Pass [MHz] f Low Pass 0 0 3 200.000000000000000 0.100000000000000 1 200.000000000000000 3000.000000000000000 5.000000000000000 1540.000000000000000 [m/s] Speed of Sound 1.000000000000000 0.100000000000000 0.500000000000000 Tukey Window α High Pass 2 1.000000000000000 0.100000000000000 0.500000000000000 Tukey Window α Low Pass If not checked, treat input as raw US/PA data with Y-axis a time coordinate [us] Assume Spatial Coordinates <html><head/><body><p align="center"><span style=" font-size:10pt; font-weight:600; color:#ff0000;">Please select an image!</span></p></body></html> Apply Bandpass Qt::Horizontal <html><head/><body><p><span style=" font-weight:600;">Crop Filter Settings</span></p></body></html> + + + + 999999999 + + + + + + + Cut Bottom + + + + + + + 999999999 + + + + + + 999999999 + + + 1 + + + 0 + + + + - 99999 + 999999999 0 - + Cut Top - - + + - Cut Bottom + Cut Right - - - - 99999 - - - 5 + + + + Cut Left - - 0 + + + + + + select slices - + false 99999 - - + + - Maximal beamformed slice + minimal beamformed slice - max + min - + false 99999 10 - - + + - minimal beamformed slice - - - min + Maximal beamformed slice - - - - - select slices + max <html><head/><body><p align="center"><span style=" font-size:10pt; font-weight:600; color:#ff0000;">Please select an image!</span></p></body></html> Apply Crop Filer Qt::Horizontal <html><head/><body><p><span style=" font-weight:600;">Beamforming Filter Settings</span></p></body></html> 5 2 - - + + + + <html><head/><body><p>Some setups' hardware produces signal delays that need to be cropped out of the image before performing beamforming. To do this, select this box.</p></body></html> + + + Consider Hardware Delay [µs] + + + false + + + + + + + Apply Beamforming + + + + + 0 0 - - 1 + + + PA Image + + + + + US Image + + + + + + + + + 0 + 0 + - 200.000000000000000 + 64 - 3000.000000000000000 + 1024 - 5.000000000000000 + 128 - 1540.000000000000000 + 128 - - - - Auto Get Depth + + + + + 0 + 0 + - - true + + 256 + + + 16384 + + + 256 + + + 2048 - + - Apply Beamforming + Reconstruction Lines - - + + - + Image Type - - - - Beamforming Method + + + + + 0 + 0 + + + + DAS + + + + + DMAS + + + + + sDMAS + + - - + + - [mm] Scan Depth + Samples - - + + 0 0 - - 3 - - - 0.010000000000000 - - - 9.000000000000000 - - - 0.050000000000000 - - - 0.300000000000000 - + + + Von Hann + + + + + Hamming + + + + + Box + + - - - - Transducer Elements + + + + true - - - - 0 0 - - 4 - - 300.000000000000000 - - - 0.100000000000000 + 100 - 50.000000000000000 + 0 - + [mm] Transducer Pitch - - + + 0 0 64 - 1024 + 2048 128 - 128 + 256 - - + + 0 0 + + 1 + - 256 + 200.000000000000000 - 16384 + 3000.000000000000000 - 256 + 5.000000000000000 - 2048 + 1540.000000000000000 - - - - - 0 - 0 - + + + + [mm] Scan Depth - - 64 + + + + + + Apodization - - 2048 + + + + + + Auto Get Depth - - 128 + + true - - 256 + + + + + + Beamforming Method - + - Samples + Transducer Elements - - + + - Reconstruction Lines + - - - - true - + + 0 0 + + 4 + - 100 + 300.000000000000000 + + + 0.100000000000000 - 0 + 50.000000000000000 - - + + 0 0 - - - DAS - - - - - DMAS - - - - - sDMAS - - - - - - - - - 0 - 0 - + + 3 - - - PA Image - - - - - US Image - - - - - - - - Image Type + + 0.010000000000000 - - - - - - - 0 - 0 - + + 9.000000000000000 - - - Von Hann - - - - - Hamming - - - - - Box - - - - - - - - Apodization + + 0.050000000000000 + + + 0.300000000000000 - - + + - [m/s] Speed of Sound + Probe Geomentry - - - - - - - + - - + + - Compute On GPU - - - true + [m/s] Speed of Sound - + 0 0 1 1.000000000000000 180.000000000000000 27.000000000000000 - + [°] Element Angle - + + + + Compute On GPU + + + true + + + + <html><head/><body><p align="center"><span style=" font-size:10pt; font-weight:600; color:#ff0000;">Please select an image!</span></p></body></html> - + [mm] Reconstruction Depth - + + + + 0.100000000000000 + + + 1.000000000000000 + + + + 4 300.000000000000000 0.100000000000000 40.000000000000000 - - - - 0.100000000000000 - - - 1.000000000000000 - + + + + + Linear + + + + + Concave + + - - - - <html><head/><body><p>Some setups' hardware produces signal delays that need to be cropped out of the image before performing beamforming. To do this, select this box.</p></body></html> + + + + 2 + + + + - Consider Hardware Delay [µs] - - - false + [mm] Radius Qt::Vertical 20 40