diff --git a/Modules/PhotoacousticsLib/include/mitkPASpectralUnmixingSO2.h b/Modules/PhotoacousticsLib/include/mitkPASpectralUnmixingSO2.h index e6c54f709e..12deffc995 100644 --- a/Modules/PhotoacousticsLib/include/mitkPASpectralUnmixingSO2.h +++ b/Modules/PhotoacousticsLib/include/mitkPASpectralUnmixingSO2.h @@ -1,51 +1,50 @@ /*=================================================================== 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 MITKPHOTOACOUSTICSPECTRALUNMIXINGSO2_H #define MITKPHOTOACOUSTICSPECTRALUNMIXINGSO2_H #include "mitkImageToImageFilter.h" #include //Includes for smart pointer usage #include "mitkCommon.h" #include "itkLightObject.h" namespace mitk { namespace pa { class MITKPHOTOACOUSTICSLIB_EXPORT SpectralUnmixingSO2 : public mitk::ImageToImageFilter { public: mitkClassMacro(SpectralUnmixingSO2, mitk::ImageToImageFilter); itkFactorylessNewMacro(Self); virtual void AddSO2Settings(float value); protected: SpectralUnmixingSO2(); virtual ~SpectralUnmixingSO2(); std::vector m_SO2Settings; - private: virtual void GenerateData() override; virtual void InitializeOutputs(); virtual void CheckPreConditions(mitk::Image::Pointer inputHbO2, mitk::Image::Pointer inputHb); float calculateSO2(float pixelHb, float pixelHbO2); }; } } #endif // MITKPHOTOACOUSTICSPECTRALUNMIXINGSO2_ diff --git a/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingSO2.cpp b/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingSO2.cpp index d9fd0d043f..514b90371a 100644 --- a/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingSO2.cpp +++ b/Modules/PhotoacousticsLib/src/SUFilter/mitkPASpectralUnmixingSO2.cpp @@ -1,154 +1,155 @@ /*=================================================================== 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 "mitkPASpectralUnmixingSO2.h" // ImageAccessor #include #include mitk::pa::SpectralUnmixingSO2::SpectralUnmixingSO2() { this->SetNumberOfIndexedInputs(2); this->SetNumberOfIndexedOutputs(1); this->SetNthOutput(0, mitk::Image::New()); } mitk::pa::SpectralUnmixingSO2::~SpectralUnmixingSO2() { } void mitk::pa::SpectralUnmixingSO2::GenerateData() { MITK_INFO << "GENERATING DATA.."; // Get input image mitk::Image::Pointer inputHbO2 = GetInput(0); mitk::Image::Pointer inputHb = GetInput(1); CheckPreConditions(inputHbO2, inputHb); unsigned int xDim = inputHbO2->GetDimensions()[0]; unsigned int yDim = inputHbO2->GetDimensions()[1]; unsigned int zDim = inputHbO2->GetDimensions()[2]; unsigned int size = xDim * yDim*zDim; MITK_INFO << "x dimension: " << xDim; MITK_INFO << "y dimension: " << yDim; MITK_INFO << "z dimension: " << zDim; InitializeOutputs(); // Copy input image into array mitk::ImageReadAccessor readAccessHbO2(inputHbO2); mitk::ImageReadAccessor readAccessHb(inputHb); const float* inputDataArrayHbO2 = ((const float*)readAccessHbO2.GetData()); const float* inputDataArrayHb = ((const float*)readAccessHb.GetData()); for (unsigned int x = 0; x < xDim; x++) { for (unsigned int y = 0; y < yDim; y++) { for (unsigned int z = 0;z < zDim; z++) { // Calculate SO2 and write to output unsigned int pixelNumber = (xDim*yDim * z) + x * yDim + y; float pixelHb = inputDataArrayHb[pixelNumber]; float pixelHbO2 = inputDataArrayHbO2[pixelNumber]; float resultSO2 = calculateSO2(pixelHb, pixelHbO2); auto output = GetOutput(0); mitk::ImageWriteAccessor writeOutput(output); float* writeBuffer = (float *)writeOutput.GetData(); writeBuffer[(xDim*yDim * z) + x * yDim + y] = resultSO2; } } } MITK_INFO << "GENERATING DATA...[DONE]"; } void mitk::pa::SpectralUnmixingSO2::CheckPreConditions(mitk::Image::Pointer inputHbO2, mitk::Image::Pointer inputHb) { unsigned int xDimHb = inputHb->GetDimensions()[0]; unsigned int yDimHb = inputHb->GetDimensions()[1]; unsigned int zDimHb = inputHb->GetDimensions()[2]; unsigned int xDimHbO2 = inputHbO2->GetDimensions()[0]; unsigned int yDimHbO2 = inputHbO2->GetDimensions()[1]; unsigned int zDimHbO2 = inputHbO2->GetDimensions()[2]; if (xDimHb != xDimHbO2 || yDimHb != yDimHbO2 || zDimHb != zDimHbO2) mitkThrow() << "DIMENTIONALITY ERROR!"; MITK_INFO << "CHECK PRECONDITIONS ...[DONE]"; } void mitk::pa::SpectralUnmixingSO2::InitializeOutputs() { unsigned int numberOfInputs = GetNumberOfIndexedInputs(); unsigned int numberOfOutputs = GetNumberOfIndexedOutputs(); // 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; dimIdxGetDimensions()[dimIdx]; } // 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); } } float mitk::pa::SpectralUnmixingSO2::calculateSO2(float Hb, float HbO2) { float MinHb = m_SO2Settings[0]; float MinHbO2 = m_SO2Settings[1]; float MinSum = m_SO2Settings[2]; float MinSO2 = m_SO2Settings[3]; float result = HbO2 / (Hb + HbO2); float zero = 0.0; if (result != result) { MITK_WARN << "SO2 VALUE NAN! WILL BE SET TO ZERO!"; return zero; } - else { + else + { if (MinHb != 0 && MinHb > Hb) return zero; else if (MinHbO2 != 0 && MinHbO2 > HbO2) return zero; else if (MinSum != 0 && MinSum > Hb + HbO2) return zero; else if (MinSO2 != 0 && 100 * result > MinSO2) return result; else if (MinSO2 != 0 && result < MinSO2) return zero; else return result; } } void mitk::pa::SpectralUnmixingSO2::AddSO2Settings(float value) { m_SO2Settings.push_back(value); } diff --git a/Modules/PhotoacousticsLib/test/mitkSpectralUnmixingTest.cpp b/Modules/PhotoacousticsLib/test/mitkSpectralUnmixingTest.cpp index 69e63fe1eb..e9ef3fd988 100644 --- a/Modules/PhotoacousticsLib/test/mitkSpectralUnmixingTest.cpp +++ b/Modules/PhotoacousticsLib/test/mitkSpectralUnmixingTest.cpp @@ -1,299 +1,328 @@ //*=================================================================== //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); //--> SO2 Settings FAILS 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]"; } void testEigenSUAlgorithm() { MITK_INFO << "START FILTER TEST ... "; // Set input image auto m_SpectralUnmixingFilter = mitk::pa::LinearSpectralUnmixingFilter::New(); 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); } //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/*mitk::pa::LinearSpectralUnmixingFilter::AlgortihmType::test*/}; 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 :)"; } void testVigraSUAlgorithm() { MITK_INFO << "START FILTER TEST ... "; // Set input image auto m_SpectralUnmixingFilter = mitk::pa::SpectralUnmixingFilterVigra::New(); m_SpectralUnmixingFilter->SetInput(inputImage); //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); //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); - m_sO2->SetInput(1, inputImage); - std::vector m_CorrectSO2Result = { 0.5, 0.5, 0.5, 0.5, 0 }; - std::vector SO2Settings = { 1, 1, 1, 49 }; // active settings (sum, fraction) are not tested (for false case) --> BUT ALL FALSE - for (int i = 0; i < SO2Settings.size(); ++i) - m_sO2->AddSO2Settings(SO2Settings[i]); - m_sO2->Update(); + 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"); - 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) + for (int k = 0; k < 5; ++k) { - auto Value = inputDataArray[Pixel]; + 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()); - myfile << "Output " << Pixel << ": " << "\n" << Value <<"\n"; - myfile << "Correct Result: " << "\n" << m_CorrectSO2Result[Pixel] << "\n"; + 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); + CPPUNIT_ASSERT(std::abs(Value - m_CorrectSO2Result[Pixel]) < threshold); + } } myfile.close(); MITK_INFO << "SO2 TEST SUCCESFULL :)"; } void tearDown() override { m_SpectralUnmixingFilter = nullptr; inputImage = nullptr; m_inputWavelengths.clear(); m_CorrectResult.clear(); MITK_INFO << "tearDown ... [DONE]"; } }; MITK_TEST_SUITE_REGISTRATION(mitkSpectralUnmixing)