diff --git a/Modules/PhotoacousticsLib/include/mitkPASpectralUnmixingFilter.h b/Modules/PhotoacousticsLib/include/mitkPASpectralUnmixingFilter.h index ebda78518c..2b3b915edb 100644 --- a/Modules/PhotoacousticsLib/include/mitkPASpectralUnmixingFilter.h +++ b/Modules/PhotoacousticsLib/include/mitkPASpectralUnmixingFilter.h @@ -1,94 +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) 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); - + Eigen::VectorXd SpectralUnmixingAlgorithms(Eigen::VectorXd inputVector); }; } } #endif // MITKPHOTOACOUSTICSPECTRALUNMIXINGFILTER_H diff --git a/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingFilter.cpp b/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingFilter.cpp index 1f7531a493..ca86333880 100644 --- a/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingFilter.cpp +++ b/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingFilter.cpp @@ -1,220 +1,207 @@ /*=================================================================== 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() { MITK_INFO << "GENERATING DATA.."; - //checking Preconditions: - unsigned int NumberOfInputImages = m_Dimensions[2]; - MITK_INFO << "NumberOfInputImages: " << NumberOfInputImages; - CheckPreConditions(NumberOfInputImages); + unsigned int xDim = m_Dimensions[0]; + unsigned int yDim = m_Dimensions[1]; + unsigned int zDim = m_Dimensions[2]; - InitializeOutputs(); - - EndmemberMatrix = AddEndmemberMatrix(); - - MITK_INFO << "GENERATING DATA...[DONE]"; + MITK_INFO << "NumberOfInputImages: " << zDim; - // code recreaction from "old" SUF.cpp see end of document - + CheckPreConditions(zDim); - - //*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++{ - - - /* - 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); + InitializeOutputs(); mitk::Image::Pointer data = GetInput(); - for (unsigned int x = 0; x < xdim; x++) + for (unsigned int x = 0; x < xDim; x++) { - for (unsigned int y = 0; y < ydim; y++) + for (unsigned int y = 0; y < yDim; y++) { - Eigen::VectorXd b(numberOfInputs); + Eigen::VectorXd inputVector(numberOfInputs); - for (unsigned int z = 0; z < zdim; z++) + for (unsigned int z = 0; z < zDim; z++) { + //inputVector[z] = wertvonpixelmitderwellenlänge; - b(z) = ; } - Eigen::Vector3d reultVector = EndmemberMatrix.householderQr().solve(b); + Eigen::VectorXd resultVector = SpectralUnmixingAlgorithms(inputVector); - for (int outpuIdx = 0; outpuIdx < GetNumberOfIndexedOutputs(); ++outpuIdx) + for (int outputIdx = 0; outputIdx < GetNumberOfIndexedOutputs(); ++outputIdx) { - auto output = GetOutput(outpuIdx); + auto output = GetOutput(outputIdx); mitk::ImageWriteAccessor writeOutput(output); float* writeBuffer = (float *)writeOutput.GetData(); - writeBuffer[y*xdim + x] = reultVector[outpuIdx]; + writeBuffer[y*xDim + x] = resultVector[outputIdx]; } } } - - - - /*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++}*/ + MITK_INFO << "GENERATING DATA...[DONE]"; } // 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]"; } // 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]); /* 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!"; return EndmemberMatrix; } +// Perform SU algorithm +Eigen::VectorXd mitk::pa::SpectralUnmixingFilter::SpectralUnmixingAlgorithms(Eigen::VectorXd inputVector) +{ + EndmemberMatrix = AddEndmemberMatrix(); + Eigen::VectorXd resultVector = EndmemberMatrix.householderQr().solve(inputVector); + return resultVector; +} + /* +++++++++++++++++++++++++++++++++++++++++ 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]; } } /**/