diff --git a/Modules/PhotoacousticsLib/include/mitkPASpectralUnmixingFilterBase.h b/Modules/PhotoacousticsLib/include/mitkPASpectralUnmixingFilterBase.h index 8c6785efcf..ed8b8b64bf 100644 --- a/Modules/PhotoacousticsLib/include/mitkPASpectralUnmixingFilterBase.h +++ b/Modules/PhotoacousticsLib/include/mitkPASpectralUnmixingFilterBase.h @@ -1,65 +1,71 @@ /*=================================================================== 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); void AddChromophore(mitk::pa::PropertyCalculator::ChromophoreType); void AddWavelength(int wavelength); + void AddWeight(int Weight); + ofstream myfile; + protected: SpectralUnmixingFilterBase(); virtual ~SpectralUnmixingFilterBase(); std::vector m_Chromophore; std::vector m_Wavelength; virtual Eigen::VectorXf SpectralUnmixingAlgorithm(Eigen::Matrix EndmemberMatrix, Eigen::VectorXf inputVector) = 0; + std::vector m_Weight; + + private: virtual void CheckPreConditions(unsigned int size, unsigned int NumberOfInputImages, const float* inputDataArray); virtual void GenerateData() override; virtual void InitializeOutputs(unsigned int TotalNumberOfSequences); 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 // MITKPHOTOACOUSTICSPECTRALUNMIXINGFILTERBASE_ diff --git a/Modules/PhotoacousticsLib/include/mitkPASpectralUnmixingFilterVigra.h b/Modules/PhotoacousticsLib/include/mitkPASpectralUnmixingFilterVigra.h index 54b2f65fd0..7ecbd0c27a 100644 --- a/Modules/PhotoacousticsLib/include/mitkPASpectralUnmixingFilterVigra.h +++ b/Modules/PhotoacousticsLib/include/mitkPASpectralUnmixingFilterVigra.h @@ -1,57 +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, LS, vigratest }; void mitk::pa::SpectralUnmixingFilterVigra::SetAlgorithm(VigraAlgortihmType SetAlgorithmIndex); - protected: SpectralUnmixingFilterVigra(); virtual ~SpectralUnmixingFilterVigra(); 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 52421e2820..127186fba1 100644 --- a/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingFilterBase.cpp +++ b/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingFilterBase.cpp @@ -1,227 +1,239 @@ /*=================================================================== 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(int Weight) +{ + m_Weight.push_back(Weight); + MITK_INFO << "Weight: " << m_Weight[m_Weight.size()-1]; +} + 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 NumberOfInputImages = input->GetDimensions()[2]; unsigned int size = xDim * yDim * NumberOfInputImages; MITK_INFO << "x dimension: " << xDim; MITK_INFO << "y dimension: " << yDim; MITK_INFO << "z dimension: " << NumberOfInputImages; unsigned int sequenceSize = m_Wavelength.size(); unsigned int TotalNumberOfSequences = NumberOfInputImages / sequenceSize; MITK_INFO << "TotalNumberOfSequences: " << TotalNumberOfSequences; InitializeOutputs(TotalNumberOfSequences); 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, NumberOfInputImages, inputDataArray); // test to see pixel values @ txt file //ofstream myfile; //myfile.open("PASpectralUnmixingPixelValues.txt"); - std::chrono::steady_clock::time_point _start(std::chrono::steady_clock::now()); + + myfile.open("SimplexNormalisation.txt"); + for (unsigned int SequenceCounter = 0; SequenceCounter < TotalNumberOfSequences;++SequenceCounter) { + MITK_INFO << "SequenceCounter: " << SequenceCounter; //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(sequenceSize); for (unsigned int z = 0; z < sequenceSize; z++) { // Get pixel value of pixel x,y @ wavelength z unsigned int pixelNumber = (xDim*yDim*(z+SequenceCounter*sequenceSize)) + x * yDim + y; auto pixel = inputDataArray[pixelNumber]; //write all wavelength absorbtion values for one(!) pixel to a vector inputVector[z] = pixel; } Eigen::VectorXf resultVector = SpectralUnmixingAlgorithm(EndmemberMatrix, inputVector); // write output 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]; } //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"; } } } 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 << "Chronotets: " << std::chrono::duration_cast>(_test - _start).count() << "s"; + + MITK_INFO << "Chrono: " << std::chrono::duration_cast>(_end - _start).count() << "s"; MITK_INFO << "GENERATING DATA...[DONE]"; - //myfile.close(); - MITK_INFO << "Chrono: " << std::chrono::duration_cast>( _end - _start).count() << "s"; + 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) MITK_WARN << "NUMBER OF WAVELENGTHS < NUMBER OF INPUT IMAGES"; if (m_Wavelength.size() > NumberOfInputImages) mitkThrow() << "ERROR! REMOVE WAVELENGTHS!"; // 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 TotalNumberOfSequences) { unsigned int numberOfInputs = GetNumberOfIndexedInputs(); unsigned int numberOfOutputs = GetNumberOfIndexedOutputs(); MITK_INFO << "Inputs: " << numberOfInputs << " Outputs: " << numberOfOutputs; // Set dimensions (3) and pixel type (float) for output 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; // 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); } } 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::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/mitkPASpectralUnmixingFilterSimplex.cpp b/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingFilterSimplex.cpp index 5283c3ef1f..7fb91f3159 100644 --- a/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingFilterSimplex.cpp +++ b/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingFilterSimplex.cpp @@ -1,169 +1,176 @@ /*=================================================================== 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 // Includes for AddEnmemberMatrix #include #include mitk::pa::SpectralUnmixingFilterSimplex::SpectralUnmixingFilterSimplex() { } mitk::pa::SpectralUnmixingFilterSimplex::~SpectralUnmixingFilterSimplex() { } Eigen::VectorXf mitk::pa::SpectralUnmixingFilterSimplex::SpectralUnmixingAlgorithm( Eigen::Matrix EndmemberMatrix, Eigen::VectorXf inputVector) { int numberOfChromophores = EndmemberMatrix.cols(); Eigen::VectorXf resultVector(numberOfChromophores); Eigen::VectorXf normalizedInputVector(EndmemberMatrix.rows()); normalizedInputVector = Normalization(EndmemberMatrix, inputVector); + //normalizedInputVector = inputVector; + float VolumeMax = simplexVolume(EndmemberMatrix); for (int i = 0; i < numberOfChromophores; ++i) { Eigen::Matrix A = GenerateA(EndmemberMatrix, normalizedInputVector, i); float Volume = simplexVolume(A); - //MITK_INFO << "VolumeMax: " << VolumeMax; - //MITK_INFO << "Volume: " << Volume; - //MITK_INFO << "Result vector[i]: " << Volume / VolumeMax; + + resultVector[i] = Volume / VolumeMax; + + myfile << "resultVector["< mitk::pa::SpectralUnmixingFilterSimplex::GenerateA (Eigen::Matrix EndmemberMatrix, Eigen::VectorXf inputVector, int i) { Eigen::Matrix A = EndmemberMatrix; int numberOfChromophores = EndmemberMatrix.cols(); for (int j = 0; j < numberOfChromophores; ++j) { A(i, j) = inputVector(j); } return A; } Eigen::Matrix mitk::pa::SpectralUnmixingFilterSimplex::GenerateD2 (Eigen::Matrix A) { int numberOfChromophores = A.cols(); Eigen::Matrix D2(numberOfChromophores, numberOfChromophores); for (int i = 0; i < numberOfChromophores; ++i) { for (int j = 0; j < numberOfChromophores; ++j) { Eigen::VectorXf x = A.col(i) - A.col(j); //MITK_INFO << "a_col_i: " < Matrix) { float Volume; int numberOfChromophores = Matrix.cols(); Eigen::Matrix C(numberOfChromophores + 1, numberOfChromophores + 1); Eigen::Matrix D2 = GenerateD2(Matrix); for (int i = 0; i < numberOfChromophores; ++i) { for (int j = 0; j < numberOfChromophores; ++j) { C(i, j) = D2(i, j); } C(i, numberOfChromophores) = 1; for (int k = 0; k < numberOfChromophores; ++k) { C(numberOfChromophores, k) = 1; } C(numberOfChromophores, numberOfChromophores) = 0; } float detC = -C.determinant();// determinate von C float denominator = (factorial(numberOfChromophores - 1)) ^ 2 * 2 ^ (numberOfChromophores - 1)*(-1) ^ numberOfChromophores; Volume = std::sqrt(detC / denominator); //MITK_INFO << "detC: " << detC; //MITK_INFO << "denominator: " << denominator; //MITK_INFO << "Volume: " << Volume; return Volume; } int mitk::pa::SpectralUnmixingFilterSimplex::factorial(int n) { if (n == 1) return 1; else return n * factorial(n - 1); } Eigen::VectorXf mitk::pa::SpectralUnmixingFilterSimplex::Normalization( Eigen::Matrix EndmemberMatrix, Eigen::VectorXf inputVector) { int numberOfWavelengths = inputVector.rows(); Eigen::VectorXf result(numberOfWavelengths); float NormalizationFactor = 1; float foo; float norm = 0; for (int i = 0; i < numberOfWavelengths; ++i) { foo = EndmemberMatrix(i, 0) - EndmemberMatrix(i, 1); if (std::abs(foo) > norm) norm = std::abs(foo); } -ofstream myfile; -myfile.open("SimplexNormalisation.txt"); +//ofstream myfile; +//myfile.open("SimplexNormalisation.txt"); + //NormalizationFactor = inputVector[0] * 2 / norm; + myfile << "Normalizationfactor " << NormalizationFactor << "\n"; for (int i = 0; i < numberOfWavelengths; ++i) { - NormalizationFactor = inputVector[i] * 2 / norm; - myfile << "Normalizationfactor " << NormalizationFactor << "\n"; + result[i]=(inputVector[i]/NormalizationFactor); } -myfile.close(); +//myfile.close(); return result; } diff --git a/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingFilterVigra.cpp b/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingFilterVigra.cpp index 345b97dd4c..287765b241 100644 --- a/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingFilterVigra.cpp +++ b/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingFilterVigra.cpp @@ -1,150 +1,157 @@ /*=================================================================== 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) { 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; - 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) 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 = { 1,1,1,1,1 }; //am besten über GUI Eingabe festlegen + 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) linearSolve(A, b, x); //mitkThrow() << "nothing implemented"; else mitkThrow() << "404 VIGRA ALGORITHM NOT FOUND"; for (int k = 0; k < numberOfChromophores; ++k) { float test = x(k, 0); resultVector[k] = test; } 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 a4aa3d27a9..a361b09c16 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,367 +1,378 @@ /*=================================================================== 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); } 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; if (Algorithm == "householderQr") { m_SpectralUnmixingFilter = mitk::pa::LinearSpectralUnmixingFilter::New(); dynamic_cast(m_SpectralUnmixingFilter.GetPointer()) ->SetAlgorithm(mitk::pa::LinearSpectralUnmixingFilter::AlgortihmType::householderQr); } else if (Algorithm == "ldlt") { m_SpectralUnmixingFilter = mitk::pa::LinearSpectralUnmixingFilter::New(); dynamic_cast(m_SpectralUnmixingFilter.GetPointer()) ->SetAlgorithm(mitk::pa::LinearSpectralUnmixingFilter::AlgortihmType::ldlt); } else if (Algorithm == "llt") { m_SpectralUnmixingFilter = mitk::pa::LinearSpectralUnmixingFilter::New(); dynamic_cast(m_SpectralUnmixingFilter.GetPointer()) ->SetAlgorithm(mitk::pa::LinearSpectralUnmixingFilter::AlgortihmType::llt); } else if (Algorithm == "colPivHouseholderQr") { m_SpectralUnmixingFilter = mitk::pa::LinearSpectralUnmixingFilter::New(); dynamic_cast(m_SpectralUnmixingFilter.GetPointer()) ->SetAlgorithm(mitk::pa::LinearSpectralUnmixingFilter::AlgortihmType::colPivHouseholderQr); } else if (Algorithm == "jacobiSvd") { m_SpectralUnmixingFilter = mitk::pa::LinearSpectralUnmixingFilter::New(); dynamic_cast(m_SpectralUnmixingFilter.GetPointer()) ->SetAlgorithm(mitk::pa::LinearSpectralUnmixingFilter::AlgortihmType::jacobiSvd); } else if (Algorithm == "fullPivLu") { m_SpectralUnmixingFilter = mitk::pa::LinearSpectralUnmixingFilter::New(); dynamic_cast(m_SpectralUnmixingFilter.GetPointer()) ->SetAlgorithm(mitk::pa::LinearSpectralUnmixingFilter::AlgortihmType::fullPivLu); } else if (Algorithm == "fullPivHouseholderQr") { m_SpectralUnmixingFilter = mitk::pa::LinearSpectralUnmixingFilter::New(); dynamic_cast(m_SpectralUnmixingFilter.GetPointer()) ->SetAlgorithm(mitk::pa::LinearSpectralUnmixingFilter::AlgortihmType::fullPivHouseholderQr); } else if (Algorithm == "test") { m_SpectralUnmixingFilter = mitk::pa::LinearSpectralUnmixingFilter::New(); dynamic_cast(m_SpectralUnmixingFilter.GetPointer()) ->SetAlgorithm(mitk::pa::LinearSpectralUnmixingFilter::AlgortihmType::test); } else if (Algorithm == "NNLARS") { m_SpectralUnmixingFilter = mitk::pa::SpectralUnmixingFilterVigra::New(); dynamic_cast(m_SpectralUnmixingFilter.GetPointer()) ->SetAlgorithm(mitk::pa::SpectralUnmixingFilterVigra::VigraAlgortihmType::LARS); } else if (Algorithm == "NNGoldfarb") { m_SpectralUnmixingFilter = mitk::pa::SpectralUnmixingFilterVigra::New(); dynamic_cast(m_SpectralUnmixingFilter.GetPointer()) ->SetAlgorithm(mitk::pa::SpectralUnmixingFilterVigra::VigraAlgortihmType::GOLDFARB); } else if (Algorithm == "weighted") { m_SpectralUnmixingFilter = mitk::pa::SpectralUnmixingFilterVigra::New(); dynamic_cast(m_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) + m_SpectralUnmixingFilter->AddWeight(Weight); + ++colunm; + } } else if (Algorithm == "LS") { m_SpectralUnmixingFilter = mitk::pa::SpectralUnmixingFilterVigra::New(); dynamic_cast(m_SpectralUnmixingFilter.GetPointer()) ->SetAlgorithm(mitk::pa::SpectralUnmixingFilterVigra::VigraAlgortihmType::LS); } else if (Algorithm == "vigratest") { m_SpectralUnmixingFilter = mitk::pa::SpectralUnmixingFilterVigra::New(); dynamic_cast(m_SpectralUnmixingFilter.GetPointer()) ->SetAlgorithm(mitk::pa::SpectralUnmixingFilterVigra::VigraAlgortihmType::vigratest); } else if (Algorithm == "SimplexMax") { m_SpectralUnmixingFilter = mitk::pa::SpectralUnmixingFilterSimplex::New(); dynamic_cast(m_SpectralUnmixingFilter.GetPointer()); } else mitkThrow() << "404 ALGORITHM NOT FOUND!"; 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); outputCounter += 1; 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); outputCounter += 1; 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); outputCounter += 1; 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"); this->GetDataStorage()->Add(dataNodesO2); MITK_INFO << "[DONE]"; } mitk::RenderingManager::GetInstance()->InitializeViewsByBoundingObjects(this->GetDataStorage()); MITK_INFO << "Adding images to DataStorage...[DONE]"; } } } diff --git a/Plugins/org.mitk.gui.qt.photoacoustics.spectralunmixing/src/internal/SpectralUnmixingControls.ui b/Plugins/org.mitk.gui.qt.photoacoustics.spectralunmixing/src/internal/SpectralUnmixingControls.ui index 8e8067b944..677fb61c98 100644 --- a/Plugins/org.mitk.gui.qt.photoacoustics.spectralunmixing/src/internal/SpectralUnmixingControls.ui +++ b/Plugins/org.mitk.gui.qt.photoacoustics.spectralunmixing/src/internal/SpectralUnmixingControls.ui @@ -1,847 +1,872 @@ SpectralUnmixingControls 0 0 398 575 0 0 QmitkTemplate QLabel { color: rgb(255, 0, 0) } Please select an image! Qt::Horizontal Do image processing Perform spectral unmixing Qt::Horizontal 75 true Wavelengths settings 0 75 16777215 75 10 true Qt::SolidLine 42 50 30 30 λ/nm Weights 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 720 755 831 906 942 + + + 60 + + + + + 85 + + + + + 90 + + + + + 69 + + + + + 52 + + Qt::Horizontal 75 true Chromophore settings Oxyhemoglobin true true Deoxygenated hemoglobin true true false Melanin Additional Endmember Qt::Horizontal 75 true Unmixing algorithm true MS Shell Dlg 2 false 21 2147483647 ==CHOSE ALGORITHM== ==QR decomposition== householderQr colPivHouseholderQr fullPivHouseholderQr ==LU decompositon fullPivLu ==Cholesky decompostion== ldlt llt .. ==Least squares== LS jacobiSvd NNLARS NNGoldfarb weighted ==Others== SimplexMax Lagrange ==Test algorithms== test vigratest Qt::Horizontal 75 true Oxygen saturation calculate sO2 false 380 82 Value Min Hb Min HbO2 Min Sum Min SO2 % 100 200 200 50 Qt::Horizontal Qt::Vertical 20 69