diff --git a/Modules/PhotoacousticsLib/MitkSpectralUnmixing/SpectralUnmixingApp.cpp b/Modules/PhotoacousticsLib/MitkSpectralUnmixing/SpectralUnmixingApp.cpp index c3537c2c01..2157b65b76 100644 --- a/Modules/PhotoacousticsLib/MitkSpectralUnmixing/SpectralUnmixingApp.cpp +++ b/Modules/PhotoacousticsLib/MitkSpectralUnmixing/SpectralUnmixingApp.cpp @@ -1,289 +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..."; + //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, "Input Filename (NAME.nrrd)", "input filename", us::Any(), false); parser.addArgument("outputFileStruct", "o", mitkCommandLineParser::OutputDirectory, - "Input save name (name without ending!)", - "input save name", + "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]; - } + //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]"; + //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_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"; + //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]; + //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.nrrd"; - std::string unmixingOutputHb = outputDir + "_Hb.nrrd"; - mitk::IOUtil::Save(output1, unmixingOutputHbO2); - mitk::IOUtil::Save(output2, unmixingOutputHb); + 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.nrrd"; + 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"; } diff --git a/Modules/PhotoacousticsLib/include/mitkPASpectralUnmixingSO2.h b/Modules/PhotoacousticsLib/include/mitkPASpectralUnmixingSO2.h index c822e102f7..63bd582815 100644 --- a/Modules/PhotoacousticsLib/include/mitkPASpectralUnmixingSO2.h +++ b/Modules/PhotoacousticsLib/include/mitkPASpectralUnmixingSO2.h @@ -1,106 +1,116 @@ /*=================================================================== 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 MITKPHOTOACOUSTICSPECTRALUNMIXINGSO2_H #define MITKPHOTOACOUSTICSPECTRALUNMIXINGSO2_H #include "mitkImageToImageFilter.h" #include //Includes for smart pointer usage #include "mitkCommon.h" #include "itkLightObject.h" namespace mitk { namespace pa { /** * \brief derives out of two identical sized MITK images the oxygen saturation and return one MITK image as result. Furthermore * it is possible to set settings that the result shows just SO2 values above a threshold, or above a input value for Hb, HbO2 to * get just a oxygen saturation image of interessting structures. * * Input: * The input has to be two 3D MITK images. The order of the inputs matter! The first input has to be the Hb image the second input * has to be the HbO2 image. The settings are integer values. The SO2 threshold therefore is percentage value. * * Output: * The output will be one MITK image. Where one can see the oxygen saturation of all pixels above the set threholds. If a pixel is * below a threhold or NAN then the value will be set to zero. + * + * UPDATE: SO2 Filter will get an second output, the total hemoglobin value with ->GetOutput(1). */ class MITKPHOTOACOUSTICSLIB_EXPORT SpectralUnmixingSO2 : public mitk::ImageToImageFilter { public: mitkClassMacro(SpectralUnmixingSO2, mitk::ImageToImageFilter); itkFactorylessNewMacro(Self); /** * \brief AddSO2Settings takes integers and writes them at the end of the m_SO2Settings vector. * @param value of the Setting */ virtual void AddSO2Settings(int value); /** * \brief Verbose gives more information to the console. Default value is false. * @param m_Verbose is the boolian to activate the MITK_INFO logged to the console */ virtual void Verbose(bool verbose); protected: /** * \brief Constructor sets number of input images to two and number of output images to one, respectively. */ SpectralUnmixingSO2(); virtual ~SpectralUnmixingSO2(); std::vector m_SO2Settings; bool m_Verbose = false; private: /** * \brief Inherit from the "ImageToImageFilter" Superclass. Herain it calls InitializeOutputs and the CheckPreConditions * methods and enables pixelwise access to the inputs to calculate the oxygen saturation via the "calculate SO2" method. */ virtual void GenerateData() override; /** * \brief Initialized output images with the same size like the input image. The pixel type is set to float. */ virtual void InitializeOutputs(); /** * \brief Checks if the dimensions of the input images are equal and made out of floating poinhts. * @throws if the inputs don't have the same size * @throws if input images don't consist of floats */ virtual void CheckPreConditions(mitk::Image::Pointer inputHbO2, mitk::Image::Pointer inputHb); /** * \brief calculates HbO2 / (Hb + HbO2) and afterwards checks if the result is significant (SO2ValueNotSiginificant method). * If not the method returns zero otherwise it returns the calculated result. * @param pixelHb is the pixel value of the Hb input. * @param pixelHb is the pixel value of the Hb input. - * @warn if the HbO2 value is NAN (in patricular if Hb == -HbO2), but result will be set to zero + * @warn if the sO2 value is NAN (in patricular if Hb == -HbO2), but result will be set to zero */ float CalculateSO2(float pixelHb, float pixelHbO2); + /** + * \brief calculates (Hb + HbO2). + * @param pixelHb is the pixel value of the Hb input. + * @param pixelHb is the pixel value of the Hb input. + * @warn if the tHb value is NAN (in patricular if Hb == -HbO2), but result will be set to zero + */ + float CalculateTHb(float pixelHb, float pixelHbO2); + /** * \brief return true if SO2 result is not significant by checking if the input values are above the threshold of the settings */ bool SO2ValueNotSiginificant(float Hb, float HbO2, float result); }; } } #endif // MITKPHOTOACOUSTICSPECTRALUNMIXINGSO2_ diff --git a/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingSO2.cpp b/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingSO2.cpp index 0c9bc04bac..85fed87b2b 100644 --- a/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingSO2.cpp +++ b/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingSO2.cpp @@ -1,165 +1,189 @@ /*=================================================================== 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 "mitkPASpectralUnmixingSO2.h" // ImageAccessor #include #include mitk::pa::SpectralUnmixingSO2::SpectralUnmixingSO2() { this->SetNumberOfIndexedInputs(2); - this->SetNumberOfIndexedOutputs(1); + this->SetNumberOfIndexedOutputs(2); this->SetNthOutput(0, mitk::Image::New()); + this->SetNthOutput(1, mitk::Image::New()); + } mitk::pa::SpectralUnmixingSO2::~SpectralUnmixingSO2() { } void mitk::pa::SpectralUnmixingSO2::Verbose(bool verbose) { m_Verbose = verbose; } void mitk::pa::SpectralUnmixingSO2::GenerateData() { MITK_INFO(m_Verbose) << "GENERATING DATA.."; // Get input image mitk::Image::Pointer inputHbO2 = GetInput(0); mitk::Image::Pointer inputHb = GetInput(1); CheckPreConditions(inputHbO2, inputHb); unsigned int xDim = inputHbO2->GetDimensions()[0]; unsigned int yDim = inputHbO2->GetDimensions()[1]; unsigned int zDim = inputHbO2->GetDimensions()[2]; InitializeOutputs(); mitk::ImageReadAccessor readAccessHbO2(inputHbO2); mitk::ImageReadAccessor readAccessHb(inputHb); const float* inputDataArrayHbO2 = ((const float*)readAccessHbO2.GetData()); const float* inputDataArrayHb = ((const float*)readAccessHb.GetData()); auto output = GetOutput(0); + auto output1 = GetOutput(1); + mitk::ImageWriteAccessor writeOutput(output); float* writeBuffer = (float *)writeOutput.GetData(); + mitk::ImageWriteAccessor writeOutput1(output1); + float* writeBuffer1 = (float *)writeOutput1.GetData(); + for (unsigned int x = 0; x < xDim; x++) { for (unsigned int y = 0; y < yDim; y++) { for (unsigned int z = 0;z < zDim; z++) { unsigned int pixelNumber = (xDim*yDim * z) + x * yDim + y; float pixelHb = inputDataArrayHb[pixelNumber]; float pixelHbO2 = inputDataArrayHbO2[pixelNumber]; float resultSO2 = CalculateSO2(pixelHb, pixelHbO2); writeBuffer[(xDim*yDim * z) + x * yDim + y] = resultSO2; + float resultTHb = CalculateTHb(pixelHb, pixelHbO2); + writeBuffer1[(xDim*yDim * z) + x * yDim + y] = resultTHb; } } } MITK_INFO(m_Verbose) << "GENERATING DATA...[DONE]"; } void mitk::pa::SpectralUnmixingSO2::CheckPreConditions(mitk::Image::Pointer inputHbO2, mitk::Image::Pointer inputHb) { unsigned int xDimHb = inputHb->GetDimensions()[0]; unsigned int yDimHb = inputHb->GetDimensions()[1]; unsigned int zDimHb = inputHb->GetDimensions()[2]; unsigned int xDimHbO2 = inputHbO2->GetDimensions()[0]; unsigned int yDimHbO2 = inputHbO2->GetDimensions()[1]; unsigned int zDimHbO2 = inputHbO2->GetDimensions()[2]; if (xDimHb != xDimHbO2 || yDimHb != yDimHbO2 || zDimHb != zDimHbO2) mitkThrow() << "DIMENTIONALITY ERROR!"; if (inputHbO2->GetPixelType() != mitk::MakeScalarPixelType()) mitkThrow() << "PIXELTYPE ERROR! FLOAT REQUIRED"; if (inputHb->GetPixelType() != mitk::MakeScalarPixelType()) mitkThrow() << "PIXELTYPE ERROR! FLOAT REQUIRED"; MITK_INFO(m_Verbose) << "CHECK PRECONDITIONS ...[DONE]"; } void mitk::pa::SpectralUnmixingSO2::InitializeOutputs() { // UNUSED unsigned int numberOfInputs = GetNumberOfIndexedInputs(); unsigned int numberOfOutputs = GetNumberOfIndexedOutputs(); mitk::PixelType pixelType = mitk::MakeScalarPixelType(); const int NUMBER_OF_SPATIAL_DIMENSIONS = 3; auto* dimensions = new unsigned int[NUMBER_OF_SPATIAL_DIMENSIONS]; for(unsigned int dimIdx=0; dimIdxGetDimensions()[dimIdx]; } for (unsigned int outputIdx = 0; outputIdx < numberOfOutputs; outputIdx++) { GetOutput(outputIdx)->Initialize(pixelType, NUMBER_OF_SPATIAL_DIMENSIONS, dimensions); } } float mitk::pa::SpectralUnmixingSO2::CalculateSO2(float Hb, float HbO2) { float result = HbO2 / (Hb + HbO2); if (result != result) { MITK_WARN(m_Verbose) << "SO2 VALUE NAN! WILL BE SET TO ZERO!"; return 0; } else { if (SO2ValueNotSiginificant(Hb, HbO2, result)) { return 0; } else return result; } } +float mitk::pa::SpectralUnmixingSO2::CalculateTHb(float Hb, float HbO2) +{ + float result = (Hb + HbO2); + + if (result != result) + { + MITK_WARN(m_Verbose) << "SO2 VALUE NAN! WILL BE SET TO ZERO!"; + return 0; + } + else + { + return result; + } +} + void mitk::pa::SpectralUnmixingSO2::AddSO2Settings(int value) { m_SO2Settings.push_back(value); } bool mitk::pa::SpectralUnmixingSO2::SO2ValueNotSiginificant(float Hb, float HbO2, float result) { std::vector significant; significant.push_back(HbO2); significant.push_back(Hb); significant.push_back(HbO2 + Hb); significant.push_back(100*(result)); for (unsigned int i = 0; i < m_SO2Settings.size(); ++i) { if (m_SO2Settings[i] != 0 && m_SO2Settings[i] > significant[i] && (std::abs(m_SO2Settings[i] - significant[i]) > 1e-7)) { return true; } } return false; }