diff --git a/Modules/PhotoacousticsLib/include/mitkPALinearSpectralUnmixingFilter.h b/Modules/PhotoacousticsLib/include/mitkPALinearSpectralUnmixingFilter.h index e120f71c91..8f3dc5cbad 100644 --- a/Modules/PhotoacousticsLib/include/mitkPALinearSpectralUnmixingFilter.h +++ b/Modules/PhotoacousticsLib/include/mitkPALinearSpectralUnmixingFilter.h @@ -1,66 +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 MITKLINEARPHOTOACOUSTICSPECTRALUNMIXINGFILTER_H #define MITKLINEARPHOTOACOUSTICSPECTRALUNMIXINGFILTER_H #include #include namespace mitk { namespace pa { class MITKPHOTOACOUSTICSLIB_EXPORT LinearSpectralUnmixingFilter : public SpectralUnmixingFilterBase { public: mitkClassMacro(LinearSpectralUnmixingFilter, SpectralUnmixingFilterBase) itkFactorylessNewMacro(Self) enum AlgortihmType { householderQr, ldlt, llt, colPivHouseholderQr, jacobiSvd, fullPivLu, fullPivHouseholderQr, test }; void mitk::pa::LinearSpectralUnmixingFilter::SetAlgorithm(AlgortihmType SetAlgorithmIndex); protected: LinearSpectralUnmixingFilter(); virtual ~LinearSpectralUnmixingFilter(); - - // Test algorithm for SU --> later a new class should be set up virtual Eigen::VectorXf SpectralUnmixingAlgorithm(Eigen::Matrix EndmemberMatrix, Eigen::VectorXf inputVector) override; - virtual Eigen::VectorXf NNLSLARS(Eigen::Matrix EndmemberMatrix, - Eigen::VectorXf inputVector); - virtual Eigen::VectorXf NNLSGoldfarb(Eigen::Matrix EndmemberMatrix, - Eigen::VectorXf inputVector); - private: mitk::pa::LinearSpectralUnmixingFilter::AlgortihmType algorithmIndex; }; } } #endif // MITKLINEARPHOTOACOUSTICSPECTRALUNMIXINGFILTER_H diff --git a/Modules/PhotoacousticsLib/include/mitkPASpectralUnmixingFilterVigra.h b/Modules/PhotoacousticsLib/include/mitkPASpectralUnmixingFilterVigra.h index b25766cd86..0d08a550c5 100644 --- a/Modules/PhotoacousticsLib/include/mitkPASpectralUnmixingFilterVigra.h +++ b/Modules/PhotoacousticsLib/include/mitkPASpectralUnmixingFilterVigra.h @@ -1,46 +1,60 @@ /*=================================================================== 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) + itkFactorylessNewMacro(Self) + + enum VigraAlgortihmType + { + LARS, + GOLDFARB, + WEIGHTED, + vigratest + }; + + void mitk::pa::SpectralUnmixingFilterVigra::SetAlgorithm(VigraAlgortihmType SetAlgorithmIndex); protected: SpectralUnmixingFilterVigra(); virtual ~SpectralUnmixingFilterVigra(); - + + virtual Eigen::VectorXf SpectralUnmixingAlgorithm(Eigen::Matrix EndmemberMatrix, + Eigen::VectorXf inputVector) override; virtual Eigen::VectorXf NNLSLARS(Eigen::Matrix EndmemberMatrix, Eigen::VectorXf inputVector); virtual Eigen::VectorXf NNLSGoldfarb(Eigen::Matrix EndmemberMatrix, Eigen::VectorXf inputVector); private: + mitk::pa::SpectralUnmixingFilterVigra::VigraAlgortihmType algorithmIndex; }; } } #endif // MITKPHOTOACOUSTICSPECTRALUNMIXINGFILTERVIGRA_H diff --git a/Modules/PhotoacousticsLib/src/SUFilter/mitkPALinearSpectralUnmixingFilter.cpp b/Modules/PhotoacousticsLib/src/SUFilter/mitkPALinearSpectralUnmixingFilter.cpp index ca5575985e..077ff8b7ab 100644 --- a/Modules/PhotoacousticsLib/src/SUFilter/mitkPALinearSpectralUnmixingFilter.cpp +++ b/Modules/PhotoacousticsLib/src/SUFilter/mitkPALinearSpectralUnmixingFilter.cpp @@ -1,247 +1,99 @@ /*=================================================================== 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 "mitkPALinearSpectralUnmixingFilter.h" // Includes for AddEnmemberMatrix #include "mitkPAPropertyCalculator.h" #include // Testing algorithms #include // ImageAccessor #include #include -//vigra -/* -The VIGRA License -================= -(identical to the MIT X11 License) - -Permission is hereby granted, free of charge, to any person -obtaining a copy of this software and associated documentation -files (the "Software"), to deal in the Software without -restriction, including without limitation the rights to use, -copy, modify, merge, publish, distribute, sublicense, and/or -sell copies of the Software, and to permit persons to whom the -Software is furnished to do so, subject to the following -conditions: - -The above copyright notice and this permission notice shall be -included in all copies or substantial portions of the -Software. - -THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND -EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES -OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND -NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT -HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, -WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING -FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR -OTHER DEALINGS IN THE SOFTWARE. -*/ -#include -#include -#include - mitk::pa::LinearSpectralUnmixingFilter::LinearSpectralUnmixingFilter() { } mitk::pa::LinearSpectralUnmixingFilter::~LinearSpectralUnmixingFilter() { } void mitk::pa::LinearSpectralUnmixingFilter::SetAlgorithm(mitk::pa::LinearSpectralUnmixingFilter::AlgortihmType SetAlgorithmIndex) { algorithmIndex = SetAlgorithmIndex; } Eigen::VectorXf mitk::pa::LinearSpectralUnmixingFilter::SpectralUnmixingAlgorithm( Eigen::Matrix EndmemberMatrix, Eigen::VectorXf inputVector) { - - //test other solvers https://eigen.tuxfamily.org/dox/group__TutorialLinearAlgebra.html - Eigen::Vector2f resultVector; if (mitk::pa::LinearSpectralUnmixingFilter::AlgortihmType::colPivHouseholderQr == algorithmIndex) - { - resultVector = EndmemberMatrix.colPivHouseholderQr().solve(inputVector); //works :) - } + resultVector = EndmemberMatrix.colPivHouseholderQr().solve(inputVector); if (mitk::pa::LinearSpectralUnmixingFilter::AlgortihmType::llt == algorithmIndex) - { - resultVector = EndmemberMatrix.llt().solve(inputVector); //works with negativ values (no correct unmixing) - } + resultVector = EndmemberMatrix.llt().solve(inputVector); + if (mitk::pa::LinearSpectralUnmixingFilter::AlgortihmType::householderQr == algorithmIndex) - { - resultVector = EndmemberMatrix.householderQr().solve(inputVector); //works :) - } + resultVector = EndmemberMatrix.householderQr().solve(inputVector); if (mitk::pa::LinearSpectralUnmixingFilter::AlgortihmType::ldlt == algorithmIndex) { mitkThrow() << "algorithm not working"; resultVector = EndmemberMatrix.ldlt().solve(inputVector); //not working because matrix not quadratic(?) } if (mitk::pa::LinearSpectralUnmixingFilter::AlgortihmType::jacobiSvd == algorithmIndex) - { resultVector = EndmemberMatrix.jacobiSvd(Eigen::ComputeFullU | Eigen::ComputeFullV).solve(inputVector); - } if (mitk::pa::LinearSpectralUnmixingFilter::AlgortihmType::fullPivLu == algorithmIndex) - { - resultVector = EndmemberMatrix.fullPivLu().solve(inputVector); //works :) - // look at: https://eigen.tuxfamily.org/dox/classEigen_1_1FullPivLU.html#af563471f6f3283fd10779ef02dd0b748 - /* ALSO FOR ALL QR METHODS - This method just tries to find as good a solution as possible. If you want to check whether a solution exists or if - it is accurate, just call this function to get a result and then compute the error of this result, or use MatrixBase::isApprox() - directly, for instance like this: bool a_solution_exists = (A*result).isApprox(b, precision); - This method avoids dividing by zero, so that the non-existence of a solution doesn't by itself mean that you'll get inf or nan values. - If there exists more than one solution, this method will arbitrarily choose one. - If you need a complete analysis of the space of solutions, take the one solution obtained by this method and add to it - elements of the kernel, as determined by kernel(). - */ - } + resultVector = EndmemberMatrix.fullPivLu().solve(inputVector); if (mitk::pa::LinearSpectralUnmixingFilter::AlgortihmType::householderQr == algorithmIndex) - { - resultVector = EndmemberMatrix.householderQr().solve(inputVector); //works :) - } + resultVector = EndmemberMatrix.householderQr().solve(inputVector); if (mitk::pa::LinearSpectralUnmixingFilter::AlgortihmType::fullPivHouseholderQr == algorithmIndex) - { - resultVector = EndmemberMatrix.fullPivHouseholderQr().solve(inputVector);//works :) - } + resultVector = EndmemberMatrix.fullPivHouseholderQr().solve(inputVector); - //testing new algorithms: if (mitk::pa::LinearSpectralUnmixingFilter::AlgortihmType::test == algorithmIndex) { - //resultVector = NNLSLARS(EndmemberMatrix, inputVector); - resultVector = NNLSGoldfarb(EndmemberMatrix, inputVector); - - //mitkThrow() << "nothing implemented"; - /* - Eigen::LLT lltOfA(EndmemberMatrix); // compute the Cholesky decomposition of A - if (lltOfA.info() == Eigen::NumericalIssue) - MITK_INFO << "not positive semi definite"; - else - MITK_INFO << "semi definite"; - resultVector = EndmemberMatrix.fullPivHouseholderQr().solve(inputVector);//works :) - //resultVector = EndmemberMatrix.completeOrthogonalDecomposition().solve(inputVector); //not a member of Eigen? - //resultVector = EndmemberMatrix.bdcSvd().solve(inputVector); //not a member of Eigen?*/ + mitkThrow() << "404"; } /*double relativeError = (EndmemberMatrix*resultVector - inputVector).norm() / inputVector.norm(); // norm() is L2 norm MITK_INFO << "relativ error: " << relativeError; float accuracyLevel = .1; bool resultIsApprox = inputVector.isApprox(EndmemberMatrix*resultVector, accuracyLevel); MITK_INFO << "IS APPROX RESULT: " << resultIsApprox;*/ - return resultVector; -} - -Eigen::VectorXf mitk::pa::LinearSpectralUnmixingFilter::NNLSLARS( - Eigen::Matrix EndmemberMatrix, Eigen::VectorXf inputVector) -{ - using namespace vigra; - using namespace vigra::linalg; - - int numberOfWavelengths = EndmemberMatrix.rows(); - int numberOfChromophores = EndmemberMatrix.cols(); - - std::vector A_data; - std::vector B_data; - for (int i = 0; i < numberOfWavelengths; ++i) - { - B_data.push_back(inputVector(i)); - for (int j = 0; j < numberOfChromophores; ++j) - { - A_data.push_back(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) - nonnegativeLeastSquares(A, b, x); - - Eigen::VectorXf result(numberOfChromophores); - for (int k = 0; k < numberOfChromophores; ++k) - { - float test = x(k, 0); - result[k] = test; - } - return result; + return resultVector; } -Eigen::VectorXf mitk::pa::LinearSpectralUnmixingFilter::NNLSGoldfarb( - Eigen::Matrix EndmemberMatrix, Eigen::VectorXf inputVector) -{ - using namespace vigra; - using namespace vigra::linalg; - - int numberOfWavelengths = EndmemberMatrix.rows(); - int numberOfChromophores = EndmemberMatrix.cols(); - - std::vector A_data; - std::vector B_data; - - for (int i = 0; i < numberOfWavelengths; ++i) - { - B_data.push_back(inputVector(i)); - for (int j = 0; j < numberOfChromophores; ++j) - { - A_data.push_back(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)); - - 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); - - Eigen::VectorXf result(numberOfChromophores); - for (int k = 0; k < numberOfChromophores; ++k) - { - float test = x(k, 0); - result[k] = test; - } - return result; -} +/* look at: https://eigen.tuxfamily.org/dox/classEigen_1_1FullPivLU.html#af563471f6f3283fd10779ef02dd0b748 +This method just tries to find as good a solution as possible. If you want to check whether a solution exists or if +it is accurate, just call this function to get a result and then compute the error of this result, or use MatrixBase::isApprox() +directly, for instance like this: bool a_solution_exists = (A*result).isApprox(b, precision); +This method avoids dividing by zero, so that the non-existence of a solution doesn't by itself mean that you'll get inf or nan values. +If there exists more than one solution, this method will arbitrarily choose one. +If you need a complete analysis of the space of solutions, take the one solution obtained by this method and add to it +elements of the kernel, as determined by kernel().*/ diff --git a/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingFilterVigra.cpp b/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingFilterVigra.cpp index 77f5b250a3..979b47deb6 100644 --- a/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingFilterVigra.cpp +++ b/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingFilterVigra.cpp @@ -1,157 +1,183 @@ /*=================================================================== 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 /* The VIGRA License ================= (identical to the MIT X11 License) Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions: The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. */ #include #include #include mitk::pa::SpectralUnmixingFilterVigra::SpectralUnmixingFilterVigra() { } mitk::pa::SpectralUnmixingFilterVigra::~SpectralUnmixingFilterVigra() { } +void mitk::pa::SpectralUnmixingFilterVigra::SetAlgorithm(mitk::pa::SpectralUnmixingFilterVigra::VigraAlgortihmType SetAlgorithmIndex) +{ + algorithmIndex = SetAlgorithmIndex; +} + +Eigen::VectorXf mitk::pa::SpectralUnmixingFilterVigra::SpectralUnmixingAlgorithm( + Eigen::Matrix EndmemberMatrix, Eigen::VectorXf inputVector) +{ + Eigen::Vector2f resultVector; + + if (mitk::pa::SpectralUnmixingFilterVigra::VigraAlgortihmType::LARS == algorithmIndex) + resultVector = NNLSLARS(EndmemberMatrix, inputVector); + + if (mitk::pa::SpectralUnmixingFilterVigra::VigraAlgortihmType::GOLDFARB == algorithmIndex) + resultVector = NNLSGoldfarb(EndmemberMatrix, inputVector); + + if (mitk::pa::SpectralUnmixingFilterVigra::VigraAlgortihmType::WEIGHTED == algorithmIndex) + mitkThrow() << "nothing implemented"; + + if (mitk::pa::SpectralUnmixingFilterVigra::VigraAlgortihmType::vigratest == algorithmIndex) + mitkThrow() << "nothing implemented"; + + return resultVector; +} + + Eigen::VectorXf mitk::pa::SpectralUnmixingFilterVigra::NNLSLARS( Eigen::Matrix EndmemberMatrix, Eigen::VectorXf inputVector) { using namespace vigra; using namespace vigra::linalg; int numberOfWavelengths = EndmemberMatrix.rows(); int numberOfChromophores = EndmemberMatrix.cols(); std::vector A_data; std::vector B_data; for (int i = 0; i < numberOfWavelengths; ++i) { B_data.push_back(inputVector(i)); for (int j = 0; j < numberOfChromophores; ++j) { A_data.push_back(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) nonnegativeLeastSquares(A, b, x); Eigen::VectorXf result(numberOfChromophores); for (int k = 0; k < numberOfChromophores; ++k) { float test = x(k, 0); result[k] = test; } return result; } Eigen::VectorXf mitk::pa::SpectralUnmixingFilterVigra::NNLSGoldfarb( Eigen::Matrix EndmemberMatrix, Eigen::VectorXf inputVector) { using namespace vigra; using namespace vigra::linalg; int numberOfWavelengths = EndmemberMatrix.rows(); int numberOfChromophores = EndmemberMatrix.cols(); std::vector A_data; std::vector B_data; for (int i = 0; i < numberOfWavelengths; ++i) { B_data.push_back(inputVector(i)); for (int j = 0; j < numberOfChromophores; ++j) { A_data.push_back(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)); 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); Eigen::VectorXf result(numberOfChromophores); for (int k = 0; k < numberOfChromophores; ++k) { float test = x(k, 0); result[k] = test; } return result; } 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 bfdf1dad0a..7646de598e 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,256 +1,319 @@ /*=================================================================== 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); connect(m_Controls.ButtonAddWavelength, &QPushButton::clicked, this, &SpectralUnmixing::Wavelength); connect(m_Controls.ButtonClearWavelength, &QPushButton::clicked, this, &SpectralUnmixing::ClearWavelength); } // Adds Wavelength @ Plugin with button void SpectralUnmixing::Wavelength() { if (m_Wavelengths.empty()) { size = 0; } unsigned int wavelength = m_Controls.spinBoxAddWavelength->value(); m_Wavelengths.push_back(wavelength); size += 1; // size implemented like this because '.size' is const MITK_INFO << "ALL WAVELENGTHS: "; for (int i = 0; i < size; ++i) { MITK_INFO << m_Wavelengths[i] << "nm"; } } 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); } 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(); - auto m_SpectralUnmixingFilter = mitk::pa::LinearSpectralUnmixingFilter::New(); - - m_SpectralUnmixingFilter->SetInput(image); - // Set Algortihm to filter auto qs = m_Controls.QComboBoxAlgorithm->currentText(); std::string Algorithm = qs.toUtf8().constData(); + auto m_SpectralUnmixingFilter = mitk::pa::LinearSpectralUnmixingFilter::New(); + m_SpectralUnmixingFilter->SetInput(image); + if (Algorithm == "householderQr") + { + //auto m_SpectralUnmixingFilter = mitk::pa::LinearSpectralUnmixingFilter::New(); + //m_SpectralUnmixingFilter->SetInput(image); m_SpectralUnmixingFilter->SetAlgorithm(mitk::pa::LinearSpectralUnmixingFilter::AlgortihmType::householderQr); + } else if (Algorithm == "ldlt") + { + //auto m_SpectralUnmixingFilter = mitk::pa::LinearSpectralUnmixingFilter::New(); + //m_SpectralUnmixingFilter->SetInput(image); m_SpectralUnmixingFilter->SetAlgorithm(mitk::pa::LinearSpectralUnmixingFilter::AlgortihmType::ldlt); + } else if (Algorithm == "llt") + { + //auto m_SpectralUnmixingFilter = mitk::pa::LinearSpectralUnmixingFilter::New(); + //m_SpectralUnmixingFilter->SetInput(image); m_SpectralUnmixingFilter->SetAlgorithm(mitk::pa::LinearSpectralUnmixingFilter::AlgortihmType::llt); + } else if (Algorithm == "colPivHouseholderQr") + { + //auto m_SpectralUnmixingFilter = mitk::pa::LinearSpectralUnmixingFilter::New(); + //m_SpectralUnmixingFilter->SetInput(image); m_SpectralUnmixingFilter->SetAlgorithm(mitk::pa::LinearSpectralUnmixingFilter::AlgortihmType::colPivHouseholderQr); + } else if (Algorithm == "jacobiSvd") + { + //auto m_SpectralUnmixingFilter = mitk::pa::LinearSpectralUnmixingFilter::New(); + //m_SpectralUnmixingFilter->SetInput(image); m_SpectralUnmixingFilter->SetAlgorithm(mitk::pa::LinearSpectralUnmixingFilter::AlgortihmType::jacobiSvd); + } else if (Algorithm == "fullPivLu") + { + //auto m_SpectralUnmixingFilter = mitk::pa::LinearSpectralUnmixingFilter::New(); + //m_SpectralUnmixingFilter->SetInput(image); m_SpectralUnmixingFilter->SetAlgorithm(mitk::pa::LinearSpectralUnmixingFilter::AlgortihmType::fullPivLu); + } else if (Algorithm == "fullPivHouseholderQr") + { + //auto m_SpectralUnmixingFilter = mitk::pa::LinearSpectralUnmixingFilter::New(); + //m_SpectralUnmixingFilter->SetInput(image); m_SpectralUnmixingFilter->SetAlgorithm(mitk::pa::LinearSpectralUnmixingFilter::AlgortihmType::fullPivHouseholderQr); + } else if (Algorithm == "test") + { + //auto m_SpectralUnmixingFilter = mitk::pa::LinearSpectralUnmixingFilter::New(); + //m_SpectralUnmixingFilter->SetInput(image); m_SpectralUnmixingFilter->SetAlgorithm(mitk::pa::LinearSpectralUnmixingFilter::AlgortihmType::test); + } + + else if (Algorithm == "Lars") + { + auto m_SpectralUnmixingFilter = mitk::pa::SpectralUnmixingFilterVigra::New(); + //m_SpectralUnmixingFilter->SetInput(image); + m_SpectralUnmixingFilter->SetAlgorithm(mitk::pa::SpectralUnmixingFilterVigra::VigraAlgortihmType::LARS); + } + + else if (Algorithm == "Goldfarb") + { + auto m_SpectralUnmixingFilter = mitk::pa::SpectralUnmixingFilterVigra::New(); + //m_SpectralUnmixingFilter->SetInput(image); + m_SpectralUnmixingFilter->SetAlgorithm(mitk::pa::SpectralUnmixingFilterVigra::VigraAlgortihmType::GOLDFARB); + } + + else if (Algorithm == "weighted") + { + auto m_SpectralUnmixingFilter = mitk::pa::SpectralUnmixingFilterVigra::New(); + //m_SpectralUnmixingFilter->SetInput(image); + m_SpectralUnmixingFilter->SetAlgorithm(mitk::pa::SpectralUnmixingFilterVigra::VigraAlgortihmType::WEIGHTED); + } + + else if (Algorithm == "vigratest") + { + auto m_SpectralUnmixingFilter = mitk::pa::SpectralUnmixingFilterVigra::New(); + //m_SpectralUnmixingFilter->SetInput(image); + m_SpectralUnmixingFilter->SetAlgorithm(mitk::pa::SpectralUnmixingFilterVigra::VigraAlgortihmType::vigratest); + } + else mitkThrow() << "ALGORITHM ERROR!"; // Wavelength implementation into filter for (unsigned int imageIndex = 0; imageIndex < m_Wavelengths.size(); imageIndex++) { unsigned int wavelength = m_Wavelengths[imageIndex]; m_SpectralUnmixingFilter->AddWavelength(wavelength); } // Checking which chromophores wanted for SU if none throw exeption! unsigned int numberofChromophores = 0; DeOxbool = m_Controls.checkBoxDeOx->isChecked(); Oxbool = m_Controls.checkBoxOx->isChecked(); if (DeOxbool || Oxbool) { MITK_INFO << "CHOSEN CHROMOPHORES:"; } if (Oxbool) { numberofChromophores += 1; MITK_INFO << "- Oxyhemoglobin"; // Set chromophore Oxyhemoglobon: m_SpectralUnmixingFilter->AddChromophore( mitk::pa::SpectralUnmixingFilterBase::ChromophoreType::OXYGENATED_HEMOGLOBIN); } if (DeOxbool) { numberofChromophores += 1; MITK_INFO << "- Deoxygenated hemoglobin"; // Set chromophore Deoxygenated hemoglobin: m_SpectralUnmixingFilter->AddChromophore( mitk::pa::SpectralUnmixingFilterBase::ChromophoreType::DEOXYGENATED_HEMOGLOBIN); } if (numberofChromophores == 0) { mitkThrow() << "PRESS 'IGNORE' AND CHOOSE A CHROMOPHORE!"; } MITK_INFO << "Updating Filter..."; m_SpectralUnmixingFilter->Update(); // Write Output images to Data Storage if (Oxbool) { mitk::Image::Pointer HbO2 = m_SpectralUnmixingFilter->GetOutput(0); 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(1); mitk::DataNode::Pointer dataNodeHb = mitk::DataNode::New(); dataNodeHb->SetData(Hb); dataNodeHb->SetName("Hb " + Algorithm); this->GetDataStorage()->Add(dataNodeHb); } //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 ..."; // Initialize pipeline from SU filter class to SO2 class auto m_sO2 = mitk::pa::SpectralUnmixingSO2::New(); 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"); this->GetDataStorage()->Add(dataNodesO2); MITK_INFO << "[DONE]"; } mitk::RenderingManager::GetInstance()->InitializeViewsByBoundingObjects(this->GetDataStorage()); MITK_INFO << "Adding images to DataStorage...[DONE]"; } } }