diff --git a/Modules/PhotoacousticsLib/include/mitkPAPropertyCalculator.h b/Modules/PhotoacousticsLib/include/mitkPAPropertyCalculator.h index 9fd1417ca7..89bf0be93d 100644 --- a/Modules/PhotoacousticsLib/include/mitkPAPropertyCalculator.h +++ b/Modules/PhotoacousticsLib/include/mitkPAPropertyCalculator.h @@ -1,79 +1,79 @@ /*=================================================================== 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 MITKPHOTOACOUSTICPROPERTYCALCULATOR_H #define MITKPHOTOACOUSTICPROPERTYCALCULATOR_H #include "MitkPhotoacousticsLibExports.h" //Includes for smart pointer usage #include "mitkCommon.h" #include "itkLightObject.h" namespace mitk { namespace pa { class MITKPHOTOACOUSTICSLIB_EXPORT PropertyCalculator : public itk::LightObject { public: mitkClassMacroItkParent(PropertyCalculator, itk::LightObject) itkFactorylessNewMacro(Self) struct Properties { double mua; double mus; double g; }; enum TissueType { AIR = 1, BLOOD = 2, EPIDERMIS = 3, FAT = 4, STANDARD_TISSUE = 5 }; enum ChromophoreType { OXYGENATED = 1, DEOXYGENATED = 2, WATER = 3, FATTY = 4, MELANIN = 5, ONEENDMEMBER = 6 }; double GetAbsorptionForWavelength( - ChromophoreType ChromophoreType, int wavelength); + ChromophoreType chromophoreType, int wavelength); Properties CalculatePropertyForSpecificWavelength( TissueType tissueType, int wavelength, double oxygenSaturatonInFraction = 0); protected: PropertyCalculator(); ~PropertyCalculator() override; bool m_Valid = false; std::map> m_SpectralLibMap; }; } } #endif // MITKPHOTOACOUSTICPROPERTYCALCULATOR_H diff --git a/Modules/PhotoacousticsLib/src/SUFilter/mitkPALinearSpectralUnmixingFilter.cpp b/Modules/PhotoacousticsLib/src/SUFilter/mitkPALinearSpectralUnmixingFilter.cpp index c526b2b9f1..4a2bf6c2b7 100644 --- a/Modules/PhotoacousticsLib/src/SUFilter/mitkPALinearSpectralUnmixingFilter.cpp +++ b/Modules/PhotoacousticsLib/src/SUFilter/mitkPALinearSpectralUnmixingFilter.cpp @@ -1,83 +1,83 @@ /*=================================================================== 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" // Testing algorithms #include // ImageAccessor #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) { - Eigen::Vector2f resultVector; + Eigen::VectorXf resultVector; if (mitk::pa::LinearSpectralUnmixingFilter::AlgortihmType::colPivHouseholderQr == algorithmIndex) resultVector = EndmemberMatrix.colPivHouseholderQr().solve(inputVector); if (mitk::pa::LinearSpectralUnmixingFilter::AlgortihmType::llt == algorithmIndex) resultVector = EndmemberMatrix.llt().solve(inputVector); if (mitk::pa::LinearSpectralUnmixingFilter::AlgortihmType::householderQr == algorithmIndex) resultVector = EndmemberMatrix.householderQr().solve(inputVector); if (mitk::pa::LinearSpectralUnmixingFilter::AlgortihmType::ldlt == algorithmIndex) { mitkThrow() << "algorithm not working"; resultVector = EndmemberMatrix.ldlt().solve(inputVector); } 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); if (mitk::pa::LinearSpectralUnmixingFilter::AlgortihmType::householderQr == algorithmIndex) resultVector = EndmemberMatrix.householderQr().solve(inputVector); if (mitk::pa::LinearSpectralUnmixingFilter::AlgortihmType::fullPivHouseholderQr == algorithmIndex) resultVector = EndmemberMatrix.fullPivHouseholderQr().solve(inputVector); if (mitk::pa::LinearSpectralUnmixingFilter::AlgortihmType::test == algorithmIndex) { mitkThrow() << "404 NO ALGORITHM IMPLEMENTED"; } /*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;*/ //MITK_INFO << "Result vector: " << resultVector; return resultVector; } diff --git a/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingFilterBase.cpp b/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingFilterBase.cpp index 468619d541..f663ac9abf 100644 --- a/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingFilterBase.cpp +++ b/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingFilterBase.cpp @@ -1,225 +1,227 @@ /*=================================================================== 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 mitk::pa::SpectralUnmixingFilterBase::SpectralUnmixingFilterBase() { this->SetNumberOfIndexedOutputs(4);// find solution 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::GenerateData() { MITK_INFO << "GENERATING DATA.."; // Get input image mitk::Image::Pointer input = GetInput(0); unsigned int xDim = input->GetDimensions()[0]; unsigned int yDim = input->GetDimensions()[1]; unsigned int zDim = input->GetDimensions()[2]; unsigned int size = xDim * yDim * zDim; MITK_INFO << "x dimension: " << xDim; MITK_INFO << "y dimension: " << yDim; MITK_INFO << "z dimension: " << zDim; InitializeOutputs(); auto EndmemberMatrix = CalculateEndmemberMatrix(m_Chromophore, m_Wavelength); //Eigen Endmember Matrix // Copy input image into array mitk::ImageReadAccessor readAccess(input); const float* inputDataArray = ((const float*)readAccess.GetData()); CheckPreConditions(size, zDim, inputDataArray); /* READ OUT INPUTARRAY MITK_INFO << "Info Array:"; int numberOfPixels= 6; for (int i=0; i< numberOfPixels; ++i) MITK_INFO << inputDataArray[i];/**/ // test to see pixel values @ txt file - //ofstream myfile; - //myfile.open("PASpectralUnmixingPixelValues.txt"); + ofstream myfile; + myfile.open("PASpectralUnmixingPixelValues.txt"); //loop over every pixel @ x,y plane for (unsigned int x = 0; x < xDim; x++) { for (unsigned int y = 0; y < yDim; y++) { Eigen::VectorXf inputVector(zDim); for (unsigned int z = 0; z < zDim; z++) { // Get pixel value of pixel x,y @ wavelength z unsigned int pixelNumber = (xDim*yDim*z) + x * yDim + y; auto pixel = inputDataArray[pixelNumber]; //MITK_INFO << "Pixel_values: " << pixel; // dummy values for pixel for testing purposes //float pixel = rand(); //write all wavelength absorbtion values for one(!) pixel to a vector inputVector[z] = pixel; } Eigen::VectorXf resultVector = SpectralUnmixingAlgorithm(EndmemberMatrix, inputVector); //Eigen::VectorXf resultVector = SpectralUnmixingAlgorithm(m_Chromophore, inputVector); // write output for (int outputIdx = 0; outputIdx < GetNumberOfIndexedOutputs(); ++outputIdx) { auto output = GetOutput(outputIdx); mitk::ImageWriteAccessor writeOutput(output); float* writeBuffer = (float *)writeOutput.GetData(); writeBuffer[x*yDim + y] = resultVector[outputIdx]; } - //myfile << "Input Pixel(x,y): " << x << "," << y << "\n" << inputVector << "\n"; - //myfile << "Result: " << "\n HbO2: " << resultVector[0] << "\n Hb: " << resultVector[1] <<"\n"; + myfile << "Input Pixel(x,y): " << x << "," << y << "\n" << inputVector << "\n"; + myfile << "Result: " << "\n HbO2: " << resultVector[0] << "\n Hb: " << resultVector[1] << "\n Mel: " + << resultVector[2] << "Result vector size" << resultVector.size() << "\n"; } } MITK_INFO << "GENERATING DATA...[DONE]"; - //myfile.close(); + myfile.close(); } void mitk::pa::SpectralUnmixingFilterBase::CheckPreConditions(unsigned int size, unsigned int NumberOfInputImages, const float* inputDataArray) { // Checking if number of Inputs == added wavelengths if (m_Wavelength.size() != NumberOfInputImages) mitkThrow() << "CHECK INPUTS! WAVELENGTHERROR"; // Checking if number of wavelengths >= number of chromophores if (m_Chromophore.size() > m_Wavelength.size()) mitkThrow() << "PRESS 'IGNORE' AND ADD MORE WAVELENGTHS!"; // Checking if pixel type is float int maxPixel = size; for (int i = 0; i < maxPixel; ++i) { if (typeid(inputDataArray[i]).name() != typeid(float).name()) { mitkThrow() << "PIXELTYPE ERROR! FLOAT 32 REQUIRED"; } else continue; } MITK_INFO << "CHECK PRECONDITIONS ...[DONE]"; } void mitk::pa::SpectralUnmixingFilterBase::InitializeOutputs() { unsigned int numberOfInputs = GetNumberOfIndexedInputs(); unsigned int numberOfOutputs = GetNumberOfIndexedOutputs(); //MITK_INFO << "Inputs: " << numberOfInputs << " Outputs: " << numberOfOutputs; // Set dimensions (2) and pixel type (float) for output mitk::PixelType pixelType = mitk::MakeScalarPixelType(); const int NUMBER_OF_SPATIAL_DIMENSIONS = 2; auto* dimensions = new unsigned int[NUMBER_OF_SPATIAL_DIMENSIONS]; for(unsigned int dimIdx=0; dimIdxGetDimensions()[dimIdx]; } // Initialize numberOfOutput pictures with pixel type and dimensions for (unsigned int outputIdx = 0; outputIdx < numberOfOutputs; outputIdx++) { GetOutput(outputIdx)->Initialize(pixelType, NUMBER_OF_SPATIAL_DIMENSIONS, dimensions); } } //Matrix with #chromophores colums and #wavelengths rows //so Matrix Element (i,j) contains the Absorbtion of chromophore j @ wavelength i 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) { if (m_Chromophore[j] == mitk::pa::PropertyCalculator::ChromophoreType::OXYGENATED) { for (unsigned int i = 0; i < numberOfWavelengths; ++i) EndmemberMatrixEigen(i, j) = propertyElement(m_Chromophore[j], m_Wavelength[i]); } else if (m_Chromophore[j] == mitk::pa::PropertyCalculator::ChromophoreType::DEOXYGENATED) { for (unsigned int i = 0; i < numberOfWavelengths; ++i) EndmemberMatrixEigen(i, j) = propertyElement(m_Chromophore[j], m_Wavelength[i]); } else if (m_Chromophore[j] == mitk::pa::PropertyCalculator::ChromophoreType::MELANIN) { for (unsigned int i = 0; i < numberOfWavelengths; ++i) { EndmemberMatrixEigen(i, j) = propertyElement(m_Chromophore[j], m_Wavelength[i]); - MITK_INFO << "MELANIN." << EndmemberMatrixEigen(i, j); + MITK_INFO << "MELANIN " << EndmemberMatrixEigen(i, j); } } else if (m_Chromophore[j] == mitk::pa::PropertyCalculator::ChromophoreType::ONEENDMEMBER) { for (unsigned int i = 0; i < numberOfWavelengths; ++i) EndmemberMatrixEigen(i, j) = 1; } else mitkThrow() << "404 CHROMOPHORE NOT FOUND!"; } - //MITK_INFO << "GENERATING ENMEMBERMATRIX [DONE]!"; + MITK_INFO << "GENERATING ENMEMBERMATRIX [DONE]!"; + MITK_INFO << "endmember matrix: " << EndmemberMatrixEigen; return EndmemberMatrixEigen; } float mitk::pa::SpectralUnmixingFilterBase::propertyElement(mitk::pa::PropertyCalculator::ChromophoreType Chromophore, int Wavelength) { float value = m_PropertyCalculatorEigen->GetAbsorptionForWavelength(Chromophore, Wavelength); if (value == 0) mitkThrow() << "WAVELENGTH " << Wavelength << "nm NOT SUPPORTED!"; return value; } diff --git a/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingFilterVigra.cpp b/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingFilterVigra.cpp index 7ea8a85be2..9637dc0832 100644 --- a/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingFilterVigra.cpp +++ b/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingFilterVigra.cpp @@ -1,266 +1,266 @@ /*=================================================================== 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; + 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; 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; 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; } Eigen::VectorXf mitk::pa::SpectralUnmixingFilterVigra::weightedLeastSquares( 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(); 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)); vigra::linalg::weightedLeastSquares(A, b, weigths, x); Eigen::VectorXf result(numberOfChromophores); for (int k = 0; k < numberOfChromophores; ++k) { float test = x(k, 0); result[k] = test; } return result; }