diff --git a/Modules/PhotoacousticsLib/include/mitkPASpectralUnmixingFilterVigra.h b/Modules/PhotoacousticsLib/include/mitkPASpectralUnmixingFilterVigra.h index 809d29e135..4f4db58c31 100644 --- a/Modules/PhotoacousticsLib/include/mitkPASpectralUnmixingFilterVigra.h +++ b/Modules/PhotoacousticsLib/include/mitkPASpectralUnmixingFilterVigra.h @@ -1,64 +1,56 @@ /*=================================================================== 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, 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); - virtual Eigen::VectorXf weightedLeastSquares(Eigen::Matrix EndmemberMatrix, - Eigen::VectorXf inputVector); - virtual Eigen::VectorXf LS(Eigen::Matrix EndmemberMatrix, - Eigen::VectorXf inputVector); private: mitk::pa::SpectralUnmixingFilterVigra::VigraAlgortihmType algorithmIndex; }; } } #endif // MITKPHOTOACOUSTICSPECTRALUNMIXINGFILTERVIGRA_H diff --git a/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingFilterVigra.cpp b/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingFilterVigra.cpp index 9637dc0832..1107d33bdb 100644 --- a/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingFilterVigra.cpp +++ b/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingFilterVigra.cpp @@ -1,266 +1,148 @@ /*=================================================================== 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::VectorXf 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) - resultVector = weightedLeastSquares(EndmemberMatrix, inputVector); - - if (mitk::pa::SpectralUnmixingFilterVigra::VigraAlgortihmType::vigratest == algorithmIndex) - resultVector = LS(EndmemberMatrix, inputVector); - //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::LS( - 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; + std::vector weightsvec = { 1,1,1,1,1 }; //am besten über GUI Eingabe festlegen 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) - linearSolve(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; + if (mitk::pa::SpectralUnmixingFilterVigra::VigraAlgortihmType::LARS == algorithmIndex) + nonnegativeLeastSquares(A, b, x); - for (int i = 0; i < numberOfWavelengths; ++i) + else if (mitk::pa::SpectralUnmixingFilterVigra::VigraAlgortihmType::GOLDFARB == algorithmIndex) { - B_data.push_back(inputVector(i)); - for (int j = 0; j < numberOfChromophores; ++j) - { - A_data.push_back(EndmemberMatrix(i, j)); - } + 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); } - 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) + else if (mitk::pa::SpectralUnmixingFilterVigra::VigraAlgortihmType::WEIGHTED == algorithmIndex) { - float test = x(k, 0); - result[k] = test; - } - return result; -} - -Eigen::VectorXf mitk::pa::SpectralUnmixingFilterVigra::weightedLeastSquares( - Eigen::Matrix EndmemberMatrix, Eigen::VectorXf inputVector) -{ - using namespace vigra; - using namespace vigra::linalg; + std::vector weightsvec = { 1,1,1,1,1 }; //am besten über GUI Eingabe festlegen + const double* weightsdat = weightsvec.data(); + vigra::Matrix weigths(Shape2(numberOfWavelengths, 1), weightsdat); - int numberOfWavelengths = EndmemberMatrix.rows(); - int numberOfChromophores = EndmemberMatrix.cols(); - - std::vector A_data; - std::vector B_data; - std::vector weightsvec = { 1,1,1,1,1 }; //am besten über GUI Eingabe festlegen - - 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)); - } + vigra::linalg::weightedLeastSquares(A, b, weigths, x); } - const double* Adat = A_data.data(); - const double* Bdat = B_data.data(); - const double* weightsdat = weightsvec.data(); - - vigra::Matrix A(Shape2(numberOfWavelengths, numberOfChromophores), Adat); - vigra::Matrix b(Shape2(numberOfWavelengths, 1), Bdat); - vigra::Matrix weigths(Shape2(numberOfWavelengths, 1), weightsdat); - vigra::Matrix x(Shape2(numberOfChromophores, 1)); + else if (mitk::pa::SpectralUnmixingFilterVigra::VigraAlgortihmType::vigratest == algorithmIndex) + linearSolve(A, b, x); + //mitkThrow() << "nothing implemented"; - vigra::linalg::weightedLeastSquares(A, b, weigths, x); + else + mitkThrow() << "404 VIGRA ALGORITHM NOT FOUND"; Eigen::VectorXf result(numberOfChromophores); for (int k = 0; k < numberOfChromophores; ++k) { float test = x(k, 0); result[k] = test; } - return result; + + return resultVector; }