diff --git a/Modules/PhotoacousticsLib/include/mitkPASpectralUnmixingFilterBase.h b/Modules/PhotoacousticsLib/include/mitkPASpectralUnmixingFilterBase.h index 9a9d41d4eb..a06439862e 100644 --- a/Modules/PhotoacousticsLib/include/mitkPASpectralUnmixingFilterBase.h +++ b/Modules/PhotoacousticsLib/include/mitkPASpectralUnmixingFilterBase.h @@ -1,194 +1,200 @@ /*=================================================================== 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 { /** * \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. Furthermore it is possible to creat an output image that contains the information about the relative error between unmixing result * and the input image. * * 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. * - "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 + * - "ADD OUTPUTS HAS TO BE LARGER THEN ZERO!" + * - "NO WAVELENGHTS/CHROMOPHORES SELECZED! + * - "INDEX ERROR! NUMBER OF OUTPUTS DOESN'T FIT TO OTHER SETTIGNS!" */ class MITKPHOTOACOUSTICSLIB_EXPORT SpectralUnmixingFilterBase : public mitk::ImageToImageFilter { public: mitkClassMacro(SpectralUnmixingFilterBase, mitk::ImageToImageFilter); /** * \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. * @param wavelength database supports integers between 300 and 1000 nm */ void AddWavelength(int wavelength); /* * \brief Verbose gives more information to the console. Default value is false. * @param verbose is the boolian to activate the MITK_INFO logged to the console */ virtual void Verbose(bool verbose); /** * \brief AddOutputs takes an integer and sets indexed outputs * @param outputs integer correponds to the number of output images + * @throws if outputs == 0 */ virtual void AddOutputs(unsigned int outputs); /* * \brief RelativeError returns a image which compare the L2 norm of the input vector with the unmixing result * @param relativeError is the boolian to activate this tool */ virtual void RelativeError(bool relativeError); /** * \brief AddRelativeErrorSettings takes integers and writes them at the end of the m_RelativeErrorSettings vector. * @param value has to be a integer */ virtual void AddRelativeErrorSettings(int value); ofstream myfile; // just for testing purposes; has to be removeed protected: /** * \brief Constructor creats proptery calculater smart pointer new() */ SpectralUnmixingFilterBase(); virtual ~SpectralUnmixingFilterBase(); /** * \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. * @throws if algorithm implementiation fails (implemented for the algorithms with critical requirements) */ virtual Eigen::VectorXf SpectralUnmixingAlgorithm(Eigen::Matrix endmemberMatrix, Eigen::VectorXf inputVector) = 0; bool m_Verbose = false; bool m_RelativeError = false; std::vector m_Chromophore; std::vector m_Wavelength; std::vector m_RelativeErrorSettings; private: /* * \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. * @throws if there are more wavelength then images * @throws if there are more chromophores then wavelengths * @throws if the pixel type is not float 32 + * @throws if no chromophores or wavelengths selected as input + * @throws if the number of indexed outputs does'nt fit to the expected number */ 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; /* * \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); /* * \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 * @throws 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); /* * \brief calculates the relative error between the input image and the unmixing result in the L2 norm * @param endmemberMatrix is a Eigen matrix containing the endmember information * @param inputVector is a Eigen vector containing the multispectral information of one pixel * @param resultVector is a Eigen vector containing the spectral unmmixing result */ float CalculateRelativeError(Eigen::Matrix endmemberMatrix, Eigen::VectorXf inputVector, Eigen::VectorXf resultVector); PropertyCalculator::Pointer m_PropertyCalculatorEigen; }; } } #endif // MITKPHOTOACOUSTICSPECTRALUNMIXINGFILTERBASE_ diff --git a/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingFilterBase.cpp b/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingFilterBase.cpp index f676c53fbc..240dc01dae 100644 --- a/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingFilterBase.cpp +++ b/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingFilterBase.cpp @@ -1,234 +1,240 @@ /*=================================================================== 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 mitk::pa::SpectralUnmixingFilterBase::SpectralUnmixingFilterBase() { m_PropertyCalculatorEigen = mitk::pa::PropertyCalculator::New(); } mitk::pa::SpectralUnmixingFilterBase::~SpectralUnmixingFilterBase() { } void mitk::pa::SpectralUnmixingFilterBase::AddOutputs(unsigned int outputs) { + if (outputs == 0) + mitkThrow() << "ADD OUTPUTS HAS TO BE LARGER THEN ZERO!"; this->SetNumberOfIndexedOutputs(outputs); for (unsigned int i = 0; iSetNthOutput(i, mitk::Image::New()); } 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::Verbose(bool verbose) { m_Verbose = verbose; } void mitk::pa::SpectralUnmixingFilterBase::RelativeError(bool relativeError) { m_RelativeError = relativeError; } void mitk::pa::SpectralUnmixingFilterBase::AddRelativeErrorSettings(int value) { m_RelativeErrorSettings.push_back(value); } void mitk::pa::SpectralUnmixingFilterBase::GenerateData() { MITK_INFO(m_Verbose) << "GENERATING DATA.."; mitk::Image::Pointer input = GetInput(0); unsigned int xDim = input->GetDimensions()[0]; unsigned int yDim = input->GetDimensions()[1]; unsigned int numberOfInputImages = input->GetDimensions()[2]; MITK_INFO(m_Verbose) << "x dimension: " << xDim; MITK_INFO(m_Verbose) << "y dimension: " << yDim; MITK_INFO(m_Verbose) << "z dimension: " << numberOfInputImages; + // Copy input image into array + mitk::ImageReadAccessor readAccess(input); + const float* inputDataArray = ((const float*)readAccess.GetData()); + + CheckPreConditions(numberOfInputImages, inputDataArray); + unsigned int sequenceSize = m_Wavelength.size(); unsigned int totalNumberOfSequences = numberOfInputImages / sequenceSize; - if (totalNumberOfSequences == 0) //means that more chromophores then wavelengths - mitkThrow() << "ERROR! REMOVE WAVELENGTHS!"; MITK_INFO(m_Verbose) << "TotalNumberOfSequences: " << totalNumberOfSequences; InitializeOutputs(totalNumberOfSequences); auto endmemberMatrix = CalculateEndmemberMatrix(m_Chromophore, m_Wavelength); - - // Copy input image into array - mitk::ImageReadAccessor readAccess(input); - const float* inputDataArray = ((const float*)readAccess.GetData()); - - CheckPreConditions(numberOfInputImages, inputDataArray); // test to see pixel values @ txt file myfile.open("SimplexNormalisation.txt"); unsigned int outputCounter = GetNumberOfIndexedOutputs(); std::vector writteBufferVector; for (unsigned int i = 0; i < outputCounter; ++i) { auto output = GetOutput(i); mitk::ImageWriteAccessor writeOutput(output); float* writeBuffer = (float *)writeOutput.GetData(); writteBufferVector.push_back(writeBuffer); } if (m_RelativeError == true) { // -1 because rel error is output[IndexedOutputs() - 1] and loop over chromophore outputs has to end at [IndexedOutputs() - 2] outputCounter -= 1; } for (unsigned int sequenceCounter = 0; sequenceCounter < totalNumberOfSequences; ++sequenceCounter) { 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++) { /** * '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]; inputVector[z] = pixel; } Eigen::VectorXf resultVector = SpectralUnmixingAlgorithm(endmemberMatrix, inputVector); if (m_RelativeError == true) { float relativeError = CalculateRelativeError(endmemberMatrix, inputVector, resultVector); writteBufferVector[outputCounter][(xDim*yDim * sequenceCounter) + x * yDim + y] = relativeError; } for (unsigned int outputIdx = 0; outputIdx < outputCounter; ++outputIdx) { writteBufferVector[outputIdx][(xDim*yDim * sequenceCounter) + x * yDim + y] = resultVector[outputIdx]; } } } } MITK_INFO(m_Verbose) << "GENERATING DATA...[DONE]"; myfile.close(); } void mitk::pa::SpectralUnmixingFilterBase::CheckPreConditions(unsigned int numberOfInputImages, const float* inputDataArray) { MITK_INFO(m_Verbose) << "CHECK PRECONDITIONS ..."; + if (m_Chromophore.size() == 0 || m_Wavelength.size() == 0) + mitkThrow() << "NO WAVELENGHTS/CHROMOPHORES SELECTED!"; + if (m_Wavelength.size() < numberOfInputImages) MITK_WARN << "NUMBER OF WAVELENGTHS < NUMBER OF INPUT IMAGES"; if (m_Chromophore.size() > m_Wavelength.size()) mitkThrow() << "ADD MORE WAVELENGTHS OR REMOVE ENDMEMBERS!"; if (typeid(inputDataArray[0]).name() != typeid(float).name()) mitkThrow() << "PIXELTYPE ERROR! FLOAT 32 REQUIRED"; + if ((m_Chromophore.size()+ m_RelativeError )!= GetNumberOfIndexedOutputs()) + mitkThrow() << "INDEX ERROR! NUMBER OF OUTPUTS DOESN'T FIT TO OTHER SETTIGNS!"; + MITK_INFO(m_Verbose) << "...[DONE]"; } void mitk::pa::SpectralUnmixingFilterBase::InitializeOutputs(unsigned int totalNumberOfSequences) { MITK_INFO(m_Verbose) << "Initialize Outputs ..."; unsigned int numberOfInputs = GetNumberOfIndexedInputs(); unsigned int numberOfOutputs = GetNumberOfIndexedOutputs(); MITK_INFO(m_Verbose) << "Inputs: " << numberOfInputs << " Outputs: " << numberOfOutputs; 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; for (unsigned int outputIdx = 0; outputIdx < numberOfOutputs; outputIdx++) GetOutput(outputIdx)->Initialize(pixelType, NUMBER_OF_SPATIAL_DIMENSIONS, dimensions); MITK_INFO(m_Verbose) << "...[DONE]"; } 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) { for (unsigned int i = 0; i < numberOfWavelengths; ++i) endmemberMatrixEigen(i, j) = PropertyElement(m_Chromophore[j], m_Wavelength[i]); } MITK_INFO(m_Verbose) << "GENERATING ENMEMBERMATRIX [DONE]"; return endmemberMatrixEigen; } float mitk::pa::SpectralUnmixingFilterBase::PropertyElement(mitk::pa::PropertyCalculator::ChromophoreType chromophore, int wavelength) { 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; } } float mitk::pa::SpectralUnmixingFilterBase::CalculateRelativeError(Eigen::Matrix endmemberMatrix, Eigen::VectorXf inputVector, Eigen::VectorXf resultVector) { float relativeError = (endmemberMatrix*resultVector - inputVector).norm() / inputVector.norm(); for (int i = 0; i < 2; ++i) { if (resultVector[i] < m_RelativeErrorSettings[i]) return 0; } return relativeError; } diff --git a/Modules/PhotoacousticsLib/test/mitkSpectralUnmixingTest.cpp b/Modules/PhotoacousticsLib/test/mitkSpectralUnmixingTest.cpp index 2980b69d05..c766aad867 100644 --- a/Modules/PhotoacousticsLib/test/mitkSpectralUnmixingTest.cpp +++ b/Modules/PhotoacousticsLib/test/mitkSpectralUnmixingTest.cpp @@ -1,336 +1,614 @@ //*=================================================================== //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 #include #include #include #include #include #include #include class mitkSpectralUnmixingTestSuite : public mitk::TestFixture { CPPUNIT_TEST_SUITE(mitkSpectralUnmixingTestSuite); MITK_TEST(testEigenSUAlgorithm); MITK_TEST(testVigraSUAlgorithm); //MITK_TEST(testSimplexSUAlgorithm);// --> RESULT FAILS MITK_TEST(testSO2); + MITK_TEST(testWavelengthExceptions); + MITK_TEST(testNoChromophoresSelected); + MITK_TEST(testInputImage); + MITK_TEST(testAddOutput); + MITK_TEST(testWeightsError); CPPUNIT_TEST_SUITE_END(); private: mitk::pa::SpectralUnmixingFilterBase::Pointer m_SpectralUnmixingFilter; mitk::Image::Pointer inputImage; std::vector m_inputWavelengths; std::vector m_inputWeights; std::vector m_CorrectResult; float threshold; public: void setUp() override { MITK_INFO << "setUp ... "; //Set empty input image: inputImage = mitk::Image::New(); mitk::PixelType pixelType = mitk::MakeScalarPixelType(); const int NUMBER_OF_SPATIAL_DIMENSIONS = 3; auto* dimensions = new unsigned int[NUMBER_OF_SPATIAL_DIMENSIONS]; dimensions[0] = 1; dimensions[1] = 1; dimensions[2] = 5; inputImage->Initialize(pixelType, NUMBER_OF_SPATIAL_DIMENSIONS, dimensions); //Set wavelengths for unmixing: m_inputWavelengths.push_back(750); m_inputWavelengths.push_back(800); m_inputWeights.push_back(50); m_inputWeights.push_back(100); //Set fraction of Hb and HbO2 to unmix: float fracHb = 100; float fracHbO2 = 300; m_CorrectResult.push_back(fracHbO2); m_CorrectResult.push_back(fracHb); m_CorrectResult.push_back(fracHbO2 + 10); m_CorrectResult.push_back(fracHb - 10); threshold = 0.01; //Multiply values of wavelengths (750,800,850 nm) with fractions to get pixel values: float px1 = fracHb * 7.52 + fracHbO2 * 2.77; float px2 = fracHb * 4.08 + fracHbO2 * 4.37; float px3 = (fracHb - 10) * 7.52 + (fracHbO2 + 10) * 2.77; float px4 = (fracHb - 10) * 4.08 + (fracHbO2 + 10) * 4.37; float* data = new float[3]; data[0] = px1; data[1] = px2; data[2] = px3; data[3] = px4; data[5] = 0; inputImage->SetImportVolume(data, mitk::Image::ImportMemoryManagementType::CopyMemory); delete[] data; MITK_INFO << "[DONE]"; } + // Tests implemented EIGEN algortihms with correct inputs void testEigenSUAlgorithm() { MITK_INFO << "START FILTER TEST ... "; // Set input image auto m_SpectralUnmixingFilter = mitk::pa::LinearSpectralUnmixingFilter::New(); - m_SpectralUnmixingFilter->Verbose(true); + m_SpectralUnmixingFilter->Verbose(false); m_SpectralUnmixingFilter->RelativeError(false); m_SpectralUnmixingFilter->SetInput(inputImage); m_SpectralUnmixingFilter->AddOutputs(2); //Set wavelengths to filter for (unsigned int imageIndex = 0; imageIndex < m_inputWavelengths.size(); imageIndex++) { unsigned int wavelength = m_inputWavelengths[imageIndex]; m_SpectralUnmixingFilter->AddWavelength(wavelength); } //Set Chromophores to filter m_SpectralUnmixingFilter->AddChromophore( mitk::pa::PropertyCalculator::ChromophoreType::OXYGENATED); m_SpectralUnmixingFilter->AddChromophore( mitk::pa::PropertyCalculator::ChromophoreType::DEOXYGENATED); ofstream myfile; myfile.open("EigenTestResult.txt"); std::vector m_Eigen = { mitk::pa::LinearSpectralUnmixingFilter::AlgortihmType::HOUSEHOLDERQR, /* mitk::pa::LinearSpectralUnmixingFilter::AlgortihmType::LDLT, mitk::pa::LinearSpectralUnmixingFilter::AlgortihmType::LLT,*/ mitk::pa::LinearSpectralUnmixingFilter::AlgortihmType::COLPIVHOUSEHOLDERQR, mitk::pa::LinearSpectralUnmixingFilter::AlgortihmType::JACOBISVD, mitk::pa::LinearSpectralUnmixingFilter::AlgortihmType::FULLPIVLU, mitk::pa::LinearSpectralUnmixingFilter::AlgortihmType::FULLPIVHOUSEHOLDERQR}; for (int Algorithmidx = 0; Algorithmidx < m_Eigen.size();++Algorithmidx) { m_SpectralUnmixingFilter->SetAlgorithm(m_Eigen[Algorithmidx]); m_SpectralUnmixingFilter->Update(); /*For printed pixel values and results look at: [...]\mitk-superbuild\MITK-build\Modules\PhotoacousticsLib\test\*/ for (int i = 0; i < 2; ++i) { mitk::Image::Pointer output = m_SpectralUnmixingFilter->GetOutput(i); mitk::ImageReadAccessor readAccess(output); const float* inputDataArray = ((const float*)readAccess.GetData()); auto pixel = inputDataArray[0]; auto pixel2 = inputDataArray[1]; myfile << "Algorithmidx: " << Algorithmidx << "\n Output " << i << ": " << "\n" << inputDataArray[0] << "\n" << inputDataArray[1] << "\n"; myfile << "Correct Result: " << "\n" << m_CorrectResult[i] << "\n" << m_CorrectResult[i + 2] << "\n"; CPPUNIT_ASSERT(std::abs(pixel - m_CorrectResult[i]) < threshold); CPPUNIT_ASSERT(std::abs(pixel2 - m_CorrectResult[i + 2]) < threshold); } } myfile.close(); MITK_INFO << "EIGEN FILTER TEST SUCCESFULL :)"; } + // Tests implemented VIGRA algortihms with correct inputs void testVigraSUAlgorithm() { MITK_INFO << "START FILTER TEST ... "; // Set input image auto m_SpectralUnmixingFilter = mitk::pa::SpectralUnmixingFilterVigra::New(); m_SpectralUnmixingFilter->SetInput(inputImage); m_SpectralUnmixingFilter->AddOutputs(2); - m_SpectralUnmixingFilter->Verbose(true); + m_SpectralUnmixingFilter->Verbose(false); m_SpectralUnmixingFilter->RelativeError(false); //Set wavelengths to filter for (unsigned int imageIndex = 0; imageIndex < m_inputWavelengths.size(); imageIndex++) { unsigned int wavelength = m_inputWavelengths[imageIndex]; double Weight = m_inputWeights[imageIndex]; m_SpectralUnmixingFilter->AddWavelength(wavelength); m_SpectralUnmixingFilter->AddWeight(Weight); } //Set Chromophores to filter m_SpectralUnmixingFilter->AddChromophore( mitk::pa::PropertyCalculator::ChromophoreType::OXYGENATED); m_SpectralUnmixingFilter->AddChromophore( mitk::pa::PropertyCalculator::ChromophoreType::DEOXYGENATED); std::vector Vigra = { mitk::pa::SpectralUnmixingFilterVigra::VigraAlgortihmType::LARS, mitk::pa::SpectralUnmixingFilterVigra::VigraAlgortihmType::GOLDFARB, mitk::pa::SpectralUnmixingFilterVigra::VigraAlgortihmType::WEIGHTED, mitk::pa::SpectralUnmixingFilterVigra::VigraAlgortihmType::LS/*, mitk::pa::SpectralUnmixingFilterVigra::VigraAlgortihmType::vigratest*/}; ofstream myfile; myfile.open("VigraTestResult.txt"); for (int Algorithmidx = 0; Algorithmidx < Vigra.size();++Algorithmidx) { m_SpectralUnmixingFilter->SetAlgorithm(Vigra[0]); m_SpectralUnmixingFilter->Update(); /*For printed pixel values and results look at: [...]\mitk-superbuild\MITK-build\Modules\PhotoacousticsLib\test\*/ for (int i = 0; i < 2; ++i) { mitk::Image::Pointer output = m_SpectralUnmixingFilter->GetOutput(i); mitk::ImageReadAccessor readAccess(output); const float* inputDataArray = ((const float*)readAccess.GetData()); auto pixel = inputDataArray[0]; auto pixel2 = inputDataArray[1]; myfile << "Algorithmidx: " << Algorithmidx << "\n Output " << i << ": " << "\n" << inputDataArray[0] << "\n" << inputDataArray[1] << "\n"; myfile << "Correct Result: " << "\n" << m_CorrectResult[i] << "\n" << m_CorrectResult[i + 2] << "\n"; CPPUNIT_ASSERT(std::abs(pixel - m_CorrectResult[i]) < threshold); CPPUNIT_ASSERT(std::abs(pixel2 - m_CorrectResult[i + 2]) < threshold); } } myfile.close(); MITK_INFO << "VIGRA FILTER TEST SUCCESFULL :)"; } void testSimplexSUAlgorithm() { MITK_INFO << "START FILTER TEST ... "; // Set input image auto m_SpectralUnmixingFilter = mitk::pa::SpectralUnmixingFilterSimplex::New(); m_SpectralUnmixingFilter->SetInput(inputImage); m_SpectralUnmixingFilter->AddOutputs(2); - m_SpectralUnmixingFilter->Verbose(true); + m_SpectralUnmixingFilter->Verbose(false); m_SpectralUnmixingFilter->RelativeError(false); //Set wavelengths to filter for (unsigned int imageIndex = 0; imageIndex < m_inputWavelengths.size(); imageIndex++) { unsigned int wavelength = m_inputWavelengths[imageIndex]; m_SpectralUnmixingFilter->AddWavelength(wavelength); } //Set Chromophores to filter m_SpectralUnmixingFilter->AddChromophore( mitk::pa::PropertyCalculator::ChromophoreType::OXYGENATED); m_SpectralUnmixingFilter->AddChromophore( mitk::pa::PropertyCalculator::ChromophoreType::DEOXYGENATED); m_SpectralUnmixingFilter->Update(); /*For printed pixel values and results look at: [...]\mitk-superbuild\MITK-build\Modules\PhotoacousticsLib\test\*/ ofstream myfile; myfile.open("SimplexTestResult.txt"); for (int i = 0; i < 2; ++i) { mitk::Image::Pointer output = m_SpectralUnmixingFilter->GetOutput(i); mitk::ImageReadAccessor readAccess(output); const float* inputDataArray = ((const float*)readAccess.GetData()); auto pixel = inputDataArray[0]; auto pixel2 = inputDataArray[1]; myfile << "Output " << i << ": " << "\n" << inputDataArray[0] << "\n" << inputDataArray[1] << "\n"; myfile << "Correct Result: " << "\n" << m_CorrectResult[i] << "\n" << m_CorrectResult[i + 2] << "\n"; CPPUNIT_ASSERT(std::abs(pixel - m_CorrectResult[i]) CorrectSO2Result1 = { 0, 0, 0, 0, 0 }; std::vector Test1 = { 0,0,0,51 }; std::vector CorrectSO2Result2 = { 0, 0.5, 0, 0.5, 0 }; std::vector Test2 = { 1584, 0, 0, 0 }; std::vector CorrectSO2Result3 = { 0.5, 0.5, 0, 0.5, 0 }; std::vector Test3 = { 0, 1536, 0, 0 }; std::vector CorrectSO2Result4 = { 0.5, 0.5, 0, 0.5, 0 }; std::vector Test4 = { 0, 0, 3072, 49 }; std::vector CorrectSO2Result5 = { 0.5, 0.5, 0.5, 0.5, 0 }; std::vector Test5 = { 1, 1, 1, 49 }; std::vector> TestList; std::vector> ResultList; TestList.push_back(Test1); TestList.push_back(Test2); TestList.push_back(Test3); TestList.push_back(Test4); TestList.push_back(Test5); ResultList.push_back(CorrectSO2Result1); ResultList.push_back(CorrectSO2Result2); ResultList.push_back(CorrectSO2Result3); ResultList.push_back(CorrectSO2Result4); ResultList.push_back(CorrectSO2Result5); /*For printed pixel values and results look at: [...]\mitk-superbuild\MITK-build\Modules\PhotoacousticsLib\test\*/ ofstream myfile; myfile.open("SO2TestResult.txt"); for (int k = 0; k < 5; ++k) { std::vector SO2Settings = TestList[k]; std::vector m_CorrectSO2Result = ResultList[k]; auto m_sO2 = mitk::pa::SpectralUnmixingSO2::New(); m_sO2->SetInput(0, inputImage); m_sO2->SetInput(1, inputImage); for (int i = 0; i < SO2Settings.size(); ++i) m_sO2->AddSO2Settings(SO2Settings[i]); m_sO2->Update(); mitk::Image::Pointer output = m_sO2->GetOutput(0); mitk::ImageReadAccessor readAccess(output); const float* inputDataArray = ((const float*)readAccess.GetData()); for (unsigned int Pixel = 0; Pixel < inputImage->GetDimensions()[2]; ++Pixel) { auto Value = inputDataArray[Pixel]; myfile << "Output(Test "<< k <<") "<< Pixel << ": " << "\n" << Value << "\n"; myfile << "Correct Result: " << "\n" << m_CorrectSO2Result[Pixel] << "\n"; CPPUNIT_ASSERT(std::abs(Value - m_CorrectSO2Result[Pixel]) < threshold); } } myfile.close(); MITK_INFO << "SO2 TEST SUCCESFULL :)"; } + // Test exceptions for wrong wavelength inputs + void testWavelengthExceptions() + { + MITK_INFO << "START WavelengthExceptions TEST ... "; + // Set input image + auto m_SpectralUnmixingFilter = mitk::pa::LinearSpectralUnmixingFilter::New(); + m_SpectralUnmixingFilter->Verbose(false); + m_SpectralUnmixingFilter->RelativeError(false); + m_SpectralUnmixingFilter->SetInput(inputImage); + m_SpectralUnmixingFilter->AddOutputs(2); + + m_SpectralUnmixingFilter->AddChromophore( + mitk::pa::PropertyCalculator::ChromophoreType::OXYGENATED); + m_SpectralUnmixingFilter->AddChromophore( + mitk::pa::PropertyCalculator::ChromophoreType::DEOXYGENATED); + + m_SpectralUnmixingFilter->SetAlgorithm(mitk::pa::LinearSpectralUnmixingFilter::AlgortihmType::HOUSEHOLDERQR); + + MITK_TEST_FOR_EXCEPTION_BEGIN(itk::ExceptionObject) + m_SpectralUnmixingFilter->Update(); + MITK_TEST_FOR_EXCEPTION_END(itk::ExceptionObject) + + m_SpectralUnmixingFilter->AddWavelength(300); + + MITK_TEST_FOR_EXCEPTION_BEGIN(itk::ExceptionObject) + m_SpectralUnmixingFilter->Update(); + MITK_TEST_FOR_EXCEPTION_END(itk::ExceptionObject) + + m_SpectralUnmixingFilter->AddWavelength(299); + + MITK_TEST_FOR_EXCEPTION_BEGIN(itk::ExceptionObject) + m_SpectralUnmixingFilter->Update(); + MITK_TEST_FOR_EXCEPTION_END(itk::ExceptionObject) + MITK_INFO << "DONE"; + } + + // Test exceptions for wrong chromophore inputs + void testNoChromophoresSelected() + { + MITK_INFO << "testNoChromophoresSelected"; + // Set input image + auto m_SpectralUnmixingFilter = mitk::pa::LinearSpectralUnmixingFilter::New(); + m_SpectralUnmixingFilter->Verbose(false); + m_SpectralUnmixingFilter->RelativeError(false); + m_SpectralUnmixingFilter->SetInput(inputImage); + m_SpectralUnmixingFilter->AddOutputs(2); + + //Set wavelengths to filter + for (unsigned int imageIndex = 0; imageIndex < m_inputWavelengths.size(); imageIndex++) + { + unsigned int wavelength = m_inputWavelengths[imageIndex]; + m_SpectralUnmixingFilter->AddWavelength(wavelength); + } + + m_SpectralUnmixingFilter->SetAlgorithm(mitk::pa::LinearSpectralUnmixingFilter::AlgortihmType::HOUSEHOLDERQR); + + MITK_TEST_FOR_EXCEPTION_BEGIN(itk::ExceptionObject) + m_SpectralUnmixingFilter->Update(); + MITK_TEST_FOR_EXCEPTION_END(itk::ExceptionObject) + } + + // Test exceptions for wrong input image + void testInputImage() + { + MITK_INFO << "INPUT IMAGE TEST"; + inputImage = nullptr; + // Set input image + auto m_SpectralUnmixingFilter = mitk::pa::LinearSpectralUnmixingFilter::New(); + m_SpectralUnmixingFilter->Verbose(false); + m_SpectralUnmixingFilter->RelativeError(false); + //m_SpectralUnmixingFilter->SetInput(inputImage); + m_SpectralUnmixingFilter->AddOutputs(2); + + //Set wavelengths to filter + for (unsigned int imageIndex = 0; imageIndex < m_inputWavelengths.size(); imageIndex++) + { + unsigned int wavelength = m_inputWavelengths[imageIndex]; + m_SpectralUnmixingFilter->AddWavelength(wavelength); + } + + m_SpectralUnmixingFilter->AddChromophore( + mitk::pa::PropertyCalculator::ChromophoreType::OXYGENATED); + m_SpectralUnmixingFilter->AddChromophore( + mitk::pa::PropertyCalculator::ChromophoreType::DEOXYGENATED); + + m_SpectralUnmixingFilter->SetAlgorithm(mitk::pa::LinearSpectralUnmixingFilter::AlgortihmType::HOUSEHOLDERQR); + + MITK_TEST_FOR_EXCEPTION_BEGIN(itk::ExceptionObject) + m_SpectralUnmixingFilter->Update(); + MITK_TEST_FOR_EXCEPTION_END(itk::ExceptionObject) + + inputImage = mitk::Image::New(); + mitk::PixelType pixelType = mitk::MakeScalarPixelType(); + const int NUMBER_OF_SPATIAL_DIMENSIONS = 3; + auto* dimensions = new unsigned int[NUMBER_OF_SPATIAL_DIMENSIONS]; + + dimensions[0] = 1; + dimensions[1] = 1; + dimensions[2] = 5; + + inputImage->Initialize(pixelType, NUMBER_OF_SPATIAL_DIMENSIONS, dimensions); + + double* data = new double[3]; + data[0] = 1; + data[1] = 2; + data[2] = 3; + data[3] = 4; + data[5] = 0; + + inputImage->SetImportVolume(data, mitk::Image::ImportMemoryManagementType::CopyMemory); + delete[] data; + + m_SpectralUnmixingFilter->SetInput(inputImage); + + MITK_TEST_FOR_EXCEPTION_BEGIN(itk::ExceptionObject) + m_SpectralUnmixingFilter->Update(); + MITK_TEST_FOR_EXCEPTION_END(itk::ExceptionObject) + + for (int i = 0; i < 2; ++i) + { + mitk::Image::Pointer output = m_SpectralUnmixingFilter->GetOutput(i); + mitk::ImageReadAccessor readAccess(output); + const float* inputDataArray = ((const float*)readAccess.GetData()); + auto pixel = inputDataArray[0]; + auto pixel2 = inputDataArray[1]; + + CPPUNIT_ASSERT(std::abs(pixel - m_CorrectResult[i]) < threshold); + CPPUNIT_ASSERT(std::abs(pixel2 - m_CorrectResult[i + 2]) < threshold); + } + } + + // Test exceptions for addOutputs method + void testAddOutput() + { + MITK_INFO << "addOutputs TEST"; + + // Set input image + auto m_SpectralUnmixingFilter = mitk::pa::LinearSpectralUnmixingFilter::New(); + m_SpectralUnmixingFilter->Verbose(false); + m_SpectralUnmixingFilter->RelativeError(false); + m_SpectralUnmixingFilter->SetInput(inputImage); + + //Set wavelengths to filter + for (unsigned int imageIndex = 0; imageIndex < m_inputWavelengths.size(); imageIndex++) + { + unsigned int wavelength = m_inputWavelengths[imageIndex]; + m_SpectralUnmixingFilter->AddWavelength(wavelength); + } + + m_SpectralUnmixingFilter->AddChromophore( + mitk::pa::PropertyCalculator::ChromophoreType::OXYGENATED); + m_SpectralUnmixingFilter->AddChromophore( + mitk::pa::PropertyCalculator::ChromophoreType::DEOXYGENATED); + + m_SpectralUnmixingFilter->SetAlgorithm(mitk::pa::LinearSpectralUnmixingFilter::AlgortihmType::HOUSEHOLDERQR); + + for (int i = 0; i < 4; ++i) + { + MITK_INFO << "i: " << i; + if (i != 2) + { + MITK_TEST_FOR_EXCEPTION_BEGIN(itk::ExceptionObject) + m_SpectralUnmixingFilter->AddOutputs(i); + m_SpectralUnmixingFilter->Update(); + MITK_TEST_FOR_EXCEPTION_END(itk::ExceptionObject) + } + else + { + m_SpectralUnmixingFilter->AddOutputs(2); + m_SpectralUnmixingFilter->Update(); + for (int i = 0; i < 2; ++i) + { + mitk::Image::Pointer output = m_SpectralUnmixingFilter->GetOutput(i); + mitk::ImageReadAccessor readAccess(output); + const float* inputDataArray = ((const float*)readAccess.GetData()); + auto pixel = inputDataArray[0]; + auto pixel2 = inputDataArray[1]; + + CPPUNIT_ASSERT(std::abs(pixel - m_CorrectResult[i]) < threshold); + CPPUNIT_ASSERT(std::abs(pixel2 - m_CorrectResult[i + 2]) < threshold); + } + } + } + } + + // Test exceptions for weights error + void testWeightsError() + { + MITK_INFO << "testWeightsError"; + + // Set input image + auto m_SpectralUnmixingFilter = mitk::pa::SpectralUnmixingFilterVigra::New(); + m_SpectralUnmixingFilter->Verbose(false); + m_SpectralUnmixingFilter->RelativeError(false); + m_SpectralUnmixingFilter->SetInput(inputImage); + m_SpectralUnmixingFilter->AddOutputs(2); + + //Set wavelengths to filter + for (unsigned int imageIndex = 0; imageIndex < m_inputWavelengths.size(); imageIndex++) + { + unsigned int wavelength = m_inputWavelengths[imageIndex]; + m_SpectralUnmixingFilter->AddWavelength(wavelength); + } + + m_SpectralUnmixingFilter->AddChromophore( + mitk::pa::PropertyCalculator::ChromophoreType::OXYGENATED); + m_SpectralUnmixingFilter->AddChromophore( + mitk::pa::PropertyCalculator::ChromophoreType::DEOXYGENATED); + + m_SpectralUnmixingFilter->SetAlgorithm(mitk::pa::SpectralUnmixingFilterVigra::VigraAlgortihmType::WEIGHTED); + + MITK_TEST_FOR_EXCEPTION_BEGIN(itk::ExceptionObject) + m_SpectralUnmixingFilter->Update(); + MITK_TEST_FOR_EXCEPTION_END(itk::ExceptionObject) + + m_SpectralUnmixingFilter->AddWeight(50); + MITK_TEST_FOR_EXCEPTION_BEGIN(itk::ExceptionObject) + m_SpectralUnmixingFilter->Update(); + MITK_TEST_FOR_EXCEPTION_END(itk::ExceptionObject) + + m_SpectralUnmixingFilter->AddWeight(50); + m_SpectralUnmixingFilter->Update(); + + for (int i = 0; i < 2; ++i) + { + mitk::Image::Pointer output = m_SpectralUnmixingFilter->GetOutput(i); + mitk::ImageReadAccessor readAccess(output); + const float* inputDataArray = ((const float*)readAccess.GetData()); + auto pixel = inputDataArray[0]; + auto pixel2 = inputDataArray[1]; + + CPPUNIT_ASSERT(std::abs(pixel - m_CorrectResult[i]) < threshold); + CPPUNIT_ASSERT(std::abs(pixel2 - m_CorrectResult[i + 2]) < threshold); + } + } + + // TEST TEMPLATE: + /* + // Test exceptions for + void test() + { + MITK_INFO << "TEST"; + + // Set input image + auto m_SpectralUnmixingFilter = mitk::pa::LinearSpectralUnmixingFilter::New(); + m_SpectralUnmixingFilter->Verbose(false); + m_SpectralUnmixingFilter->RelativeError(false); + m_SpectralUnmixingFilter->SetInput(inputImage); + m_SpectralUnmixingFilter->AddOutputs(2); + + //Set wavelengths to filter + for (unsigned int imageIndex = 0; imageIndex < m_inputWavelengths.size(); imageIndex++) + { + unsigned int wavelength = m_inputWavelengths[imageIndex]; + m_SpectralUnmixingFilter->AddWavelength(wavelength); + } + + m_SpectralUnmixingFilter->AddChromophore( + mitk::pa::PropertyCalculator::ChromophoreType::OXYGENATED); + m_SpectralUnmixingFilter->AddChromophore( + mitk::pa::PropertyCalculator::ChromophoreType::DEOXYGENATED); + + m_SpectralUnmixingFilter->SetAlgorithm(mitk::pa::LinearSpectralUnmixingFilter::AlgortihmType::HOUSEHOLDERQR); + + //MITK_TEST_FOR_EXCEPTION_BEGIN(itk::ExceptionObject) + m_SpectralUnmixingFilter->Update(); + //MITK_TEST_FOR_EXCEPTION_END(itk::ExceptionObject) + + }*/ + void tearDown() override { m_SpectralUnmixingFilter = nullptr; inputImage = nullptr; m_inputWavelengths.clear(); m_CorrectResult.clear(); MITK_INFO << "tearDown ... [DONE]"; } }; MITK_TEST_SUITE_REGISTRATION(mitkSpectralUnmixing)