diff --git a/Modules/PhotoacousticsLib/include/mitkPASpectralUnmixingFilterBase.h b/Modules/PhotoacousticsLib/include/mitkPASpectralUnmixingFilterBase.h index 076445a0e2..c762eb6da4 100644 --- a/Modules/PhotoacousticsLib/include/mitkPASpectralUnmixingFilterBase.h +++ b/Modules/PhotoacousticsLib/include/mitkPASpectralUnmixingFilterBase.h @@ -1,165 +1,163 @@ /*=================================================================== 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 MITKPHOTOACOUSTICSPECTRALUNMIXINGFILTERBASE_H #define MITKPHOTOACOUSTICSPECTRALUNMIXINGFILTERBASE_H #include "mitkImageToImageFilter.h" #include //Includes for smart pointer usage #include "mitkCommon.h" #include "itkLightObject.h" // Includes for AddEnmemberMatrix #include "mitkPAPropertyCalculator.h" #include namespace mitk { namespace pa { /** * \brief The spectral unmixing filter base is designed as superclass for several spectral unmixing filter eg. Eigen- or Vigrabased ones. * One can add wavelengths and chromophores and get a unmixed output images out of one MITK input image using algorithms from the subclasses. * * Input: * The unmixing input has to be a 3D MITK image where the XY-plane is a image and the Z-direction represents recordings for different * wavelengths. Herein a XY-plane for a specific Z-direction will be called "image". Every image has to be assigned to a certain wavelength. * The "AddWavelength" uses "push_back" to write float values into a vector. The first wavelength will correspond to the first image!!! * If there a more input images 'I' then added wavelengths 'w' the filter base interprets the next x images as repetition sequence of the same * wavelengths. If I % w !=0 the surplus image(s) will be dropped. * Addtionaly one has to add chromophores from the property calculator class enum "ChromophoreType" with the "AddChromophore" method. * This method as well uses "push_back" but the chosen (arbitary) order will be the order of the outputs. * * Output: * The output will be one MITK image per chosen chromophore. Where the XY-plane is a image and the Z-direction represents recordings for different * sequences. * * Subclasses: * - mitkPASpectralUnmixingFilterVigra * - mitkPALinearSpectralUnmixingFilter (uses Eigen algorithms) * - mitkPASpectralUnmixingFilterSimplex * * Possible exceptions: * - "PIXELTYPE ERROR! FLOAT 32 REQUIRED": The MITK input image has to consist out of floats. * - "ERROR! REMOVE WAVELENGTHS!": One needs at least the same amount of images (z-dimension) then added wavelengths. * - "ADD MORE WAVELENGTHS!": One needs at least the same amount of wavelengths then added chromophores. * - "404 CHROMOPHORE NOT FOUND!": The chromophore is not a memeber of the enum "ChromophoreType" OR not supported by the filter. * - "WAVELENGTH XXX nm NOT SUPPORTED!": The wavelength is not part of the proptery calculater data base. The data base can be found @ * [...]\mitk\Modules\PhotoacousticsLib\Resources\spectralLIB.dat */ class MITKPHOTOACOUSTICSLIB_EXPORT SpectralUnmixingFilterBase : public mitk::ImageToImageFilter { public: mitkClassMacro(SpectralUnmixingFilterBase, mitk::ImageToImageFilter); /** * \brief AddChromophore takes mitk::pa::PropertyCalculator::ChromophoreType and writes them at the end of the m_Chromophore vector. * The call of the method sets the order of the GetOutput method! * * @param chromophore has to be element of porperty calculater enum chromophore type * @return for wavelength smaller then 300nm and larger then 1000nm the return will be 0, because not at the data base (2018/06/19) */ void AddChromophore(mitk::pa::PropertyCalculator::ChromophoreType chromophore); /** * \brief AddWavelength takes integers and writes them at the end of the m_Wavelength vector. The first call of the method then * corresponds to the first input image and so on. * @wavelength database supports integers between 300 and 1000 nm */ void AddWavelength(int wavelength); - void AddWeight(unsigned int weight); // has to be a part of the vigra class ofstream myfile; // just for testing purposes; has to be removeed protected: /** * \brief Constructor creats proptery calculater smart pointer new() */ SpectralUnmixingFilterBase(); virtual ~SpectralUnmixingFilterBase(); /** * \brief The subclasses will override the mehtod to calculate the spectral unmixing result vector. * @param endmemberMatrix Matrix with number of chromophores colums and number of wavelengths rows so matrix element (i,j) contains * the absorbtion of chromophore j @ wavelength i taken from the database by PropertyElement method. * @param inputVector Vector containing values of one pixel of XY-plane image with number of wavelength rows (z-dimension of a sequenece) * so the pixelvalue of the first wavelength is stored in inputVector[0] and so on. */ virtual Eigen::VectorXf SpectralUnmixingAlgorithm(Eigen::Matrix endmemberMatrix, Eigen::VectorXf inputVector) = 0; /* * \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); bool m_Verbose = false; std::vector m_Chromophore; std::vector m_Wavelength; - std::vector m_Weight; // has to be a part of the vigra class private: /* * \brief Initialized output images with the same XY-plane size like the input image and total size in z-direction equals number of sequences. * The pixel type is set to float. * @param totalNumberOfSequences = (unsigned int) (numberOfInputImages / numberOfWavelength) >= 1 */ virtual void InitializeOutputs(unsigned int totalNumberOfSequences); /* * \brief Checks if there are a suitable amount of wavelengths and if the input image consists of floats. * @param numberOfInputImages corresponds to the z-dimension of the MITK input image * @param inputDataArray contains every pixel of the MITK input image. Accessable like pseudo multidimensional array, that means * pixel(x,y,z) = (xdim * yDim * z) + x * yDim + y where eg. "xdim" means the total size of the MITK input image in x direction. * @throw if there are more wavelength then images * @throw if there are more chromophores then wavelengths * @throw if the pixel type is not float 32 */ virtual void CheckPreConditions(unsigned int numberOfInputImages, const float* inputDataArray); /* * \brief Inherit from the "ImageToImageFilter" Superclass. Herain it calls InitializeOutputs, CalculateEndmemberMatrix and * CheckPreConditions methods and enables pixelwise access to do spectral unmixing with the "SpectralUnmixingAlgorithm" * method. In the end the method writes the results into the new MITK output images. */ virtual void GenerateData() override; /* * \brief Creats a Matrix with number of chromophores colums and number of wavelengths rows so matrix element (i,j) contains * the absorbtion of chromophore j @ wavelength i. The absorbtion values are taken from the "PropertyElement" method. * @param m_Chromophore is a vector of "PropertyCalculator::ChromophoreType" containing all selected chromophores for the unmixing * @param m_Wavelength is a vector of integers containing all wavelengths of one sequence */ virtual Eigen::Matrix CalculateEndmemberMatrix( std::vector m_Chromophore, std::vector m_Wavelength); /* * \brief "PropertyElement" is the tool to access the absorbtion values out of the database using mitk::pa::PropertyCalculator::GetAbsorptionForWavelengt * and checks if the requested wavelength is part of the database (not zero values). The "ONEENDMEMBER" is a pseudo absorber with static absorbtion 1 * at every wavelength and is therefor not part of the database. If this one is the selected chromophore the return value is 1 for every wavelength. * @param wavelength has to be integer between 300 and 1000 nm * @param chromophore has to be a mitk::pa::PropertyCalculator::ChromophoreType * @throw if mitk::pa::PropertyCalculator::GetAbsorptionForWavelengt returns 0, because this means that the delivered wavelength is not in the database. */ virtual float PropertyElement(mitk::pa::PropertyCalculator::ChromophoreType, int wavelength); PropertyCalculator::Pointer m_PropertyCalculatorEigen; }; } } #endif // MITKPHOTOACOUSTICSPECTRALUNMIXINGFILTERBASE_ diff --git a/Modules/PhotoacousticsLib/include/mitkPASpectralUnmixingFilterVigra.h b/Modules/PhotoacousticsLib/include/mitkPASpectralUnmixingFilterVigra.h index 7ecbd0c27a..02fa4a9c53 100644 --- a/Modules/PhotoacousticsLib/include/mitkPASpectralUnmixingFilterVigra.h +++ b/Modules/PhotoacousticsLib/include/mitkPASpectralUnmixingFilterVigra.h @@ -1,56 +1,59 @@ /*=================================================================== 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 MITKPHOTOACOUSTICSPECTRALUNMIXINGFILTERVIGRA_H #define MITKPHOTOACOUSTICSPECTRALUNMIXINGFILTERVIGRA_H #include #include #include namespace mitk { namespace pa { class MITKPHOTOACOUSTICSLIB_EXPORT SpectralUnmixingFilterVigra : public SpectralUnmixingFilterBase { public: mitkClassMacro(SpectralUnmixingFilterVigra, SpectralUnmixingFilterBase) itkFactorylessNewMacro(Self) enum VigraAlgortihmType { LARS, GOLDFARB, WEIGHTED, LS, vigratest }; + void AddWeight(unsigned int weight); + void mitk::pa::SpectralUnmixingFilterVigra::SetAlgorithm(VigraAlgortihmType SetAlgorithmIndex); protected: SpectralUnmixingFilterVigra(); virtual ~SpectralUnmixingFilterVigra(); + std::vector m_Weight; virtual Eigen::VectorXf SpectralUnmixingAlgorithm(Eigen::Matrix EndmemberMatrix, Eigen::VectorXf inputVector) override; private: mitk::pa::SpectralUnmixingFilterVigra::VigraAlgortihmType algorithmIndex; }; } } #endif // MITKPHOTOACOUSTICSPECTRALUNMIXINGFILTERVIGRA_H diff --git a/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingFilterBase.cpp b/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingFilterBase.cpp index e1e3e1d576..2e30222e60 100644 --- a/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingFilterBase.cpp +++ b/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingFilterBase.cpp @@ -1,206 +1,201 @@ /*=================================================================== 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 "mitkPASpectralUnmixingFilterBase.h" // Includes for AddEnmemberMatrix #include "mitkPAPropertyCalculator.h" #include // ImageAccessor #include #include #include mitk::pa::SpectralUnmixingFilterBase::SpectralUnmixingFilterBase() { this->SetNumberOfIndexedOutputs(3);// find solution --> 4 is max outputs for (unsigned int i = 0; iSetNthOutput(i, mitk::Image::New()); } m_PropertyCalculatorEigen = mitk::pa::PropertyCalculator::New(); } mitk::pa::SpectralUnmixingFilterBase::~SpectralUnmixingFilterBase() { } void mitk::pa::SpectralUnmixingFilterBase::AddWavelength(int wavelength) { m_Wavelength.push_back(wavelength); } void mitk::pa::SpectralUnmixingFilterBase::AddChromophore(mitk::pa::PropertyCalculator::ChromophoreType chromophore) { m_Chromophore.push_back(chromophore); } -void mitk::pa::SpectralUnmixingFilterBase::AddWeight(unsigned int weight) -{ - m_Weight.push_back(weight); -} - void mitk::pa::SpectralUnmixingFilterBase::Verbose(bool verbose) { m_Verbose = verbose; } void mitk::pa::SpectralUnmixingFilterBase::GenerateData() { MITK_INFO(m_Verbose) << "GENERATING DATA.."; mitk::Image::Pointer input = GetInput(0); unsigned int xDim = input->GetDimensions()[0]; unsigned int yDim = input->GetDimensions()[1]; unsigned int numberOfInputImages = input->GetDimensions()[2]; MITK_INFO(m_Verbose) << "x dimension: " << xDim; MITK_INFO(m_Verbose) << "y dimension: " << yDim; MITK_INFO(m_Verbose) << "z dimension: " << numberOfInputImages; unsigned int sequenceSize = m_Wavelength.size(); unsigned int totalNumberOfSequences = numberOfInputImages / sequenceSize; MITK_INFO(m_Verbose) << "TotalNumberOfSequences: " << totalNumberOfSequences; InitializeOutputs(totalNumberOfSequences); auto endmemberMatrix = CalculateEndmemberMatrix(m_Chromophore, m_Wavelength); // Copy input image into array mitk::ImageReadAccessor readAccess(input); const float* inputDataArray = ((const float*)readAccess.GetData()); CheckPreConditions(numberOfInputImages, inputDataArray); // test to see pixel values @ txt file myfile.open("SimplexNormalisation.txt"); for (unsigned int sequenceCounter = 0; sequenceCounter < totalNumberOfSequences; ++sequenceCounter) { MITK_INFO(m_Verbose) << "SequenceCounter: " << sequenceCounter; //loop over every pixel in XY-plane for (unsigned int x = 0; x < xDim; x++) { for (unsigned int y = 0; y < yDim; y++) { Eigen::VectorXf inputVector(sequenceSize); for (unsigned int z = 0; z < sequenceSize; z++) { /** * 'sequenceCounter*sequenceSize' has to be added to 'z' to ensure that one accesses the * correct pixel, because the inputDataArray contains the information of all sequences and * not just the one of the current sequence. */ unsigned int pixelNumber = (xDim*yDim*(z+sequenceCounter*sequenceSize)) + x * yDim + y; auto pixel = inputDataArray[pixelNumber]; //write all wavelength absorbtion values for one(!) pixel of a sequence to a vector inputVector[z] = pixel; } Eigen::VectorXf resultVector = SpectralUnmixingAlgorithm(endmemberMatrix, inputVector); for (int outputIdx = 0; outputIdx < GetNumberOfIndexedOutputs(); ++outputIdx) { auto output = GetOutput(outputIdx); mitk::ImageWriteAccessor writeOutput(output); float* writeBuffer = (float *)writeOutput.GetData(); writeBuffer[(xDim*yDim * sequenceCounter) + x * yDim + y] = resultVector[outputIdx]; } } } } std::chrono::steady_clock::time_point _end(std::chrono::steady_clock::now()); std::chrono::steady_clock::time_point _test(std::chrono::steady_clock::now()); MITK_INFO(m_Verbose) << "GENERATING DATA...[DONE]"; myfile.close(); } void mitk::pa::SpectralUnmixingFilterBase::CheckPreConditions(unsigned int numberOfInputImages, const float* inputDataArray) { MITK_INFO(m_Verbose) << "CHECK PRECONDITIONS ..."; if (m_Wavelength.size() < numberOfInputImages) MITK_WARN << "NUMBER OF WAVELENGTHS < NUMBER OF INPUT IMAGES"; if (m_Wavelength.size() > numberOfInputImages) mitkThrow() << "ERROR! REMOVE WAVELENGTHS!"; if (m_Chromophore.size() > m_Wavelength.size()) mitkThrow() << "ADD MORE WAVELENGTHS!"; if (typeid(inputDataArray[0]).name() != typeid(float).name()) mitkThrow() << "PIXELTYPE ERROR! FLOAT 32 REQUIRED"; MITK_INFO(m_Verbose) << "...[DONE]"; } void mitk::pa::SpectralUnmixingFilterBase::InitializeOutputs(unsigned int totalNumberOfSequences) { MITK_INFO(m_Verbose) << "Initialize Outputs ..."; unsigned int numberOfInputs = GetNumberOfIndexedInputs(); unsigned int numberOfOutputs = GetNumberOfIndexedOutputs(); MITK_INFO(m_Verbose) << "Inputs: " << numberOfInputs << " Outputs: " << numberOfOutputs; 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; dimIdx < 2; dimIdx++) { dimensions[dimIdx] = GetInput()->GetDimensions()[dimIdx]; } dimensions[2] = totalNumberOfSequences; for (unsigned int outputIdx = 0; outputIdx < numberOfOutputs; outputIdx++) GetOutput(outputIdx)->Initialize(pixelType, NUMBER_OF_SPATIAL_DIMENSIONS, dimensions); MITK_INFO(m_Verbose) << "...[DONE]"; } Eigen::Matrix mitk::pa::SpectralUnmixingFilterBase::CalculateEndmemberMatrix( std::vector m_Chromophore, std::vector m_Wavelength) { unsigned int numberOfChromophores = m_Chromophore.size(); //columns unsigned int numberOfWavelengths = m_Wavelength.size(); //rows Eigen::Matrix endmemberMatrixEigen(numberOfWavelengths, numberOfChromophores); for (unsigned int j = 0; j < numberOfChromophores; ++j) { for (unsigned int i = 0; i < numberOfWavelengths; ++i) endmemberMatrixEigen(i, j) = PropertyElement(m_Chromophore[j], m_Wavelength[i]); } MITK_INFO(m_Verbose) << "GENERATING ENMEMBERMATRIX [DONE]"; return endmemberMatrixEigen; } float mitk::pa::SpectralUnmixingFilterBase::PropertyElement(mitk::pa::PropertyCalculator::ChromophoreType chromophore, int wavelength) { if (chromophore == mitk::pa::PropertyCalculator::ChromophoreType::ONEENDMEMBER) return 1; else { float value = m_PropertyCalculatorEigen->GetAbsorptionForWavelength(chromophore, wavelength); if (value == 0) mitkThrow() << "WAVELENGTH " << wavelength << "nm NOT SUPPORTED!"; else return value; } } diff --git a/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingFilterVigra.cpp b/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingFilterVigra.cpp index e4f9fd94d7..48975b48d2 100644 --- a/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingFilterVigra.cpp +++ b/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingFilterVigra.cpp @@ -1,126 +1,130 @@ /*=================================================================== 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 "mitkPASpectralUnmixingFilterVigra.h" // Includes for AddEnmemberMatrix #include "mitkPAPropertyCalculator.h" #include // Testing algorithms #include // ImageAccessor #include #include //vigra #include #include #include mitk::pa::SpectralUnmixingFilterVigra::SpectralUnmixingFilterVigra() { } mitk::pa::SpectralUnmixingFilterVigra::~SpectralUnmixingFilterVigra() { } void mitk::pa::SpectralUnmixingFilterVigra::SetAlgorithm(mitk::pa::SpectralUnmixingFilterVigra::VigraAlgortihmType SetAlgorithmIndex) { algorithmIndex = SetAlgorithmIndex; } +void mitk::pa::SpectralUnmixingFilterVigra::AddWeight(unsigned int weight) +{ + m_Weight.push_back(weight); +} Eigen::VectorXf mitk::pa::SpectralUnmixingFilterVigra::SpectralUnmixingAlgorithm( Eigen::Matrix EndmemberMatrix, Eigen::VectorXf inputVector) { int numberOfWavelengths = EndmemberMatrix.rows(); int numberOfChromophores = EndmemberMatrix.cols(); Eigen::VectorXf resultVector(numberOfChromophores); using namespace vigra; using namespace vigra::linalg; std::vector A_data; std::vector B_data; for (int i = 0; i < numberOfWavelengths; ++i) { B_data.push_back((double)inputVector(i)); for (int j = 0; j < numberOfChromophores; ++j) { A_data.push_back((double)EndmemberMatrix(i, j)); } } const double* Adat = A_data.data(); const double* Bdat = B_data.data(); vigra::Matrix A(Shape2(numberOfWavelengths, numberOfChromophores), Adat); vigra::Matrix b(Shape2(numberOfWavelengths, 1), Bdat); vigra::Matrix x(Shape2(numberOfChromophores, 1)); // minimize (A*x-b)^2 s.t. x>=0 using least angle regression (LARS algorithm) if (mitk::pa::SpectralUnmixingFilterVigra::VigraAlgortihmType::LARS == algorithmIndex) nonnegativeLeastSquares(A, b, x); else if (mitk::pa::SpectralUnmixingFilterVigra::VigraAlgortihmType::GOLDFARB == algorithmIndex) { vigra::Matrix eye(identityMatrix(numberOfChromophores)), zeros(Shape2(numberOfChromophores, 1)), empty, U = transpose(A)*A, v = -transpose(A)*b; x = 0; // minimize (A*x-b)^2 s.t. x>=0 using the Goldfarb-Idnani algorithm quadraticProgramming(U, v, empty, empty, eye, zeros, x); } else if (mitk::pa::SpectralUnmixingFilterVigra::VigraAlgortihmType::WEIGHTED == algorithmIndex) { std::vector weightsvec; for (int i = 0; i < m_Weight.size(); ++i) { double value = double (m_Weight[i]) / 100.0; weightsvec.push_back(value); } const double* weightsdat = weightsvec.data(); vigra::Matrix weigths(Shape2(numberOfWavelengths, 1), weightsdat); vigra::linalg::weightedLeastSquares(A, b, weigths, x); } else if (mitk::pa::SpectralUnmixingFilterVigra::VigraAlgortihmType::LS == algorithmIndex) linearSolve(A, b, x); else if (mitk::pa::SpectralUnmixingFilterVigra::VigraAlgortihmType::vigratest == algorithmIndex) mitkThrow() << "nothing implemented"; else mitkThrow() << "404 VIGRA ALGORITHM NOT FOUND"; for (int k = 0; k < numberOfChromophores; ++k) resultVector[k] = (float)x(k, 0); return resultVector; } diff --git a/Plugins/org.mitk.gui.qt.photoacoustics.spectralunmixing/src/internal/SpectralUnmixing.cpp b/Plugins/org.mitk.gui.qt.photoacoustics.spectralunmixing/src/internal/SpectralUnmixing.cpp index 4e9e4be3da..a4ef67f8dc 100644 --- a/Plugins/org.mitk.gui.qt.photoacoustics.spectralunmixing/src/internal/SpectralUnmixing.cpp +++ b/Plugins/org.mitk.gui.qt.photoacoustics.spectralunmixing/src/internal/SpectralUnmixing.cpp @@ -1,384 +1,385 @@ /*=================================================================== 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 "SpectralUnmixing.h" // Qt #include // mitk image #include // Include to perform Spectral Unmixing #include "mitkPASpectralUnmixingFilterBase.h" #include "mitkPALinearSpectralUnmixingFilter.h" #include "mitkPASpectralUnmixingSO2.h" #include "mitkPASpectralUnmixingFilterVigra.h" #include "mitkPASpectralUnmixingFilterLagrange.h" #include "mitkPASpectralUnmixingFilterSimplex.h" const std::string SpectralUnmixing::VIEW_ID = "org.mitk.views.spectralunmixing"; void SpectralUnmixing::SetFocus() { m_Controls.buttonPerformImageProcessing->setFocus(); } void SpectralUnmixing::CreateQtPartControl(QWidget *parent) { // create GUI widgets from the Qt Designer's .ui file m_Controls.setupUi(parent); connect(m_Controls.buttonPerformImageProcessing, &QPushButton::clicked, this, &SpectralUnmixing::DoImageProcessing); } void SpectralUnmixing::ClearWavelength() { m_Wavelengths.clear(); } void SpectralUnmixing::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.buttonPerformImageProcessing->setEnabled(true); return; } } m_Controls.labelWarning->setVisible(true); m_Controls.buttonPerformImageProcessing->setEnabled(false); } mitk::pa::SpectralUnmixingFilterBase::Pointer SpectralUnmixing::GetFilterInstance(std::string algorithm) { mitk::pa::SpectralUnmixingFilterBase::Pointer spectralUnmixingFilter; if (algorithm == "householderQr") { spectralUnmixingFilter = mitk::pa::LinearSpectralUnmixingFilter::New(); dynamic_cast(spectralUnmixingFilter.GetPointer()) ->SetAlgorithm(mitk::pa::LinearSpectralUnmixingFilter::AlgortihmType::householderQr); } else if (algorithm == "ldlt") { spectralUnmixingFilter = mitk::pa::LinearSpectralUnmixingFilter::New(); dynamic_cast(spectralUnmixingFilter.GetPointer()) ->SetAlgorithm(mitk::pa::LinearSpectralUnmixingFilter::AlgortihmType::ldlt); } else if (algorithm == "llt") { spectralUnmixingFilter = mitk::pa::LinearSpectralUnmixingFilter::New(); dynamic_cast(spectralUnmixingFilter.GetPointer()) ->SetAlgorithm(mitk::pa::LinearSpectralUnmixingFilter::AlgortihmType::llt); } else if (algorithm == "colPivHouseholderQr") { spectralUnmixingFilter = mitk::pa::LinearSpectralUnmixingFilter::New(); dynamic_cast(spectralUnmixingFilter.GetPointer()) ->SetAlgorithm(mitk::pa::LinearSpectralUnmixingFilter::AlgortihmType::colPivHouseholderQr); } else if (algorithm == "jacobiSvd") { spectralUnmixingFilter = mitk::pa::LinearSpectralUnmixingFilter::New(); dynamic_cast(spectralUnmixingFilter.GetPointer()) ->SetAlgorithm(mitk::pa::LinearSpectralUnmixingFilter::AlgortihmType::jacobiSvd); } else if (algorithm == "fullPivLu") { spectralUnmixingFilter = mitk::pa::LinearSpectralUnmixingFilter::New(); dynamic_cast(spectralUnmixingFilter.GetPointer()) ->SetAlgorithm(mitk::pa::LinearSpectralUnmixingFilter::AlgortihmType::fullPivLu); } else if (algorithm == "fullPivHouseholderQr") { spectralUnmixingFilter = mitk::pa::LinearSpectralUnmixingFilter::New(); dynamic_cast(spectralUnmixingFilter.GetPointer()) ->SetAlgorithm(mitk::pa::LinearSpectralUnmixingFilter::AlgortihmType::fullPivHouseholderQr); } else if (algorithm == "test") { spectralUnmixingFilter = mitk::pa::LinearSpectralUnmixingFilter::New(); dynamic_cast(spectralUnmixingFilter.GetPointer()) ->SetAlgorithm(mitk::pa::LinearSpectralUnmixingFilter::AlgortihmType::test); } else if (algorithm == "NNLARS") { spectralUnmixingFilter = mitk::pa::SpectralUnmixingFilterVigra::New(); dynamic_cast(spectralUnmixingFilter.GetPointer()) ->SetAlgorithm(mitk::pa::SpectralUnmixingFilterVigra::VigraAlgortihmType::LARS); } else if (algorithm == "NNGoldfarb") { spectralUnmixingFilter = mitk::pa::SpectralUnmixingFilterVigra::New(); dynamic_cast(spectralUnmixingFilter.GetPointer()) ->SetAlgorithm(mitk::pa::SpectralUnmixingFilterVigra::VigraAlgortihmType::GOLDFARB); } else if (algorithm == "weighted") { spectralUnmixingFilter = mitk::pa::SpectralUnmixingFilterVigra::New(); dynamic_cast(spectralUnmixingFilter.GetPointer()) ->SetAlgorithm(mitk::pa::SpectralUnmixingFilterVigra::VigraAlgortihmType::WEIGHTED); //Tranfer GUI information(Weights) to filter int colunm = 0; int Weight = 1; while (m_Controls.inputtable->item(1, colunm) && Weight > 0) { QString Text = m_Controls.inputtable->item(1, colunm)->text(); Weight = Text.toInt(); if (Weight > 0) - spectralUnmixingFilter->AddWeight(Weight); + dynamic_cast(spectralUnmixingFilter.GetPointer()) + ->AddWeight(Weight); ++colunm; } } else if (algorithm == "LS") { spectralUnmixingFilter = mitk::pa::SpectralUnmixingFilterVigra::New(); dynamic_cast(spectralUnmixingFilter.GetPointer()) ->SetAlgorithm(mitk::pa::SpectralUnmixingFilterVigra::VigraAlgortihmType::LS); } else if (algorithm == "vigratest") { spectralUnmixingFilter = mitk::pa::SpectralUnmixingFilterVigra::New(); dynamic_cast(spectralUnmixingFilter.GetPointer()) ->SetAlgorithm(mitk::pa::SpectralUnmixingFilterVigra::VigraAlgortihmType::vigratest); } else if (algorithm == "SimplexMax") { spectralUnmixingFilter = mitk::pa::SpectralUnmixingFilterSimplex::New(); } else mitkThrow() << "404 ALGORITHM NOT FOUND!"; return spectralUnmixingFilter; } void SpectralUnmixing::DoImageProcessing() { QList nodes = this->GetDataManagerSelection(); if (nodes.empty()) return; mitk::DataNode *node = nodes.front(); if (!node) { // Nothing selected. Inform the user and return QMessageBox::information(nullptr, "Template", "Please load and select an image before starting image processing."); return; } // here we have a valid mitk::DataNode // a node itself is not very useful, we need its data item (the image) 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 SPECTRAL UNMIXING "; if (node->GetName(name)) { // a property called "name" was found for this DataNode message << "'" << name << "'"; } message << "."; MITK_INFO << message.str(); //Tranfer GUI information(Algortihm) to filter auto qs = m_Controls.QComboBoxAlgorithm->currentText(); std::string Algorithm = qs.toUtf8().constData(); mitk::pa::SpectralUnmixingFilterBase::Pointer m_SpectralUnmixingFilter = GetFilterInstance(Algorithm); m_SpectralUnmixingFilter->SetInput(image); //Tranfer GUI information(Wavelength) to filter ClearWavelength(); int col = 0; int Wavelength = 1; while (m_Controls.inputtable->item(0, col) && Wavelength > 0) { QString Text = m_Controls.inputtable->item(0, col)->text(); Wavelength = Text.toInt(); if (Wavelength > 0) m_SpectralUnmixingFilter->AddWavelength(Wavelength); ++col; } //Tranfer GUI information(Chromophores) to filter unsigned int numberofChromophores = 0; DeOxbool = m_Controls.checkBoxDeOx->isChecked(); Oxbool = m_Controls.checkBoxOx->isChecked(); bool Melaninbool = m_Controls.checkBoxMelanin->isChecked(); bool Onebool = m_Controls.checkBoxAdd->isChecked(); if (DeOxbool || Oxbool) { MITK_INFO << "CHOSEN CHROMOPHORES:"; } if (Oxbool) { numberofChromophores += 1; MITK_INFO << "- Oxyhemoglobin"; m_SpectralUnmixingFilter->AddChromophore( mitk::pa::PropertyCalculator::ChromophoreType::OXYGENATED); } if (DeOxbool) { numberofChromophores += 1; MITK_INFO << "- Deoxygenated hemoglobin"; m_SpectralUnmixingFilter->AddChromophore( mitk::pa::PropertyCalculator::ChromophoreType::DEOXYGENATED); } if (Melaninbool) { numberofChromophores += 1; MITK_INFO << "- Melanin"; m_SpectralUnmixingFilter->AddChromophore( mitk::pa::PropertyCalculator::ChromophoreType::MELANIN); } if (Onebool) { numberofChromophores += 1; MITK_INFO << "- Additional Chromophore"; m_SpectralUnmixingFilter->AddChromophore( mitk::pa::PropertyCalculator::ChromophoreType::ONEENDMEMBER); } if (numberofChromophores == 0) { mitkThrow() << "PRESS 'IGNORE' AND CHOOSE A CHROMOPHORE!"; } MITK_INFO << "Updating Filter..."; m_SpectralUnmixingFilter->Update(); int outputCounter = 0; // Write Output images to Data Storage if (Oxbool) { mitk::Image::Pointer HbO2 = m_SpectralUnmixingFilter->GetOutput(outputCounter++); mitk::DataNode::Pointer dataNodeHbO2 = mitk::DataNode::New(); dataNodeHbO2->SetData(HbO2); dataNodeHbO2->SetName("HbO2 " + Algorithm); this->GetDataStorage()->Add(dataNodeHbO2); } if (DeOxbool) { mitk::Image::Pointer Hb = m_SpectralUnmixingFilter->GetOutput(outputCounter++); mitk::DataNode::Pointer dataNodeHb = mitk::DataNode::New(); dataNodeHb->SetData(Hb); dataNodeHb->SetName("Hb " + Algorithm); this->GetDataStorage()->Add(dataNodeHb); } if (Melaninbool) { mitk::Image::Pointer Melanin = m_SpectralUnmixingFilter->GetOutput(outputCounter++); mitk::DataNode::Pointer dataNodeMelanin = mitk::DataNode::New(); dataNodeMelanin->SetData(Melanin); dataNodeMelanin->SetName("Melanin " + Algorithm); this->GetDataStorage()->Add(dataNodeMelanin); } if (Onebool) { mitk::Image::Pointer One = m_SpectralUnmixingFilter->GetOutput(outputCounter++); mitk::DataNode::Pointer dataNodeOne = mitk::DataNode::New(); dataNodeOne->SetData(One); dataNodeOne->SetName("One " + Algorithm); this->GetDataStorage()->Add(dataNodeOne); } //Calculate oxygen saturation bool sO2bool = m_Controls.checkBoxsO2->isChecked(); if (sO2bool) { if (!DeOxbool) mitkThrow() << "SELECT CHROMOPHORE DEOXYHEMOGLOBIN!"; if (!Oxbool) mitkThrow() << "SELECT CHROMOPHORE OXYHEMOGLOBIN!"; MITK_INFO << "CALCULATE OXYGEN SATURATION ..."; auto m_sO2 = mitk::pa::SpectralUnmixingSO2::New(); // Oxygen Saturation Setting for (int i = 0; i < 4; ++i) { if (m_Controls.inputtable->item(0, i)) { QString Text = m_Controls.tableSO2->item(0, i)->text(); float value = Text.toFloat(); MITK_INFO << "value: " << value; m_sO2->AddSO2Settings(value); } else m_sO2->AddSO2Settings(0); } // Initialize pipeline from SU filter class to SO2 class auto output1 = m_SpectralUnmixingFilter->GetOutput(0); auto output2 = m_SpectralUnmixingFilter->GetOutput(1); m_sO2->SetInput(0, output1); m_sO2->SetInput(1, output2); m_sO2->Update(); // Write Output images to Data Storage mitk::Image::Pointer sO2 = m_sO2->GetOutput(0); mitk::DataNode::Pointer dataNodesO2 = mitk::DataNode::New(); dataNodesO2->SetData(sO2); dataNodesO2->SetName("sO2" + Algorithm); this->GetDataStorage()->Add(dataNodesO2); MITK_INFO << "[DONE]"; } mitk::RenderingManager::GetInstance()->InitializeViewsByBoundingObjects(this->GetDataStorage()); MITK_INFO << "Adding images to DataStorage...[DONE]"; } } }