diff --git a/Modules/PhotoacousticsLib/include/mitkPALinearSpectralUnmixingFilter.h b/Modules/PhotoacousticsLib/include/mitkPALinearSpectralUnmixingFilter.h index 3374d20478..d4b9c6d2a0 100644 --- a/Modules/PhotoacousticsLib/include/mitkPALinearSpectralUnmixingFilter.h +++ b/Modules/PhotoacousticsLib/include/mitkPALinearSpectralUnmixingFilter.h @@ -1,67 +1,58 @@ /*=================================================================== 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 -#include "mitkPAPropertyCalculator.h" - 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(); virtual Eigen::VectorXf SpectralUnmixingAlgorithm(Eigen::Matrix EndmemberMatrix, Eigen::VectorXf inputVector) override; private: mitk::pa::LinearSpectralUnmixingFilter::AlgortihmType algorithmIndex; - - virtual Eigen::Matrix CalculateEndmemberMatrix( - std::vector m_Chromophore, std::vector m_Wavelength); - PropertyCalculator::Pointer m_PropertyCalculatorEigen; - - virtual float propertyElement(mitk::pa::PropertyCalculator::ChromophoreType, int wavelength); - }; } } #endif // MITKLINEARPHOTOACOUSTICSPECTRALUNMIXINGFILTER_H diff --git a/Modules/PhotoacousticsLib/include/mitkPASpectralUnmixingFilterBase.h b/Modules/PhotoacousticsLib/include/mitkPASpectralUnmixingFilterBase.h index 5c84cb9776..5b3658ccef 100644 --- a/Modules/PhotoacousticsLib/include/mitkPASpectralUnmixingFilterBase.h +++ b/Modules/PhotoacousticsLib/include/mitkPASpectralUnmixingFilterBase.h @@ -1,76 +1,76 @@ /*=================================================================== 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 { class MITKPHOTOACOUSTICSLIB_EXPORT SpectralUnmixingFilterBase : public mitk::ImageToImageFilter { public: mitkClassMacro(SpectralUnmixingFilterBase, mitk::ImageToImageFilter); // Contains all (so far) possible chromophores // Void to creat m_vector of all chosen chromophores with push back method void AddChromophore(mitk::pa::PropertyCalculator::ChromophoreType); // Void to creat m_vector of all wavelengths with push back method void AddWavelength(int wavelength); protected: SpectralUnmixingFilterBase(); virtual ~SpectralUnmixingFilterBase(); std::vector m_Chromophore; std::vector m_Wavelength; virtual Eigen::VectorXf SpectralUnmixingAlgorithm(Eigen::Matrix EndmemberMatrix, Eigen::VectorXf inputVector) = 0; private: // Void checking precondtions possibly throwing exeptions virtual void CheckPreConditions(unsigned int size, unsigned int NumberOfInputImages, const float* inputDataArray); virtual void GenerateData() override; virtual void InitializeOutputs(); - // Void to creat Eigen::Matrix of all absorbtions - // @ specific wavelength (columns) of chromophores (rows) - Eigen::Matrix AddEndmemberMatrix(); + virtual Eigen::Matrix CalculateEndmemberMatrix( + std::vector m_Chromophore, std::vector m_Wavelength); + PropertyCalculator::Pointer m_PropertyCalculatorEigen; - PropertyCalculator::Pointer m_PropertyCalculator; + virtual float propertyElement(mitk::pa::PropertyCalculator::ChromophoreType, int wavelength); }; } } #endif // MITKPHOTOACOUSTICSPECTRALUNMIXINGFILTERBASE_ diff --git a/Modules/PhotoacousticsLib/src/SUFilter/mitkPALinearSpectralUnmixingFilter.cpp b/Modules/PhotoacousticsLib/src/SUFilter/mitkPALinearSpectralUnmixingFilter.cpp index 991d88a192..c526b2b9f1 100644 --- a/Modules/PhotoacousticsLib/src/SUFilter/mitkPALinearSpectralUnmixingFilter.cpp +++ b/Modules/PhotoacousticsLib/src/SUFilter/mitkPALinearSpectralUnmixingFilter.cpp @@ -1,147 +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" -// Includes for AddEnmemberMatrix -#include "mitkPAPropertyCalculator.h" -#include - // Testing algorithms #include // ImageAccessor #include #include - mitk::pa::LinearSpectralUnmixingFilter::LinearSpectralUnmixingFilter() { - m_PropertyCalculatorEigen = mitk::pa::PropertyCalculator::New(); } 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; 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); //not working because matrix not quadratic(?) + 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) { - Eigen::Matrix Test = CalculateEndmemberMatrix(m_Chromophore, m_Wavelength); - resultVector = Test.householderQr().solve(inputVector); - //mitkThrow() << "404"; + 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; } - -/* 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().*/ - -//Matrix with #chromophores colums and #wavelengths rows -//so Matrix Element (i,j) contains the Absorbtion of chromophore j @ wavelength i -Eigen::Matrix mitk::pa::LinearSpectralUnmixingFilter::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]); - } - 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]!"; - return EndmemberMatrixEigen; -} - -float mitk::pa::LinearSpectralUnmixingFilter::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/mitkPASpectralUnmixingFilterBase.cpp b/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingFilterBase.cpp index 64342d05f4..82e29a6ced 100644 --- a/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingFilterBase.cpp +++ b/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingFilterBase.cpp @@ -1,218 +1,225 @@ /*=================================================================== 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 // For testing reasons #include mitk::pa::SpectralUnmixingFilterBase::SpectralUnmixingFilterBase() { this->SetNumberOfIndexedOutputs(2); for (unsigned int i = 0; iSetNthOutput(i, mitk::Image::New()); } - m_PropertyCalculator = mitk::pa::PropertyCalculator::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 = AddEndmemberMatrix(); + 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"); //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"; } } MITK_INFO << "GENERATING DATA...[DONE]"; //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::AddEndmemberMatrix() +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); - Eigen::Matrix EndmemberMatrix(numberOfWavelengths, numberOfChromophores); - - //loop over i rows (Wavelengths) - for(unsigned int i = 0; i < numberOfWavelengths; ++i) + for (unsigned int j = 0; j < numberOfChromophores; ++j) { - //loop over j columns (Chromophores) - for (unsigned int j = 0; j < numberOfChromophores; ++j) + if (m_Chromophore[j] == mitk::pa::PropertyCalculator::ChromophoreType::OXYGENATED) { - //writes @ Matrix element (i,j) the absorbtion wavelength of the propertycalculator.cpp - //'i+1' because MapType is defined different then ChromophoreType - EndmemberMatrix(i,j)= m_PropertyCalculator->GetAbsorptionForWavelength( - static_cast(m_Chromophore[j]+1), m_Wavelength[i]); - - //* Tests to see what gets written in the Matrix: - //auto testtype = m_PropertyCalculator->GetAbsorptionForWavelength( - //static_cast(m_Chromophore[j]+1), m_Wavelength[i]); - //MITK_INFO << "TEST_TYPE Matrix: " << typeid(EndmemberMatrix(i,j)).name(); - MITK_INFO << "ChromophoreType: " << static_cast(m_Chromophore[j]+1); - MITK_INFO << "wavelength: " << m_Wavelength[i]; - //MITK_INFO << "Matrixelement: (" << i << ", " << j << ") Absorbtion: " << EndmemberMatrix(i, j);/**/ - - if (EndmemberMatrix(i, j) == 0) - { - m_Wavelength.clear(); - mitkThrow() << "WAVELENGTH "<< m_Wavelength[i] << "nm NOT SUPPORTED!"; - } + 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]); + } + 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]!"; - return EndmemberMatrix; + //MITK_INFO << "GENERATING ENMEMBERMATRIX [DONE]!"; + 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; }