diff --git a/Modules/PhotoacousticsLib/include/mitkPASpectralUnmixingFilter.h b/Modules/PhotoacousticsLib/include/mitkPASpectralUnmixingFilter.h index e29763febd..5e79e1d997 100644 --- a/Modules/PhotoacousticsLib/include/mitkPASpectralUnmixingFilter.h +++ b/Modules/PhotoacousticsLib/include/mitkPASpectralUnmixingFilter.h @@ -1,78 +1,78 @@ /*=================================================================== 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 MITKPHOTOACOUSTICSPECTRALUNMIXINGFILTER_H #define MITKPHOTOACOUSTICSPECTRALUNMIXINGFILTER_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 SpectralUnmixingFilter : public mitk::ImageToImageFilter { public: mitkClassMacro(SpectralUnmixingFilter, mitk::ImageToImageFilter) itkFactorylessNewMacro(Self) // Contains all (so far) possible chromophores enum ChromophoreType { OXYGENATED_HEMOGLOBIN = 1, DEOXYGENATED_HEMOGLOBIN = 2 }; // Void to creat m_vector of all chosen chromophores with push back method void AddChromophore(ChromophoreType); std::vector m_Chromophore; // Void to creat m_vector of all wavelengths with push back method void AddWavelength(int wavelength); std::vector m_Wavelength; - + protected: SpectralUnmixingFilter(); virtual ~SpectralUnmixingFilter(); private: // Void checking precondtions possibly throwing exeptions - virtual void CheckPreConditions(unsigned int NumberOfInputImages); + virtual void CheckPreConditions(unsigned int size, unsigned int NumberOfInputImages, const float* inputDataArray); virtual void GenerateData(); virtual void InitializeOutputs(); // Void to creat Eigen::Matrix of all absorbtions // @ specific wavelength (columns) of chromophores (rows) Eigen::Matrix AddEndmemberMatrix(); Eigen::Matrix EndmemberMatrix; PropertyCalculator::Pointer m_PropertyCalculator; // Test algorithm for SU --> later a new class should be set up Eigen::VectorXf SpectralUnmixingTestAlgorithm(Eigen::Matrix EndmemberMatrix, Eigen::VectorXf inputVector); }; } } #endif // MITKPHOTOACOUSTICSPECTRALUNMIXINGFILTER_H diff --git a/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingFilter.cpp b/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingFilter.cpp index 2505c5ec09..e3907ac790 100644 --- a/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingFilter.cpp +++ b/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingFilter.cpp @@ -1,201 +1,245 @@ /*=================================================================== 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 "mitkPASpectralUnmixingFilter.h" // Includes for AddEnmemberMatrix #include "mitkPAPropertyCalculator.h" #include // ImageAccessor #include #include // For testing reasons #include mitk::pa::SpectralUnmixingFilter::SpectralUnmixingFilter() { this->SetNumberOfIndexedOutputs(2); for (unsigned int i = 0; iSetNthOutput(i, mitk::Image::New()); } m_PropertyCalculator = mitk::pa::PropertyCalculator::New(); } mitk::pa::SpectralUnmixingFilter::~SpectralUnmixingFilter() { } void mitk::pa::SpectralUnmixingFilter::AddWavelength(int wavelength) { m_Wavelength.push_back(wavelength); } void mitk::pa::SpectralUnmixingFilter::AddChromophore(ChromophoreType chromophore) { m_Chromophore.push_back(chromophore); + } void mitk::pa::SpectralUnmixingFilter::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; - CheckPreConditions(zDim); InitializeOutputs(); + + EndmemberMatrix = AddEndmemberMatrix(); // 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 = SpectralUnmixingTestAlgorithm(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[x*yDim + y] = resultVector[outputIdx]; } + myfile << "Input Pixel(x,y): " << x << "," << y << "\n" << inputVector << "\n"; + myfile << "Result: " << "\n Hb: " << resultVector[0] << "\n HbO2: " << resultVector[1] <<"\n"; } } MITK_INFO << "GENERATING DATA...[DONE]"; + myfile.close(); } -void mitk::pa::SpectralUnmixingFilter::CheckPreConditions(unsigned int NumberOfInputImages) +void mitk::pa::SpectralUnmixingFilter::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::SpectralUnmixingFilter::InitializeOutputs() { unsigned int numberOfInputs = GetNumberOfIndexedInputs(); unsigned int numberOfOutputs = GetNumberOfIndexedOutputs(); - MITK_INFO << "Inputs: " << numberOfInputs << " Outputs: " << numberOfOutputs; + //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 rows and #wavelengths columns //so Matrix Element (i,j) contains the Absorbtion of chromophore i @ wavelength j Eigen::Matrix mitk::pa::SpectralUnmixingFilter::AddEndmemberMatrix() { unsigned int numberOfChromophores = m_Chromophore.size(); //rows unsigned int numberOfWavelengths = m_Wavelength.size(); //columns Eigen::Matrix EndmemberMatrix(numberOfChromophores, numberOfWavelengths); //loop over i rows (Chromophores) for(unsigned int i = 0; i < numberOfChromophores; ++i) { //loop over j columns (Wavelengths) for (unsigned int j = 0; j < numberOfWavelengths; ++j) { //writes @ Matrix element (i,j) the absorbtion wavelength of the propertycalculator.cpp EndmemberMatrix(i,j)= m_PropertyCalculator->GetAbsorptionForWavelength( static_cast(m_Chromophore[i]), m_Wavelength[j]); + auto testtype = m_PropertyCalculator->GetAbsorptionForWavelength( + static_cast(m_Chromophore[i]), m_Wavelength[j]); /* Test to see what gets written in the Matrix: + MITK_INFO << "TEST_TYPE Matrix: " << typeid(EndmemberMatrix(i,j)).name(); MITK_INFO << "map type: " << static_cast(m_Chromophore[i]); MITK_INFO << "wavelength: " << m_Wavelength[j]; MITK_INFO << "Matrixelement: (" << i << ", " << j << ") Absorbtion: " << EndmemberMatrix(i, j);/**/ if (EndmemberMatrix(i, j) == 0) { m_Wavelength.clear(); mitkThrow() << "WAVELENGTH "<< m_Wavelength[j] << "nm NOT SUPPORTED!"; } } } MITK_INFO << "GENERATING ENMEMBERMATRIX [DONE]!"; return EndmemberMatrix; } Eigen::VectorXf mitk::pa::SpectralUnmixingFilter::SpectralUnmixingTestAlgorithm( Eigen::Matrix EndmemberMatrix, Eigen::VectorXf inputVector) { - //test "solver" (input = output) - Eigen::VectorXf resultVector = inputVector; - return resultVector; - - //llt solver - //Eigen::VectorXf resultVector = EndmemberMatrix.llt().solve(inputVector); - //return resultVector; + int treshold = 0; + int defaultValue = -1; - //householderqr solver - //Eigen::VectorXf resultVector =EndmemberMatrix.householderQr().solve(inputVector); - //return resultVector; + //*llt solver + //Eigen::Vector2f resultVector = EndmemberMatrix.llt().solve(inputVector); + //double relative_error = (EndmemberMatrix*inputVector - lltVector).norm() / lltVector.norm(); // norm() is L2 norm + /**/ + + //*householderqr solver + Eigen::Vector2f resultVector = EndmemberMatrix.householderQr().solve(inputVector);/**/ + + + /*Set threshold and replace with default value if under threshold + for (int i = 0; i < 2; ++i) + { + if (resultVector[i] < treshold) + { + resultVector[i] = defaultValue; + MITK_INFO << "UNMIXING RESULT N/A"; + } + }/**/ //test other solvers https://eigen.tuxfamily.org/dox/group__TutorialLinearAlgebra.html - // define treshold for relativ error and set value to eg. -1 ;) - /*double relative_error = (EndmemberMatrix*inputVector - resultVector).norm() / resultVector.norm(); // norm() is L2 norm - MITK_INFO << relative_error;/**/ + + 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 7235869b08..417f2b2bed 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,200 +1,188 @@ /*=================================================================== 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 -/*//test mitk image -#include -#include -#include */ - // Include to perform Spectral Unmixing #include "mitkPASpectralUnmixingFilter.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"; } } -// Checking which chromophores wanted for SU if none throw exeption! -void SpectralUnmixing::numberOfChromophores() -{ - auto m_SpectralUnmixingFilter = mitk::pa::SpectralUnmixingFilter::New(); - - 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::SpectralUnmixingFilter::ChromophoreType::OXYGENATED_HEMOGLOBIN); - } - if (DeOxbool) - { - numberofChromophores += 1; - MITK_INFO << "- Deoxygenated hemoglobin"; - // Set chromophore Deoxygenated hemoglobin: - m_SpectralUnmixingFilter->AddChromophore( - mitk::pa::SpectralUnmixingFilter::ChromophoreType::DEOXYGENATED_HEMOGLOBIN); - } - if (numberofChromophores == 0) - { - mitkThrow() << "PRESS 'IGNORE' AND CHOOSE A CHROMOPHORE!"; - } -} - 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(); - MITK_INFO << "GENERATING DATA..."; - auto m_SpectralUnmixingFilter = mitk::pa::SpectralUnmixingFilter::New(); m_SpectralUnmixingFilter->SetInput(image); // Wavelength implementation into fiter 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::SpectralUnmixingFilter::ChromophoreType::OXYGENATED_HEMOGLOBIN); + } + if (DeOxbool) + { + numberofChromophores += 1; + MITK_INFO << "- Deoxygenated hemoglobin"; + // Set chromophore Deoxygenated hemoglobin: + m_SpectralUnmixingFilter->AddChromophore( + mitk::pa::SpectralUnmixingFilter::ChromophoreType::DEOXYGENATED_HEMOGLOBIN); + } + if (numberofChromophores == 0) + { + mitkThrow() << "PRESS 'IGNORE' AND CHOOSE A CHROMOPHORE!"; + } MITK_INFO << "Updating Filter..."; m_SpectralUnmixingFilter->Update(); - - mitk::Image::Pointer HbO2 = m_SpectralUnmixingFilter->GetOutput(0); - mitk::DataNode::Pointer dataNodeHbO2 = mitk::DataNode::New(); - dataNodeHbO2->SetData(HbO2); - dataNodeHbO2->SetName("HbO2"); - this->GetDataStorage()->Add(dataNodeHbO2); - - mitk::Image::Pointer Hb = m_SpectralUnmixingFilter->GetOutput(1); + + mitk::Image::Pointer Hb = m_SpectralUnmixingFilter->GetOutput(0); mitk::DataNode::Pointer dataNodeHb = mitk::DataNode::New(); dataNodeHb->SetData(Hb); dataNodeHb->SetName("Hb"); this->GetDataStorage()->Add(dataNodeHb); - + + mitk::Image::Pointer HbO2 = m_SpectralUnmixingFilter->GetOutput(1); + mitk::DataNode::Pointer dataNodeHbO2 = mitk::DataNode::New(); + dataNodeHbO2->SetData(HbO2); + dataNodeHbO2->SetName("HbO2"); + this->GetDataStorage()->Add(dataNodeHbO2); + 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/SpectralUnmixing.h b/Plugins/org.mitk.gui.qt.photoacoustics.spectralunmixing/src/internal/SpectralUnmixing.h index 9b09bfd53f..8236ff39b8 100644 --- a/Plugins/org.mitk.gui.qt.photoacoustics.spectralunmixing/src/internal/SpectralUnmixing.h +++ b/Plugins/org.mitk.gui.qt.photoacoustics.spectralunmixing/src/internal/SpectralUnmixing.h @@ -1,64 +1,63 @@ /*=================================================================== 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 SpectralUnmixing_h #define SpectralUnmixing_h #include #include #include "ui_SpectralUnmixingControls.h" #include "mitkPAPropertyCalculator.h" class SpectralUnmixing : public QmitkAbstractView { // this is needed for all Qt objects that should have a Qt meta-object // (everything that derives from QObject and wants to have signal/slots) Q_OBJECT public: static const std::string VIEW_ID; protected: virtual void CreateQtPartControl(QWidget *parent) override; virtual void SetFocus() override; /// \brief called by QmitkFunctionality when DataManager's selection has changed virtual void OnSelectionChanged(berry::IWorkbenchPart::Pointer source, const QList &nodes) override; /// \brief Called when the user clicks the GUI button void DoImageProcessing(); virtual void Wavelength(); virtual void ClearWavelength(); - virtual void numberOfChromophores(); Ui::SpectralUnmixingControls m_Controls; mitk::pa::PropertyCalculator::Pointer m_PropertyCalculator; // Add Wavelengths with button: int size = 0; std::vector m_Wavelengths; // Selection of Chromophores bool DeOxbool; bool Oxbool; }; #endif // SpectralUnmixing_h