diff --git a/Modules/PhotoacousticsLib/MitkSpectralUnmixing/SpectralUnmixingApp.cpp b/Modules/PhotoacousticsLib/MitkSpectralUnmixing/SpectralUnmixingApp.cpp index ab34a53d69..537fc52ed8 100644 --- a/Modules/PhotoacousticsLib/MitkSpectralUnmixing/SpectralUnmixingApp.cpp +++ b/Modules/PhotoacousticsLib/MitkSpectralUnmixing/SpectralUnmixingApp.cpp @@ -1,313 +1,327 @@ /*=================================================================== 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 "mitkPALinearSpectralUnmixingFilter.h" #include "mitkPASpectralUnmixingFilterBase.h" #include "mitkPASpectralUnmixingFilterVigra.h" #include "mitkPASpectralUnmixingSO2.h" #include #include #include #include #include #include "mitkPreferenceListReaderOptionsFunctor.h" struct InputParameters { std::string inputPath; std::string outputPath; - int numberOfInputs; + std::string inputWavelengths; + std::string inputAlg; + std::string inputWeights; + mitkCommandLineParser::StringContainerType test; }; + 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("inputPath", "i", mitkCommandLineParser::InputDirectory, "Input folder (directory)", "input folder", us::Any(), false); parser.addArgument("outputPath", "o", mitkCommandLineParser::OutputDirectory, "Input save folder (directory)", "input save folder", us::Any(), false); - parser.addArgument("numberOfInputs", - "n", - mitkCommandLineParser::Int, - "Number of Input files", - "number of inputs", + parser.addArgument("inputWavelengths", + "l", + mitkCommandLineParser::String, + "Input wavelengths (string)", + "input wavelengths", + us::Any(), + false); + parser.addArgument("inputAlg", + "a", + mitkCommandLineParser::String, + "Input algorithm (string)", + "input algorithm", us::Any(), false); + parser.addArgument("inputWeights", + "w", + mitkCommandLineParser::String, + "Input weights (string)", + "input weights", + us::Any(), + true); + parser.addArgument("test", + "t", + mitkCommandLineParser::StringList, + "Input tes (string)", + "input tes", + 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("inputPath")) { input.inputPath = us::any_cast(parsedArgs["inputPath"]); } else { MITK_ERROR << "Error: No inputPath"; mitkThrow() << "Error: No inputPath"; } if (parsedArgs.count("outputPath")) { input.outputPath = us::any_cast(parsedArgs["outputPath"]); } else { MITK_ERROR << "Error: No outputPath"; mitkThrow() << "Error: No outputPath"; } - if (parsedArgs.count("numberOfInputs")) + + if (parsedArgs.count("inputWavelengths")) { - input.numberOfInputs = us::any_cast(parsedArgs["numberOfInputs"]); + input.inputWavelengths = us::any_cast(parsedArgs["inputWavelengths"]); } else { - MITK_ERROR << "Error: No number of Inputs"; - mitkThrow() << "Error: No number of Inputs"; + 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"]); + } + + if (parsedArgs.count("test")) + { + input.test = us::any_cast(parsedArgs["test"]); + } + MITK_INFO << "Parsing arguments...[Done]"; return input; } -mitk::pa::SpectralUnmixingFilterBase::Pointer GetFilterInstance(std::string algorithm) +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 weigthVec = {40, 45, 47}; + std::vector weigthVec = weights; for (int i = 0; i < 3; ++i) { dynamic_cast(spectralUnmixingFilter.GetPointer()) - ->AddWeight(weigthVec[i]); + ->AddWeight(weigthVec[i]/100); } } return spectralUnmixingFilter; } -std::string numberGenerator(int input) +std::vector inputWavelengthGenerator(std::string input) { - std::string number; - int counter = input + 1; - if (counter < 10) - number = "00" + std::to_string(counter); - else if (counter < 100) - number = "0" + std::to_string(counter); - else if (counter < 1000) - number = std::to_string(counter); - else - number = "NAN"; - return number; -} + std::vector output; -std::string smileyGenerator(float input) -{ - std::string number; - if (input < 10) - number = ":'("; - else if (input < 40) - number = ":("; - else if (input < 60) - number = ":/"; - else if (input < 85) - number = ":)"; - else if (input <= 100) - number = ":D"; - else - number = "NAN"; - return number; + + for (int i = 0; i < input.size(); ++i) + { + if (input[i] == ',') + { + std::string foo; + foo.push_back(input[i-3]); + foo.push_back(input[i-2]); + foo.push_back(input[i-1]); + output.push_back(std::stoi(foo)); + } + } + std::string foo; + foo.push_back(input[input.size() - 3]); + foo.push_back(input[input.size() - 2]); + foo.push_back(input[input.size()-1]); + output.push_back(std::stoi(foo)); + + return output; } int main(int argc, char *argv[]) { auto input = parseInput(argc, argv); - std::string inputDir = input.inputPath; - std::string outputDir = input.outputPath; - unsigned int N = input.numberOfInputs; -/* - //maybee try with "itk system tools" - - //auto test = itksys::SystemTools::GetFilenameName(argv[0]).c_str(); - - //MITK_INFO << "test: " << test; - + auto tst = input.test; - / +++ temporary solution BEGIN +++ - std::vector files; - std::string file; - for (int i = 1; i < 34; ++i) + for (int s = 0; s < tst.size(); ++s) { - - if (i < 10) - { - file = "E:/NHDATA/sdmas_beamformed/merged/static-oxy_sdmas_00" + std::to_string(i) + "_merged.nrrd"; - } - else - { - file = "E:/NHDATA/sdmas_beamformed/merged/static-oxy_sdmas_0" + std::to_string(i) + "_merged.nrrd"; - } - files.push_back(file); + MITK_INFO << " s:" << s << "trs: " << tst[s] << "\n"; + int sdh = std::stoi(tst[s]); + MITK_INFO << "int: " << sdh << "\n"; } - / +++ temporary solution END +++ - std::vector files; - std::string file; - for (int i = 0; i < 7; ++i) - { - file = "E:/NHCAMI/cami-experimental/PAI/spectralUnmixing/inSilico/paImages/selection/noiselevel1_rep1000_wavelength_selction_data_" + - std::to_string(i) + ".nrrd"; - files.push_back(file); - } - std::vector files; - std::string file; - file = "E:/NHCAMI/cami-experimental/PAI/spectralUnmixing/inSilico/paImages/selection/noiselevel1_rep1000_wavelength_selction_data.nrrd"; - files.push_back(file);*/ - - std::vector algorithms = { "QR", "LU", "SVD", "NNLS", "WLS" }; - std::vector files; + std::string inputImage = input.inputPath; + std::string outputDir = input.outputPath; + std::string inputWls = input.inputWavelengths; + std::vector wavelengths = inputWavelengthGenerator(inputWls); + std::string algo = input.inputAlg; - for (unsigned file_idx = 0; file_idx < N; ++file_idx) - { - std::string file = inputDir + "/static-oxy_das_" + numberGenerator(file_idx) + "_merged.nrrd"; - files.push_back(file); - } - - for (unsigned image_idx = 0; image_idx < N; ++image_idx) - { - auto m_inputImage = mitk::IOUtil::Load(files[image_idx]); - for (unsigned alg = 0; alg < 5; ++alg) - { - float percent = (image_idx*1.0 / N + 1.0 / N * (alg / 5.0))*100.0; - MITK_INFO << "IN PROGRESS: " << algorithms[alg] << " " << numberGenerator(image_idx) << " " << percent << "% " << smileyGenerator(percent); - mitk::pa::SpectralUnmixingFilterBase::Pointer m_SpectralUnmixingFilter = GetFilterInstance(algorithms[alg]); - m_SpectralUnmixingFilter->SetInput(m_inputImage); - m_SpectralUnmixingFilter->AddOutputs(2); - 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->AddWavelength(760); - m_SpectralUnmixingFilter->AddWavelength(798); - m_SpectralUnmixingFilter->AddWavelength(858); - - m_SpectralUnmixingFilter->Update(); - - auto output1 = m_SpectralUnmixingFilter->GetOutput(0); - auto output2 = m_SpectralUnmixingFilter->GetOutput(1); - - std::string unmixingOutputHbO2 = outputDir + "/SUOutput/" + algorithms[alg] + "/static_oxy_" + algorithms[alg] + "_HbO2_SU_" + numberGenerator(image_idx) + ".nrrd"; - std::string unmixingOutputHb = outputDir + "/SUOutput/" + algorithms[alg] + "/static_oxy_" + algorithms[alg] + "_Hb_SU_" + numberGenerator(image_idx) + ".nrrd"; - mitk::IOUtil::Save(output1, unmixingOutputHbO2); - mitk::IOUtil::Save(output2, unmixingOutputHb); - - auto m_sO2 = mitk::pa::SpectralUnmixingSO2::New(); - m_sO2->Verbose(false); + + + auto m_inputImage = mitk::IOUtil::Load(inputImage); + + mitk::pa::SpectralUnmixingFilterBase::Pointer m_SpectralUnmixingFilter; + + if (algo == "WLS") + { + std::vector Weights = inputWavelengthGenerator(input.inputWeights); + m_SpectralUnmixingFilter = GetFilterInstance(algo, Weights); + } + else + { + m_SpectralUnmixingFilter = GetFilterInstance(algo); + } + + m_SpectralUnmixingFilter->SetInput(m_inputImage); + + 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 (int wIdx = 0; wIdx < wavelengths.size(); ++wIdx) + { + m_SpectralUnmixingFilter->AddWavelength(wavelengths[wIdx]); + MITK_INFO << wavelengths[wIdx]; + } + + 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_SU_.nrrd"; + std::string unmixingOutputHb = outputDir + "_Hb_SU_.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->SetInput(0, output1); + m_sO2->SetInput(1, output2); - m_sO2->Update(); + m_sO2->Update(); - mitk::Image::Pointer sO2 = m_sO2->GetOutput(0); - sO2->SetSpacing(output1->GetGeometry()->GetSpacing()); + mitk::Image::Pointer sO2 = m_sO2->GetOutput(0); + sO2->SetSpacing(m_inputImage->GetGeometry()->GetSpacing()); - std::string outputSo2 = outputDir + "/sO2/" + algorithms[alg] + "/static_oxy_" + algorithms[alg] + "_sO2_" + numberGenerator(image_idx) + ".nrrd"; - mitk::IOUtil::Save(sO2, outputSo2); + std::string outputSo2 = outputDir + "_sO2_.nrrd"; + mitk::IOUtil::Save(sO2, outputSo2); - m_SpectralUnmixingFilter = nullptr; - m_sO2 = nullptr; - } - } - MITK_INFO << "Spectral Unmixing DONE"; + m_sO2 = nullptr; + m_SpectralUnmixingFilter = nullptr; + + MITK_INFO << "Spectral Unmixing DONE"; } diff --git a/Plugins/org.mitk.gui.qt.photoacoustics.spectralunmixing/src/internal/SpectralUnmixingControls.ui b/Plugins/org.mitk.gui.qt.photoacoustics.spectralunmixing/src/internal/SpectralUnmixingControls.ui index 16095add30..aa90868457 100644 --- a/Plugins/org.mitk.gui.qt.photoacoustics.spectralunmixing/src/internal/SpectralUnmixingControls.ui +++ b/Plugins/org.mitk.gui.qt.photoacoustics.spectralunmixing/src/internal/SpectralUnmixingControls.ui @@ -1,1232 +1,1232 @@ SpectralUnmixingControls 0 0 325 742 0 0 QmitkTemplate QLabel { color: rgb(255, 0, 0) } Please select an image! Qt::Horizontal 75 true Wavelengths settings 0 75 16777215 75 10 <html><head/><body><p>* \brief All values have to be intergers. They need to have the same order than the wavelength at the image.</p></body></html> true Qt::SolidLine 42 50 30 30 λ [nm] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 - 700 + 760 - 710 + 798 - 720 + 858 - 730 + - 740 + - 750 + - 760 + - 770 + - 780 + 790 800 810 820 830 840 850 860 870 880 890 900 910 920 930 940 950 Qt::Horizontal 75 true <html><head/><body><p><span style=" font-weight:400;">* \brief Select as least one Absorber. It's not possible to select more absorbers then wavelengths.</span></p></body></html> Chromophore selection <html><head/><body><p>* \brief One of the main absorbers in near infrared spectrum.</p></body></html> Oxyhemoglobin true true <html><head/><body><p>* \brief One of the main absorbers in near infrared spectrum.</p></body></html> Deoxygenated hemoglobin true true false <html><head/><body><p>* \brief One of the main absorbers in near infrared spectrum.</p></body></html> Melanin <html><head/><body><p>* \brief This endmember will be unmixed with 1 at all wavelgnths.</p></body></html> Static Endmember Qt::Horizontal 75 true <html><head/><body><p><span style=" font-weight:400;">* \brief One needs to choose an spectral unmixing algorithm.</span></p></body></html> Unmixing algorithm true MS Shell Dlg 2 <html><head/><body><p>* \brief For detailed information about the algorithms please have a look at the documentation.</p></body></html> false 21 2147483647 ==CHOSE ALGORITHM== ==QR decomposition== householderQr colPivHouseholderQr fullPivHouseholderQr ==LU decompositon fullPivLu ==Cholesky decompostion== ldlt llt .. ==Least squares== LS jacobiSvd NNLARS NNGoldfarb weighted ==Others== SimplexMax true 0 75 380 81 <html><head/><body><p>* \brief the weights are at the same order as the wavelength</p></body></html> false false Weights [%] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 11 10 9 8 7 7 6 6 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 6 6 6 Qt::Horizontal 75 true Oxygen saturation <html><head/><body><p>* \brief calculates HbO2/(HbO2+Hb) if De- and oxyhemoglobin are selected.</p></body></html> calculate sO2 false 380 82 <html><head/><body><p>* \brief below threshold calculated sO2 value will set to zero</p></body></html> Threshold HbO2 Hb Sum SO2 % Qt::Horizontal 75 true Additional Settings <html><head/><body><p>* \brief This mode will give additional console outputs for debugging.</p></body></html> Verbose Mode (additional console outputs) false <html><head/><body><p>* \brief This checkbox will start Chrono and takes the time between clicking of the &quot;Perform spectral unmixing&quot; button until the GUI enables again.</p></body></html> Chrono <html><head/><body><p>* \brief Calculates the realtive error between unmixing result and the input image in the L2 norm.</p></body></html> Relative error image 307 70 <html><head/><body><p>* \brief below the threshold calculated relative error will set to zero</p></body></html> Threshold HbO2 Hb Do image processing Perform spectral unmixing Qt::Horizontal Qt::Vertical 17 54