diff --git a/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingFilterBase.cpp b/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingFilterBase.cpp index 240dc01dae..3d1177112b 100644 --- a/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingFilterBase.cpp +++ b/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingFilterBase.cpp @@ -1,240 +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; MITK_INFO(m_Verbose) << "TotalNumberOfSequences: " << totalNumberOfSequences; InitializeOutputs(totalNumberOfSequences); auto endmemberMatrix = CalculateEndmemberMatrix(m_Chromophore, m_Wavelength); // 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"; + MITK_WARN(m_Verbose) << "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()) + if ((m_Chromophore.size()+ m_RelativeError )!= GetNumberOfIndexedOutputs() || numberOfInputImages < 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 c766aad867..ca4b1f1fee 100644 --- a/Modules/PhotoacousticsLib/test/mitkSpectralUnmixingTest.cpp +++ b/Modules/PhotoacousticsLib/test/mitkSpectralUnmixingTest.cpp @@ -1,614 +1,648 @@ //*=================================================================== //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(testSimplexSUAlgorithm); MITK_TEST(testSO2); + MITK_TEST(testExceptionSO2); 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(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(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(false); + m_SpectralUnmixingFilter->Verbose(true); 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])SetInput(0, inputImage); + + inputImage = nullptr; + 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] = 4; + + inputImage->Initialize(pixelType, NUMBER_OF_SPATIAL_DIMENSIONS, dimensions); + + float* data = new float[3]; + data[0] = 1; + data[1] = 2; + data[2] = 3; + data[3] = 4; + + inputImage->SetImportVolume(data, mitk::Image::ImportMemoryManagementType::CopyMemory); + delete[] data; + + m_sO2->SetInput(1, inputImage); + + MITK_TEST_FOR_EXCEPTION_BEGIN(itk::ExceptionObject) + m_sO2->Update(); + MITK_TEST_FOR_EXCEPTION_END(itk::ExceptionObject) + } + // Tests SO2 Filter with correct inputs void testSO2() { MITK_INFO << "START SO2 TEST ... "; std::vector 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 << "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) 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 3a1ded9072..f2fe2cbd40 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,1094 +1,1097 @@ SpectralUnmixingControls 0 0 325 742 0 0 QmitkTemplate QLabel { color: rgb(255, 0, 0) } Please select an image! Qt::Horizontal 75 true Wavelengths settings 0 75 16777215 75 10 <html><head/><body><p>* \brief All values have to be intergers. They need to have the same order than the wavelength at the image.</p></body></html> true Qt::SolidLine 42 50 30 30 λ [nm] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 720 755 831 906 942 Qt::Horizontal 75 true <html><head/><body><p><span style=" font-weight:400;">* \brief Select as least one Absorber. It's not possible to select more absorbers then wavelengths.</span></p></body></html> Chromophore selection <html><head/><body><p>* \brief One of the main absorbers in near infrared spectrum.</p></body></html> Oxyhemoglobin true true <html><head/><body><p>* \brief One of the main absorbers in near infrared spectrum.</p></body></html> Deoxygenated hemoglobin true true false <html><head/><body><p>* \brief One of the main absorbers in near infrared spectrum.</p></body></html> Melanin <html><head/><body><p>* \brief This endmember will be unmixed with 1 at all wavelgnths.</p></body></html> Static Endmember Qt::Horizontal 75 true <html><head/><body><p><span style=" font-weight:400;">* \brief One needs to choose an spectral unmixing algorithm.</span></p></body></html> Unmixing algorithm true MS Shell Dlg 2 <html><head/><body><p>* \brief For detailed information about the algorithms please have a look at the documentation.</p></body></html> false 21 2147483647 ==CHOSE ALGORITHM== ==QR decomposition== householderQr colPivHouseholderQr fullPivHouseholderQr ==LU decompositon fullPivLu ==Cholesky decompostion== ldlt llt .. ==Least squares== LS jacobiSvd NNLARS NNGoldfarb weighted ==Others== SimplexMax true 0 75 380 81 <html><head/><body><p>* \brief the weights are at the same order as the wavelength</p></body></html> false false Weights [%] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 20 60 85 90 69 52 Qt::Horizontal 75 true Oxygen saturation <html><head/><body><p>* \brief calculates HbO2/(HbO2+Hb) if De- and oxyhemoglobin are selected.</p></body></html> calculate sO2 false 380 82 <html><head/><body><p>* \brief below threshold calculated sO2 value will set to zero</p></body></html> Threshold HbO2 Hb Sum SO2 % 200 100 200 50 Qt::Horizontal 75 true Additional Settings <html><head/><body><p>* \brief This mode will give additional console outputs for debugging.</p></body></html> Verbose Mode (additional console outputs) + + false + <html><head/><body><p>* \brief This checkbox will start Chrono and takes the time between clicking of the &quot;Perform spectral unmixing&quot; button until the GUI enables again.</p></body></html> Chrono <html><head/><body><p>* \brief Calculates the realtive error between unmixing result and the input image in the L2 norm.</p></body></html> Relative error image 307 70 <html><head/><body><p>* \brief below the threshold calculated relative error will set to zero</p></body></html> Threshold HbO2 Hb 150 200 Do image processing Perform spectral unmixing Qt::Horizontal Qt::Vertical 17 54