diff --git a/Modules/PhotoacousticsLib/include/mitkPASpectralUnmixingFilter.h b/Modules/PhotoacousticsLib/include/mitkPASpectralUnmixingFilter.h index de994dd393..ebda78518c 100644 --- a/Modules/PhotoacousticsLib/include/mitkPASpectralUnmixingFilter.h +++ b/Modules/PhotoacousticsLib/include/mitkPASpectralUnmixingFilter.h @@ -1,76 +1,94 @@ /*=================================================================== 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" //Includes for smart pointer usage #include "mitkCommon.h" #include "itkLightObject.h" #include #include "mitkPAPropertyCalculator.h" #include +// Test with ImagePixelAccessor +#include +#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 = 1, DEOXYGENATED = 2 }; //Void to creat m_vector of all chosen chromophores with push back method void SetChromophores(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; + void SetDimensions(int dimension); + std::vector m_Dimensions; + + //Void to creat Eigen::Matrix of all absorbtions //@ specific wavelength (columns) of chromophores (rows) - virtual void AddEndmemberMatrix(); - + Eigen::Matrix AddEndmemberMatrix(); + + protected: SpectralUnmixingFilter(); virtual ~SpectralUnmixingFilter(); private: virtual void GenerateData(); PropertyCalculator::Pointer m_PropertyCalculator; unsigned int numberofchromophores; unsigned int numberofwavelengths; unsigned int numberOfInputs; unsigned int numberOfOutputs; long length; Eigen::Matrix EndmemberMatrix; + + + //'New': + virtual void CheckPreConditions(unsigned int NumberOfInputImages); + + virtual void InitializeOutputs(); + //Eigen::VectorXd OutputVector = (Eigen::VectorXd, Eigen::Matrix); + + }; } } #endif // MITKPHOTOACOUSTICSPECTRALUNMIXINGFILTER_H diff --git a/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingFilter.cpp b/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingFilter.cpp index 02c382298f..1f7531a493 100644 --- a/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingFilter.cpp +++ b/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingFilter.cpp @@ -1,121 +1,220 @@ /*=================================================================== 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" #include #include // Includes for AddEnmemberMatrix #include "mitkPAPropertyCalculator.h" #include +// Test with ImagePixelAccessor +#include +#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::SetChromophores(ChromophoreType chromophore) { m_Chromophore.push_back(chromophore); } void mitk::pa::SpectralUnmixingFilter::GenerateData() { - //code recreaction from "old" SUF.cpp + MITK_INFO << "GENERATING DATA.."; - // now creats to identical output images like the input for hb and hb02 + //checking Preconditions: + unsigned int NumberOfInputImages = m_Dimensions[2]; + MITK_INFO << "NumberOfInputImages: " << NumberOfInputImages; + CheckPreConditions(NumberOfInputImages); - MITK_INFO << "GENERATING DATA.."; - numberOfInputs = GetNumberOfIndexedInputs(); - numberOfOutputs = GetNumberOfIndexedOutputs(); + InitializeOutputs(); + + EndmemberMatrix = AddEndmemberMatrix(); + + MITK_INFO << "GENERATING DATA...[DONE]"; + + // code recreaction from "old" SUF.cpp see end of document + + + + //*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++{ - for (unsigned int outputIdx = 0; outputIdx < numberOfOutputs; outputIdx++) - { - GetOutput(outputIdx)->Initialize(GetInput(0)); - } - length = GetOutput(0)->GetDimension(0)*GetOutput(0)->GetDimension(1)*GetOutput(0)->GetDimension(2); - - for (int i = 0; i < length; i++) + /* + MITK_INFO << "dim0: " << GetOutput(0)->GetDimension(0) //128 + << " dim1: " << GetOutput(0)->GetDimension(1) //xxxx for header.nnrd: 3919 + << " dim2: " << GetOutput(0)->GetDimension(2);//1 + Dimensions of picture + + int xdim = GetOutput(0)->GetDimension(0); + int ydim = GetOutput(0)->GetDimension(1); + int zdim = GetOutput(0)->GetDimension(2); + + mitk::Image::Pointer data = GetInput(); + + for (unsigned int x = 0; x < xdim; x++) { - for (unsigned int j = 0; j < numberOfInputs; j++) + for (unsigned int y = 0; y < ydim; y++) { - mitk::Image::Pointer input = GetInput(j); - mitk::ImageReadAccessor readAccess(input, input->GetVolumeData()); - } - for (unsigned int outputIdx = 0; outputIdx < numberOfOutputs; outputIdx++) - { - mitk::Image::Pointer output = GetOutput(outputIdx); - mitk::ImageWriteAccessor writeOutput(output, output->GetVolumeData()); - double* outputArray = (double *)writeOutput.GetData(); + Eigen::VectorXd b(numberOfInputs); + + for (unsigned int z = 0; z < zdim; z++) + { + + b(z) = ; + + } + + Eigen::Vector3d reultVector = EndmemberMatrix.householderQr().solve(b); + + for (int outpuIdx = 0; outpuIdx < GetNumberOfIndexedOutputs(); ++outpuIdx) + { + auto output = GetOutput(outpuIdx); + mitk::ImageWriteAccessor writeOutput(output); + float* writeBuffer = (float *)writeOutput.GetData(); + + writeBuffer[y*xdim + x] = reultVector[outpuIdx]; + } + } } - MITK_INFO << "GENERATING DATA...[DONE]"; - - AddEndmemberMatrix(); + + + + /*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++}*/ +} + +// creats vector with x, y, z dimensions as entries +void mitk::pa::SpectralUnmixingFilter::SetDimensions(int Dimension) +{ + m_Dimensions.push_back(Dimension); } +// checks if number of Inputs == added wavelengths +void mitk::pa::SpectralUnmixingFilter::CheckPreConditions(unsigned int NumberOfInputImages) +{ + if (m_Wavelength.size() != NumberOfInputImages) + mitkThrow() << "CHECK INPUTS! WAVELENGTHERROR"; + MITK_INFO << "CHECK PRECONDITIONS ...[DONE]"; +} -//Void creats Matrix with #chromophores rows and #wavelengths columns -//so Matrix Element (i,j) contains the Absorbtion of chromophore j @ wavelength i -void mitk::pa::SpectralUnmixingFilter::AddEndmemberMatrix() +// Initialize Outputs +void mitk::pa::SpectralUnmixingFilter::InitializeOutputs() +{ + numberOfInputs = m_Dimensions[2]; + numberOfOutputs = GetNumberOfIndexedOutputs(); + for (unsigned int outputIdx = 0; outputIdx < numberOfOutputs; outputIdx++) + { + GetOutput(outputIdx)->Initialize(GetInput(0)); + } +} + +//Matrix with #chromophores rows and #wavelengths columns +//so Matrix Element (j,i) contains the Absorbtion of chromophore j @ wavelength i +Eigen::Matrix mitk::pa::SpectralUnmixingFilter::AddEndmemberMatrix() { numberofchromophores = m_Chromophore.size(); numberofwavelengths = m_Wavelength.size(); Eigen::Matrix EndmemberMatrix(numberofchromophores, numberofwavelengths); //loop over j rows (Chromophores) for(unsigned int j =0; j < numberofchromophores; ++j) { //loop over i columns (Wavelength) for (unsigned int i = 0; i < numberofwavelengths; ++i) { //writes @ Matrix element (i,j) the absorbtion wavelength of the propertycalculator.cpp EndmemberMatrix(j,i)= m_PropertyCalculator->GetAbsorptionForWavelength( static_cast(m_Chromophore[j]), m_Wavelength[i]); - if (EndmemberMatrix(i, j) == 0) - { - m_Wavelength.clear(); - MITK_INFO << "size: "<< m_Wavelength.size(); - mitkThrow() << "WAVELENGTH NOT SUPPORTED!"; - } /* Test to see what gets written in the Matrix: MITK_INFO << "map type: " << static_cast(m_Chromophore[j]); MITK_INFO << "wavelength: " << m_Wavelength[i]; MITK_INFO << "Matrixelement: (" << j << ", " << i << ") Absorbtion: " << EndmemberMatrix(j, i);/**/ + + if (EndmemberMatrix(j, i) == 0) + { + m_Wavelength.clear(); + mitkThrow() << "WAVELENGTH "<< m_Wavelength[i] << "nm NOT SUPPORTED!"; + } + } } - MITK_INFO << "GENERATING ENMEMBERMATRIX SUCCESSFUL!"; + //MITK_INFO << "GENERATING ENMEMBERMATRIX SUCCESSFUL!"; + return EndmemberMatrix; +} + +/* +++++++++++++++++++++++++++++++++++++++++ OLD CODE: +++++++++++++++++++++++++++++++++++++++++++++++++++ +was @ generate data + +length = GetOutput(0)->GetDimension(0)*GetOutput(0)->GetDimension(1)*GetOutput(0)->GetDimension(2); + +for (int i = 0; i < length; i++) +{ +Eigen::VectorXd b(numberOfInputs); +for (unsigned int j = 0; j < numberOfInputs; j++) +{ +mitk::Image::Pointer input = GetInput(j); +mitk::ImageReadAccessor readAccess(input, input->GetVolumeData()); +b(j) = ((const float*)readAccess.GetData())[i]; +} + +Eigen::Vector3d x = EndmemberMatrix.householderQr().solve(b); + + +if (i == 0) +{ +MITK_INFO << "for i=0 b(#Inputs)"; +MITK_INFO << b; +MITK_INFO << "EndmemberMatrix"; +MITK_INFO << EndmemberMatrix; +MITK_INFO << "x"; +MITK_INFO << x; +} + + +for (unsigned int outputIdx = 0; outputIdx < numberOfOutputs; outputIdx++) +{ +mitk::Image::Pointer output = GetOutput(outputIdx); +mitk::ImageWriteAccessor writeOutput(output, output->GetVolumeData()); +double* outputArray = (double *)writeOutput.GetData(); +outputArray[i] = x[outputIdx]; +} } +/**/ 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 f1f6f28f96..ffa0e37da8 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,243 +1,249 @@ /*=================================================================== 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); } // Add Wavelength is working, BUT in the Plugin! Not with same implementation at the filter // probably because the m_wavelengths vector is created new every time and not 'saved'. -// Alternativ as comment here and at the filter. void SpectralUnmixing::Wavelength() { if (m_Wavelengths.empty()) { size = 0; } wavelength = m_Controls.spinBoxAddWavelength->value(); m_Wavelengths.push_back(wavelength); MITK_INFO << "ADD WAVELENGTH: " << wavelength << "nm"; 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"; } } 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 (SOON)"; if (node->GetName(name)) { // a property called "name" was found for this DataNode message << "'" << name << "'"; } message << "."; MITK_INFO << message.str(); - //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ auto m_SpectralUnmixingFilter = mitk::pa::SpectralUnmixingFilter::New(); + for (unsigned int i = 0; i < 3; i++) + { + m_SpectralUnmixingFilter->SetDimensions(image->GetDimension(i)); + } + + // Checking which chromophores wanted for SU if none throw exeption! 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->SetChromophores( mitk::pa::SpectralUnmixingFilter::ChromophoreType(1)); } if (DeOxbool) { numberofChromophores += 1; MITK_INFO << "- Deoxygenated hemoglobin"; // Set chromophore Deoxygenated hemoglobin: m_SpectralUnmixingFilter->SetChromophores( mitk::pa::SpectralUnmixingFilter::ChromophoreType(2)); } if (numberofChromophores == 0) { mitkThrow() << "PRESS 'IGNORE' AND CHOOSE A CHROMOPHORE!"; } // Checking if number of wavelengths >= number of chromophores if (numberofChromophores > size) { mitkThrow() << "PRESS 'IGNORE' AND ADD MORE WAVELENGTHS!"; } - + //code recreaction from "old" SU.cpp MITK_INFO << "GENERATING DATA..."; unsigned int numberOfInputs = size; unsigned int numberOfOutputs = numberofChromophores; MITK_INFO << "Inputs: " << numberOfInputs << " Outputs: " << numberOfOutputs; mitk::DataStorage::Pointer storage = this->GetDataStorage(); mitk::DataNode::Pointer node = nodes.front(); mitk::Image* inputImage = dynamic_cast(data); const unsigned int dimensions[]{ inputImage->GetDimension(0), inputImage->GetDimension(1), 1 }; // just the first sequence for starters TODO GENERALIZE unsigned int dimension = 3; mitk::PixelType TPixel = mitk::MakeScalarPixelType(); - // another wavelength impementation for fiter needs "above" one + // additional wavelength implementation into fiter vector; needs "above" one for (unsigned int imageIndex = 0; imageIndex < numberOfInputs; imageIndex++) { mitk::Image::Pointer fooImage = mitk::Image::New(); fooImage->Initialize(TPixel, dimension, dimensions); wavelength = m_Wavelengths[imageIndex]; MITK_INFO << wavelength; m_SpectralUnmixingFilter->AddWavelength(wavelength); mitk::ImageReadAccessor inputAcc(inputImage, inputImage->GetSliceData(imageIndex)); fooImage->SetSlice(inputAcc.GetData(), 0); m_SpectralUnmixingFilter->SetInput(imageIndex, fooImage); } - - + MITK_INFO << "Updating Filter..."; m_SpectralUnmixingFilter->Update(); - + mitk::Image::Pointer HbO2 = m_SpectralUnmixingFilter->GetOutput(0); - + + //* HbO2->GetGeometry()->SetIndexToWorldTransform(inputImage->GetGeometry()->GetIndexToWorldTransform()); mitk::ImageWriteAccessor writeOutputHbO2(HbO2, HbO2->GetVolumeData()); - double* outputArrayHbO2 = (double *)writeOutputHbO2.GetData(); + double* outputArrayHbO2 = (double *)writeOutputHbO2.GetData();/**/ 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); + //* Hb->GetGeometry()->SetIndexToWorldTransform(inputImage->GetGeometry()->GetIndexToWorldTransform()); mitk::ImageWriteAccessor writeOutputHb(Hb, Hb->GetVolumeData()); - double* outputArrayHb = (double *)writeOutputHb.GetData(); + double* outputArrayHb = (double *)writeOutputHb.GetData();/**/ mitk::DataNode::Pointer dataNodeHb = mitk::DataNode::New(); dataNodeHb->SetData(Hb); dataNodeHb->SetName("Hb"); this->GetDataStorage()->Add(dataNodeHb); 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 987bdd90b0..73093ff3d6 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,123 +1,151 @@ SpectralUnmixingControls 0 0 222 399 0 0 QmitkTemplate - - + + + + true + - Oxyhemoglobin + Deoxygenated hemoglobin + + + true true + + false + - - - - Do image processing + + + + QLabel { color: rgb(255, 0, 0) } - Do Something + Please select an image! - + Qt::Vertical 20 40 - - - - QLabel { color: rgb(255, 0, 0) } + + + + Do image processing - Please select an image! + Perform spectral unmixing - - - - true - + + - Deoxygenated hemoglobin - - - true + Oxyhemoglobin true - - false - - + Add Wavelength [nm] + + + + Qt::Vertical + + + 2000 - + Clear all Wavelengths + + + + Qt::Horizontal + + + + + + + Qt::Horizontal + + + + + + + Qt::Horizontal + + +