diff --git a/Modules/PhotoacousticsLib/include/mitkPASpectralUnmixingFilterBase.h b/Modules/PhotoacousticsLib/include/mitkPASpectralUnmixingFilterBase.h index e9be406c52..076445a0e2 100644 --- a/Modules/PhotoacousticsLib/include/mitkPASpectralUnmixingFilterBase.h +++ b/Modules/PhotoacousticsLib/include/mitkPASpectralUnmixingFilterBase.h @@ -1,75 +1,165 @@ /*=================================================================== 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 { /** - * documentation - *\brief ,,,, + * \brief The spectral unmixing filter base is designed as superclass for several spectral unmixing filter eg. Eigen- or Vigrabased ones. + * One can add wavelengths and chromophores and get a unmixed output images out of one MITK input image using algorithms from the subclasses. + * + * Input: + * The unmixing input has to be a 3D MITK image where the XY-plane is a image and the Z-direction represents recordings for different + * wavelengths. Herein a XY-plane for a specific Z-direction will be called "image". Every image has to be assigned to a certain wavelength. + * The "AddWavelength" uses "push_back" to write float values into a vector. The first wavelength will correspond to the first image!!! + * If there a more input images 'I' then added wavelengths 'w' the filter base interprets the next x images as repetition sequence of the same + * wavelengths. If I % w !=0 the surplus image(s) will be dropped. + * Addtionaly one has to add chromophores from the property calculator class enum "ChromophoreType" with the "AddChromophore" method. + * This method as well uses "push_back" but the chosen (arbitary) order will be the order of the outputs. + * + * Output: + * The output will be one MITK image per chosen chromophore. Where the XY-plane is a image and the Z-direction represents recordings for different + * sequences. + * + * Subclasses: + * - mitkPASpectralUnmixingFilterVigra + * - mitkPALinearSpectralUnmixingFilter (uses Eigen algorithms) + * - mitkPASpectralUnmixingFilterSimplex + * + * Possible exceptions: + * - "PIXELTYPE ERROR! FLOAT 32 REQUIRED": The MITK input image has to consist out of floats. + * - "ERROR! REMOVE WAVELENGTHS!": One needs at least the same amount of images (z-dimension) then added wavelengths. + * - "ADD MORE WAVELENGTHS!": One needs at least the same amount of wavelengths then added chromophores. + * - "404 CHROMOPHORE NOT FOUND!": The chromophore is not a memeber of the enum "ChromophoreType" OR not supported by the filter. + * - "WAVELENGTH XXX nm NOT SUPPORTED!": The wavelength is not part of the proptery calculater data base. The data base can be found @ + * [...]\mitk\Modules\PhotoacousticsLib\Resources\spectralLIB.dat */ + class MITKPHOTOACOUSTICSLIB_EXPORT SpectralUnmixingFilterBase : public mitk::ImageToImageFilter { public: mitkClassMacro(SpectralUnmixingFilterBase, mitk::ImageToImageFilter); - void AddChromophore(mitk::pa::PropertyCalculator::ChromophoreType); + /** + * \brief AddChromophore takes mitk::pa::PropertyCalculator::ChromophoreType and writes them at the end of the m_Chromophore vector. + * The call of the method sets the order of the GetOutput method! + * + * @param chromophore has to be element of porperty calculater enum chromophore type + * @return for wavelength smaller then 300nm and larger then 1000nm the return will be 0, because not at the data base (2018/06/19) + */ + void AddChromophore(mitk::pa::PropertyCalculator::ChromophoreType chromophore); + + /** + * \brief AddWavelength takes integers and writes them at the end of the m_Wavelength vector. The first call of the method then + * corresponds to the first input image and so on. + * @wavelength database supports integers between 300 and 1000 nm + */ void AddWavelength(int wavelength); - void AddWeight(int Weight); - ofstream myfile; + void AddWeight(unsigned int weight); // has to be a part of the vigra class + ofstream myfile; // just for testing purposes; has to be removeed protected: + /** + * \brief Constructor creats proptery calculater smart pointer new() + */ SpectralUnmixingFilterBase(); virtual ~SpectralUnmixingFilterBase(); - std::vector m_Chromophore; - std::vector m_Wavelength; - - virtual Eigen::VectorXf SpectralUnmixingAlgorithm(Eigen::Matrix EndmemberMatrix, + /** + * \brief The subclasses will override the mehtod to calculate the spectral unmixing result vector. + * @param endmemberMatrix Matrix with number of chromophores colums and number of wavelengths rows so matrix element (i,j) contains + * the absorbtion of chromophore j @ wavelength i taken from the database by PropertyElement method. + * @param inputVector Vector containing values of one pixel of XY-plane image with number of wavelength rows (z-dimension of a sequenece) + * so the pixelvalue of the first wavelength is stored in inputVector[0] and so on. + */ + virtual Eigen::VectorXf SpectralUnmixingAlgorithm(Eigen::Matrix endmemberMatrix, Eigen::VectorXf inputVector) = 0; - std::vector m_Weight; + /* + * \brief Verbose gives more information to the console. Default value is false. + * @param m_Verbose is the boolian to activate the MITK_INFO logged to the console + */ + virtual void Verbose(bool verbose); + bool m_Verbose = false; + std::vector m_Chromophore; + std::vector m_Wavelength; + std::vector m_Weight; // has to be a part of the vigra class private: - virtual void CheckPreConditions(unsigned int size, unsigned int NumberOfInputImages, const float* inputDataArray); - + /* + * \brief Initialized output images with the same XY-plane size like the input image and total size in z-direction equals number of sequences. + * The pixel type is set to float. + * @param totalNumberOfSequences = (unsigned int) (numberOfInputImages / numberOfWavelength) >= 1 + */ + virtual void InitializeOutputs(unsigned int totalNumberOfSequences); + + /* + * \brief Checks if there are a suitable amount of wavelengths and if the input image consists of floats. + * @param numberOfInputImages corresponds to the z-dimension of the MITK input image + * @param inputDataArray contains every pixel of the MITK input image. Accessable like pseudo multidimensional array, that means + * pixel(x,y,z) = (xdim * yDim * z) + x * yDim + y where eg. "xdim" means the total size of the MITK input image in x direction. + * @throw if there are more wavelength then images + * @throw if there are more chromophores then wavelengths + * @throw if the pixel type is not float 32 + */ + virtual void CheckPreConditions(unsigned int numberOfInputImages, const float* inputDataArray); + + /* + * \brief Inherit from the "ImageToImageFilter" Superclass. Herain it calls InitializeOutputs, CalculateEndmemberMatrix and + * CheckPreConditions methods and enables pixelwise access to do spectral unmixing with the "SpectralUnmixingAlgorithm" + * method. In the end the method writes the results into the new MITK output images. + */ virtual void GenerateData() override; - virtual void InitializeOutputs(unsigned int TotalNumberOfSequences); + /* + * \brief Creats a Matrix with number of chromophores colums and number of wavelengths rows so matrix element (i,j) contains + * the absorbtion of chromophore j @ wavelength i. The absorbtion values are taken from the "PropertyElement" method. + * @param m_Chromophore is a vector of "PropertyCalculator::ChromophoreType" containing all selected chromophores for the unmixing + * @param m_Wavelength is a vector of integers containing all wavelengths of one sequence + */ virtual Eigen::Matrix CalculateEndmemberMatrix( std::vector m_Chromophore, std::vector m_Wavelength); - PropertyCalculator::Pointer m_PropertyCalculatorEigen; + /* + * \brief "PropertyElement" is the tool to access the absorbtion values out of the database using mitk::pa::PropertyCalculator::GetAbsorptionForWavelengt + * and checks if the requested wavelength is part of the database (not zero values). The "ONEENDMEMBER" is a pseudo absorber with static absorbtion 1 + * at every wavelength and is therefor not part of the database. If this one is the selected chromophore the return value is 1 for every wavelength. + * @param wavelength has to be integer between 300 and 1000 nm + * @param chromophore has to be a mitk::pa::PropertyCalculator::ChromophoreType + * @throw if mitk::pa::PropertyCalculator::GetAbsorptionForWavelengt returns 0, because this means that the delivered wavelength is not in the database. + */ virtual float PropertyElement(mitk::pa::PropertyCalculator::ChromophoreType, int wavelength); + + PropertyCalculator::Pointer m_PropertyCalculatorEigen; }; } } #endif // MITKPHOTOACOUSTICSPECTRALUNMIXINGFILTERBASE_ diff --git a/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingFilterBase.cpp b/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingFilterBase.cpp index d70b156fe9..e1e3e1d576 100644 --- a/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingFilterBase.cpp +++ b/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingFilterBase.cpp @@ -1,232 +1,206 @@ /*=================================================================== 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) +void mitk::pa::SpectralUnmixingFilterBase::AddWeight(unsigned int weight) { - m_Weight.push_back(Weight); + m_Weight.push_back(weight); +} + +void mitk::pa::SpectralUnmixingFilterBase::Verbose(bool verbose) +{ + m_Verbose = verbose; } void mitk::pa::SpectralUnmixingFilterBase::GenerateData() { - MITK_INFO << "GENERATING DATA.."; + MITK_INFO(m_Verbose) << "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; + unsigned int numberOfInputImages = input->GetDimensions()[2]; - MITK_INFO << "x dimension: " << xDim; - MITK_INFO << "y dimension: " << yDim; - MITK_INFO << "z dimension: " << NumberOfInputImages; + MITK_INFO(m_Verbose) << "x dimension: " << xDim; + MITK_INFO(m_Verbose) << "y dimension: " << yDim; + MITK_INFO(m_Verbose) << "z dimension: " << numberOfInputImages; unsigned int sequenceSize = m_Wavelength.size(); - unsigned int TotalNumberOfSequences = NumberOfInputImages / sequenceSize; - MITK_INFO << "TotalNumberOfSequences: " << TotalNumberOfSequences; + unsigned int totalNumberOfSequences = numberOfInputImages / sequenceSize; + MITK_INFO(m_Verbose) << "TotalNumberOfSequences: " << totalNumberOfSequences; - InitializeOutputs(TotalNumberOfSequences); + InitializeOutputs(totalNumberOfSequences); - auto EndmemberMatrix = CalculateEndmemberMatrix(m_Chromophore, m_Wavelength); //Eigen Endmember Matrix + auto endmemberMatrix = CalculateEndmemberMatrix(m_Chromophore, m_Wavelength); // Copy input image into array mitk::ImageReadAccessor readAccess(input); const float* inputDataArray = ((const float*)readAccess.GetData()); - CheckPreConditions(size, NumberOfInputImages, inputDataArray); + CheckPreConditions(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) + for (unsigned int sequenceCounter = 0; sequenceCounter < totalNumberOfSequences; ++sequenceCounter) { - - MITK_INFO << "SequenceCounter: " << SequenceCounter; - //loop over every pixel @ x,y plane + MITK_INFO(m_Verbose) << "SequenceCounter: " << sequenceCounter; + //loop over every pixel in XY-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; + /** + * 'sequenceCounter*sequenceSize' has to be added to 'z' to ensure that one accesses the + * correct pixel, because the inputDataArray contains the information of all sequences and + * not just the one of the current sequence. + */ + 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 + //write all wavelength absorbtion values for one(!) pixel of a sequence to a vector inputVector[z] = pixel; } - Eigen::VectorXf resultVector = SpectralUnmixingAlgorithm(EndmemberMatrix, inputVector); + 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]; + 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]"; + MITK_INFO(m_Verbose) << "GENERATING DATA...[DONE]"; myfile.close(); } -void mitk::pa::SpectralUnmixingFilterBase::CheckPreConditions(unsigned int size, unsigned int NumberOfInputImages, const float* inputDataArray) +void mitk::pa::SpectralUnmixingFilterBase::CheckPreConditions(unsigned int numberOfInputImages, const float* inputDataArray) { - if (m_Wavelength.size() < NumberOfInputImages) + MITK_INFO(m_Verbose) << "CHECK PRECONDITIONS ..."; + + if (m_Wavelength.size() < numberOfInputImages) MITK_WARN << "NUMBER OF WAVELENGTHS < NUMBER OF INPUT IMAGES"; - if (m_Wavelength.size() > NumberOfInputImages) + 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() << "ADD MORE WAVELENGTHS!"; - // Checking if pixel type is float if (typeid(inputDataArray[0]).name() != typeid(float).name()) mitkThrow() << "PIXELTYPE ERROR! FLOAT 32 REQUIRED"; - MITK_INFO << "CHECK PRECONDITIONS ...[DONE]"; + MITK_INFO(m_Verbose) << "...[DONE]"; } -void mitk::pa::SpectralUnmixingFilterBase::InitializeOutputs(unsigned int TotalNumberOfSequences) +void mitk::pa::SpectralUnmixingFilterBase::InitializeOutputs(unsigned int totalNumberOfSequences) { + MITK_INFO(m_Verbose) << "Initialize Outputs ..."; unsigned int numberOfInputs = GetNumberOfIndexedInputs(); unsigned int numberOfOutputs = GetNumberOfIndexedOutputs(); - MITK_INFO << "Inputs: " << numberOfInputs << " Outputs: " << numberOfOutputs; + MITK_INFO(m_Verbose) << "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; + 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); - } + MITK_INFO(m_Verbose) << "...[DONE]"; } - -//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); + 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!"; + for (unsigned int i = 0; i < numberOfWavelengths; ++i) + endmemberMatrixEigen(i, j) = PropertyElement(m_Chromophore[j], m_Wavelength[i]); } - MITK_INFO << "GENERATING ENMEMBERMATRIX [DONE]!"; - return EndmemberMatrixEigen; + MITK_INFO(m_Verbose) << "GENERATING ENMEMBERMATRIX [DONE]"; + return endmemberMatrixEigen; } -float mitk::pa::SpectralUnmixingFilterBase::PropertyElement(mitk::pa::PropertyCalculator::ChromophoreType Chromophore, int Wavelength) +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; + if (chromophore == mitk::pa::PropertyCalculator::ChromophoreType::ONEENDMEMBER) + return 1; + else + { + float value = m_PropertyCalculatorEigen->GetAbsorptionForWavelength(chromophore, wavelength); + if (value == 0) + mitkThrow() << "WAVELENGTH " << wavelength << "nm NOT SUPPORTED!"; + else + return value; + } }