diff --git a/Modules/PhotoacousticsAlgorithms/MitkPABeamformingTool/PABeamformingTool.cpp b/Modules/PhotoacousticsAlgorithms/MitkPABeamformingTool/PABeamformingTool.cpp index c79025a0f1..b849b6a95e 100644 --- a/Modules/PhotoacousticsAlgorithms/MitkPABeamformingTool/PABeamformingTool.cpp +++ b/Modules/PhotoacousticsAlgorithms/MitkPABeamformingTool/PABeamformingTool.cpp @@ -1,324 +1,349 @@ /*=================================================================== 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 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 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, 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(); } - for (TiXmlElement* elem = root->FirstChildElement(); elem != NULL; elem = elem->NextSiblingElement()) + + 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"; 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 ); } - if (elemName == "Cropping") + 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; BModeSettings bmodeSettings{ mitk::PhotoacousticFilterService::BModeMethod::EnvelopeDetection, false }; CropSettings cropSettings{ 0,0,0,0,0,0 }; ResampleSettings resSettings{ 0.15 }; ProcessSettings processSettings{ true, false, false }; MITK_INFO << "Parsing settings XML..."; try { ParseXML(input.settingsFile, input, &bfSettings, 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.DoResampling) { 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)}; 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/PhotoacousticsLib/MitkSpectralUnmixing/SpectralUnmixingApp.cpp b/Modules/PhotoacousticsLib/MitkSpectralUnmixing/SpectralUnmixingApp.cpp index 2157b65b76..2a3ad625ea 100644 --- a/Modules/PhotoacousticsLib/MitkSpectralUnmixing/SpectralUnmixingApp.cpp +++ b/Modules/PhotoacousticsLib/MitkSpectralUnmixing/SpectralUnmixingApp.cpp @@ -1,315 +1,315 @@ /*=================================================================== 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 "mitkPALinearSpectralUnmixingFilter.h" #include "mitkPASpectralUnmixingFilterBase.h" #include "mitkPASpectralUnmixingFilterVigra.h" #include "mitkPASpectralUnmixingSO2.h" #include #include #include #include #include #include "mitkPreferenceListReaderOptionsFunctor.h" /* \brief The spectral unmixing mini app (SUMA) is designed to enable batch processing for spectral unmixing. For detailed documentation look into the header files of the included spectral unmixing filters.*/ struct InputParameters { std::string inputFilename; std::string outputFileStruct; // "E:/mydata/awesome_exp_unmixed/a" will be saved as "a_HbO2_SU_.nrrd", "a_Hb_SU_.nrrd" and "a_sO2_.nrrd"; std::string inputAlg; std::string outputFileNumber; mitkCommandLineParser::StringContainerType inputWavelengths; mitkCommandLineParser::StringContainerType inputWeights; }; InputParameters parseInput(int argc, char *argv[]) { //MITK_INFO << "Parsing arguments..."; mitkCommandLineParser parser; parser.setCategory("MITK-Photoacoustics"); parser.setTitle("Mitk Spectral Unmixing App"); parser.setDescription("Batch processing for spectral unmixing."); parser.setContributor("Computer Assisted Medical Interventions, DKFZ"); parser.setArgumentPrefix("--", "-"); parser.beginGroup("Required parameters"); parser.addArgument("inputFilename", "i", - mitkCommandLineParser::InputDirectory, + mitkCommandLineParser::Directory, "Input Filename (NAME.nrrd)", "input filename", us::Any(), false); parser.addArgument("outputFileStruct", "o", - mitkCommandLineParser::OutputDirectory, + mitkCommandLineParser::Directory, "Output save name (name without ending!)", "Output save name", us::Any(), false); parser.addArgument("outputFileNumber", "n", mitkCommandLineParser::String, "Output file number", "Output save number", us::Any(), false); parser.addArgument("inputWavelengths", "l", mitkCommandLineParser::StringList, "Input wavelengths (123 124 125 ... int blank int blank)", "input wavelengths", us::Any(), false); parser.addArgument("inputAlg", "a", mitkCommandLineParser::String, "Input algorithm (string)", "input algorithm", us::Any(), false); parser.addArgument("inputWeights", "w", mitkCommandLineParser::StringList, "Input weights (123 124 125 ... int in % blank int in % blank)", "input weights", us::Any(), true); parser.endGroup(); InputParameters input; std::map parsedArgs = parser.parseArguments(argc, argv); if (argc == 0) exit(-1); //for (int i = 0; i < argc; ++i) //{ // MITK_INFO << argv[i]; //} if (parsedArgs.count("inputFilename")) { input.inputFilename = us::any_cast(parsedArgs["inputFilename"]); } else { MITK_ERROR << "Error: No input file"; mitkThrow() << "Error: No input file"; } if (parsedArgs.count("outputFileStruct")) { input.outputFileStruct = us::any_cast(parsedArgs["outputFileStruct"]); } else { MITK_ERROR << "Error: No output"; mitkThrow() << "Error: No output"; } if (parsedArgs.count("outputFileNumber")) { input.outputFileNumber = us::any_cast(parsedArgs["outputFileNumber"]); } else { MITK_ERROR << "Error: No output number"; mitkThrow() << "Error: No output number"; } if (parsedArgs.count("inputWavelengths")) { input.inputWavelengths = us::any_cast(parsedArgs["inputWavelengths"]); } else { MITK_ERROR << "Error: No wavelengths"; mitkThrow() << "Error: No wavelengths"; } if (parsedArgs.count("inputAlg")) { input.inputAlg = us::any_cast(parsedArgs["inputAlg"]); } else { MITK_ERROR << "Error: No algorithm"; mitkThrow() << "Error: No algorithm"; } if (parsedArgs.count("inputWeights")) { input.inputWeights = us::any_cast(parsedArgs["inputWeights"]); } //MITK_INFO << "Parsing arguments...[Done]"; return input; } // Class takes string and sets algorithm for spectral unmixing in the corresponding filter class mitk::pa::SpectralUnmixingFilterBase::Pointer GetFilterInstance(std::string algorithm, std::vector weights = std::vector()) { mitk::pa::SpectralUnmixingFilterBase::Pointer spectralUnmixingFilter; if (algorithm == "QR") { spectralUnmixingFilter = mitk::pa::LinearSpectralUnmixingFilter::New(); dynamic_cast(spectralUnmixingFilter.GetPointer()) ->SetAlgorithm(mitk::pa::LinearSpectralUnmixingFilter::AlgortihmType::HOUSEHOLDERQR); } else if (algorithm == "SVD") { spectralUnmixingFilter = mitk::pa::LinearSpectralUnmixingFilter::New(); dynamic_cast(spectralUnmixingFilter.GetPointer()) ->SetAlgorithm(mitk::pa::LinearSpectralUnmixingFilter::AlgortihmType::JACOBISVD); } else if (algorithm == "LU") { spectralUnmixingFilter = mitk::pa::LinearSpectralUnmixingFilter::New(); dynamic_cast(spectralUnmixingFilter.GetPointer()) ->SetAlgorithm(mitk::pa::LinearSpectralUnmixingFilter::AlgortihmType::FULLPIVLU); } else if (algorithm == "NNLS") { spectralUnmixingFilter = mitk::pa::SpectralUnmixingFilterVigra::New(); dynamic_cast(spectralUnmixingFilter.GetPointer()) ->SetAlgorithm(mitk::pa::SpectralUnmixingFilterVigra::VigraAlgortihmType::LARS); } else if (algorithm == "WLS") { spectralUnmixingFilter = mitk::pa::SpectralUnmixingFilterVigra::New(); dynamic_cast(spectralUnmixingFilter.GetPointer()) ->SetAlgorithm(mitk::pa::SpectralUnmixingFilterVigra::VigraAlgortihmType::WEIGHTED); std::vector weightVec = weights; for (unsigned int i = 0; i < weightVec.size(); ++i) { dynamic_cast(spectralUnmixingFilter.GetPointer()) ->AddWeight(weightVec[i]); } } return spectralUnmixingFilter; } int main(int argc, char *argv[]) { auto input = parseInput(argc, argv); std::string algo = input.inputAlg; std::string outputDir = input.outputFileStruct; std::string outputNumber = input.outputFileNumber; auto inputWls = input.inputWavelengths; std::vector wavelengths; for (unsigned int s = 0; s < inputWls.size(); ++s) { int wl = std::stoi(inputWls[s]); wavelengths.push_back(wl); //MITK_INFO << "Wavelength: " << wl << "\n"; } mitk::pa::SpectralUnmixingFilterBase::Pointer m_SpectralUnmixingFilter; if (algo == "WLS") { auto inputW = input.inputWeights; std::vector Weights; for (unsigned int s = 0; s < inputW.size(); ++s) { int w = std::stoi(inputW[s]); Weights.push_back(w); //MITK_INFO << "Weights: " << w << "\n"; } m_SpectralUnmixingFilter = GetFilterInstance(algo, Weights); } else { m_SpectralUnmixingFilter = GetFilterInstance(algo); } m_SpectralUnmixingFilter->Verbose(false); m_SpectralUnmixingFilter->RelativeError(false); m_SpectralUnmixingFilter->AddChromophore(mitk::pa::PropertyCalculator::ChromophoreType::OXYGENATED); m_SpectralUnmixingFilter->AddChromophore(mitk::pa::PropertyCalculator::ChromophoreType::DEOXYGENATED); m_SpectralUnmixingFilter->AddOutputs(2); for (unsigned int wIdx = 0; wIdx < wavelengths.size(); ++wIdx) { m_SpectralUnmixingFilter->AddWavelength(wavelengths[wIdx]); //MITK_INFO << wavelengths[wIdx]; } //to add a batch processing: loop for a dir start here; don't forget to set a counter to the three output savenames!!! std::string inputImage = input.inputFilename; auto m_inputImage = mitk::IOUtil::Load(inputImage); m_SpectralUnmixingFilter->SetInput(m_inputImage); m_SpectralUnmixingFilter->Update(); auto output1 = m_SpectralUnmixingFilter->GetOutput(0); auto output2 = m_SpectralUnmixingFilter->GetOutput(1); output1->SetSpacing(m_inputImage->GetGeometry()->GetSpacing()); output2->SetSpacing(m_inputImage->GetGeometry()->GetSpacing()); std::string unmixingOutputHbO2 = outputDir + "HbO2." + outputNumber + ".nrrd"; std::string unmixingOutputHb = outputDir + "Hb." + outputNumber + ".nrrd"; //mitk::IOUtil::Save(output1, unmixingOutputHbO2); //mitk::IOUtil::Save(output2, unmixingOutputHb); auto m_sO2 = mitk::pa::SpectralUnmixingSO2::New(); m_sO2->Verbose(false); m_sO2->SetInput(0, output1); m_sO2->SetInput(1, output2); m_sO2->Update(); mitk::Image::Pointer sO2 = m_sO2->GetOutput(0); mitk::Image::Pointer tHb = m_sO2->GetOutput(1); sO2->SetSpacing(m_inputImage->GetGeometry()->GetSpacing()); tHb->SetSpacing(m_inputImage->GetGeometry()->GetSpacing()); std::string outputSo2 = outputDir + "sO2." + outputNumber + ".nrrd"; mitk::IOUtil::Save(sO2, outputSo2); std::string outputTHb = outputDir + "tHb." + outputNumber + ".nrrd"; mitk::IOUtil::Save(tHb, outputTHb); m_sO2 = nullptr; m_SpectralUnmixingFilter = nullptr; //to add a batch processing: loop for a dir end here MITK_INFO << "Spectral Unmixing DONE"; }