diff --git a/Modules/PhotoacousticArtifacts/src/Algorithms/mitkPhotoacousticArtifact.cpp b/Modules/PhotoacousticArtifacts/src/Algorithms/mitkPhotoacousticArtifact.cpp index dddaeb644a..e6862b3be1 100644 --- a/Modules/PhotoacousticArtifacts/src/Algorithms/mitkPhotoacousticArtifact.cpp +++ b/Modules/PhotoacousticArtifacts/src/Algorithms/mitkPhotoacousticArtifact.cpp @@ -1,89 +1,118 @@ /*=================================================================== 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 +// step 1 in eigene klasse extrahieren + +// fertige matrix übergeben, die zuvor aus einem bild überführt wurde template static void ComputeSVD(mitk::Image::Pointer inputImage, int svdNumber, mitk::Image::Pointer outputImage) { //Step 1: inputImage in matrix m überführen unsigned int* dimensions = inputImage->GetDimensions(); //0=x, 1=y, 2=z unsigned int numberOfSlices = dimensions[2]; unsigned int vectorLength = dimensions[0] * dimensions[1]; Eigen::MatrixXd m(vectorLength, numberOfSlices); for (unsigned int i = 0; i < numberOfSlices; ++i) { mitk::ImageReadAccessor readAccessor(inputImage, inputImage->GetSliceData(i)); const TPixel* data = (const TPixel*)readAccessor.GetData(); for (unsigned int j = 0; j < vectorLength; ++j) { m(j, i) = (double)data[j]; } } - //Step 2: eigenvectoren der matrix berechnen - Eigen::JacobiSVD svdSolver(m, Eigen::ComputeFullU); + //Step 2: eigenvectoren der thin Matrix U berechnen & die singularValues + Eigen::JacobiSVD svdSolver(m, Eigen::ComputeThinU); Eigen::MatrixXd matrixU = svdSolver.matrixU(); MITK_INFO << endl << "Matrix U: " << endl << matrixU << endl; - for (unsigned int i = 0; i < vectorLength; ++i) + if (numberOfSlices > vectorLength) { - MITK_INFO << endl << "Eigenvector " << i + 1 << ": " << endl << matrixU.col(i) << endl; + for (unsigned int i = 0; i < vectorLength; ++i) + { + MITK_INFO << endl << "Eigenvector " << i + 1 << ": " << endl << matrixU.col(i) << endl; + } + } + else { + for (unsigned int i = 0; i < numberOfSlices; ++i) + { + MITK_INFO << endl << "Eigenvector " << i + 1 << ": " << endl << matrixU.col(i) << endl; + } } + MITK_INFO << endl << "Singulaerwerte: " << endl << svdSolver.singularValues() << endl; + //Step 3: eigenvectors in outputImage überführen unsigned int* dimOutputImage = new unsigned int[3]; dimOutputImage[0] = dimensions[0]; dimOutputImage[1] = dimensions[1]; - dimOutputImage[2] = vectorLength; + if (numberOfSlices > vectorLength) + { + dimOutputImage[2] = vectorLength; + } + else { + dimOutputImage[2] = numberOfSlices; + } mitk::PixelType pixelType = mitk::MakeScalarPixelType(); outputImage->Initialize(pixelType, 3, dimOutputImage); - for (unsigned int i = 0; i < numberOfSlices; ++i) + if (numberOfSlices > vectorLength) { - outputImage->SetImportSlice(matrixU.col(i).data(), i); + for (unsigned int i = 0; i < vectorLength; ++i) + { + outputImage->SetImportSlice(matrixU.col(i).data(), i); + } + } + else { + for (unsigned int i = 0; i < numberOfSlices; ++i) + { + outputImage->SetImportSlice(matrixU.col(i).data(), i); + } } } mitk::PhotoacousticArtifact::PhotoacousticArtifact() : m_SVD(1) { this->SetPrimaryOutput(mitk::Image::New()); MITK_INFO << "Hello."; this->SetNumberOfRequiredInputs(1); this->SetNumberOfRequiredOutputs(1); } mitk::PhotoacousticArtifact::~PhotoacousticArtifact() { } void mitk::PhotoacousticArtifact::GenerateData() { mitk::Image::Pointer inputImage = this->GetInput(0); mitk::Image::Pointer outputImage = this->GetOutput(); ComputeSVD(inputImage, 1, outputImage); } diff --git a/Modules/PhotoacousticArtifacts/test/mitkPhotoacousticArtifactsTest.cpp b/Modules/PhotoacousticArtifacts/test/mitkPhotoacousticArtifactsTest.cpp index 90074ef987..8895496637 100644 --- a/Modules/PhotoacousticArtifacts/test/mitkPhotoacousticArtifactsTest.cpp +++ b/Modules/PhotoacousticArtifacts/test/mitkPhotoacousticArtifactsTest.cpp @@ -1,136 +1,305 @@ /*=================================================================== 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 class mitkPhotoacousticArtifactsTestSuite : public mitk::TestFixture { CPPUNIT_TEST_SUITE(mitkPhotoacousticArtifactsTestSuite); - MITK_TEST(TestSimpleSVD); - MITK_TEST(TestComplexSVD); + MITK_TEST(Test2x2SimpleSVD); + MITK_TEST(Test2x3ComplexSVD); + MITK_TEST(Test7x3ComplexSVD); + MITK_TEST(Test7x4ComplexSVD); CPPUNIT_TEST_SUITE_END(); private: // alle Parameter: inputImage, expectedOutputImage... mitk::PhotoacousticArtifact::Pointer m_Filter; mitk::Image::Pointer m_Input; mitk::Image::Pointer m_ExpectedOutput; public: void setUp() { m_Filter = mitk::PhotoacousticArtifact::New(); m_Input = mitk::Image::New(); m_ExpectedOutput = mitk::Image::New(); } - void TestSimpleSVD() + void Test2x2SimpleSVD() { unsigned int* dimensions = new unsigned int[3]; dimensions[0] = 1; // dim[0]=x, always 1 in my case dimensions[1] = 2; // dim[1]=y, number of pixel in one imageSlice = length of my vector dimensions[2] = 2; // dim[2]=z, number of slices // inputImage initialisieren double* dataIn = new double[4]; dataIn[0] = 2; dataIn[1] = -1; dataIn[2] = 2; dataIn[3] = 1; mitk::PixelType pixelType = mitk::MakeScalarPixelType(); m_Input->Initialize(pixelType, 3, dimensions); m_Input->SetImportVolume(dataIn, mitk::Image::ImportMemoryManagementType::CopyMemory); delete dataIn; // expectedOutputImage initialisieren double* dataOut = new double[4]; dataOut[0] = 1.0; dataOut[1] = 0.0; dataOut[2] = 0.0; dataOut[3] = 1.0; m_ExpectedOutput->Initialize(pixelType, 3, dimensions); m_ExpectedOutput->SetImportVolume(dataOut, mitk::Image::ImportMemoryManagementType::CopyMemory); delete dataOut; - MITK_INFO << "StartTest: SimpleSVD"; + MITK_INFO << "StartTest: 2x2 SimpleSVD"; m_Filter->SetInput(m_Input); m_Filter->Update(); mitk::Image::Pointer outputFilter = m_Filter->GetOutput(); CPPUNIT_ASSERT(mitk::Equal(*outputFilter, *m_ExpectedOutput, mitk::eps, true)); } - void TestComplexSVD() + void Test2x3ComplexSVD() { + // inputImage initialisieren unsigned int* dimensions = new unsigned int[3]; dimensions[0] = 1; // dim[0]=x dimensions[1] = 2; // dim[1]=y dimensions[2] = 3; // dim[2]=z - // inputImage initialisieren double* dataIn = new double[6]; dataIn[0] = 0.36; dataIn[1] = 0.48; dataIn[2] = 1.60; dataIn[3] = -1.20; dataIn[4] = 0.48; dataIn[5] = 0.64; mitk::PixelType pixelType = mitk::MakeScalarPixelType(); m_Input->Initialize(pixelType, 3, dimensions); m_Input->SetImportVolume(dataIn, mitk::Image::ImportMemoryManagementType::CopyMemory); delete dataIn; + // expectedOutputImage initialisieren unsigned int* dimOfOutputImage = new unsigned int[3]; dimOfOutputImage[0] = 1; dimOfOutputImage[1] = 2; dimOfOutputImage[2] = 2; - // expectedOutputImage initialisieren + double* dataOut = new double[4]; dataOut[0] = -0.8; dataOut[1] = 0.6; dataOut[2] = 0.6; dataOut[3] = 0.8; m_ExpectedOutput->Initialize(pixelType, 3, dimOfOutputImage); m_ExpectedOutput->SetImportVolume(dataOut, mitk::Image::ImportMemoryManagementType::CopyMemory); delete dataOut; - MITK_INFO << "StartTest: ComplexSVD"; + MITK_INFO << "StartTest: 2x3 ComplexSVD"; + m_Filter->SetInput(m_Input); + m_Filter->Update(); + mitk::Image::Pointer outputFilter = m_Filter->GetOutput(); + CPPUNIT_ASSERT(mitk::Equal(*outputFilter, *m_ExpectedOutput, mitk::eps, true)); + } + + void Test7x3ComplexSVD() + { + // inputImage initialisieren + unsigned int* dimensions = new unsigned int[3]; + dimensions[0] = 1; // dim[0]=x + dimensions[1] = 7; // dim[1]=y + dimensions[2] = 3; // dim[2]=z + + double* dataIn = new double[21]; + dataIn[0] = 1.0; + dataIn[1] = -2.0; + dataIn[2] = 3.0; + dataIn[3] = -2.0; + dataIn[4] = 5.0; + dataIn[5] = -3.0; + dataIn[6] = 2.0; + dataIn[7] = 0.0; + dataIn[8] = 1.0; + dataIn[9] = 0.0; + dataIn[10] = 2.0; + dataIn[11] = 1.0; + dataIn[12] = 2.0; + dataIn[13] = 0.0; + dataIn[14] = -1.0; + dataIn[15] = 4.0; + dataIn[16] = -4.0; + dataIn[17] = 5.0; + dataIn[18] = 0.0; + dataIn[19] = 6.0; + dataIn[20] = 3.0; + mitk::PixelType pixelType = mitk::MakeScalarPixelType(); + m_Input->Initialize(pixelType, 3, dimensions); + m_Input->SetImportVolume(dataIn, mitk::Image::ImportMemoryManagementType::CopyMemory); + + delete dataIn; + + // expectedOutputImage initialisieren + unsigned int* dimOfOutputImage = new unsigned int[3]; + dimOfOutputImage[0] = 1; + dimOfOutputImage[1] = 7; + dimOfOutputImage[2] = 3; + + double* dataOut = new double[21]; + dataOut[0] = -0.115414; + dataOut[1] = 0.3953; + dataOut[2] = -0.419598; + dataOut[3] = 0.486417; + dataOut[4] = -0.192536; + dataOut[5] = 0.601832; + dataOut[6] = 0.135945; + dataOut[7] = 0.0752405; + dataOut[8] = 0.0460629; + dataOut[9] = 0.148534; + dataOut[10] = 0.165419; + dataOut[11] = 0.804309; + dataOut[12] = 0.0901785; + dataOut[13] = 0.53642; + dataOut[14] = 0.109204; // ab dem dritten Eigenvector kommen vorzeicheninvertierte Werte raus + dataOut[15] = -0.00205604; // damit der Test nicht fehlschläg sind die Vorzeichen 14-20 invertiert + dataOut[16] = 0.504987; + dataOut[17] = 0.391671; + dataOut[18] = 0.230245; + dataOut[19] = 0.282467; + dataOut[20] = -0.668469; + m_ExpectedOutput->Initialize(pixelType, 3, dimOfOutputImage); + m_ExpectedOutput->SetImportVolume(dataOut, mitk::Image::ImportMemoryManagementType::CopyMemory); + + delete dataOut; + + MITK_INFO << "StartTest: 7x3 ComplexSVD"; + m_Filter->SetInput(m_Input); + m_Filter->Update(); + mitk::Image::Pointer outputFilter = m_Filter->GetOutput(); + CPPUNIT_ASSERT(mitk::Equal(*outputFilter, *m_ExpectedOutput, mitk::eps, true)); + } + + void Test7x4ComplexSVD() + { + // inputImage initialisieren + unsigned int* dimensions = new unsigned int[3]; + dimensions[0] = 1; // dim[0]=x + dimensions[1] = 7; // dim[1]=y + dimensions[2] = 4; // dim[2]=z + + double* dataIn = new double[28]; + dataIn[0] = 3.0; + dataIn[1] = 1.0; + dataIn[2] = 0.0; + dataIn[3] = 8.0; + dataIn[4] = -1.0; + dataIn[5] = -4.0; + dataIn[6] = 7.0; + dataIn[7] = 4.0; + dataIn[8] = -2.0; + dataIn[9] = -5.0; + dataIn[10] = -6.0; + dataIn[11] = -2.0; + dataIn[12] = 7.0; + dataIn[13] = 7.0; + dataIn[14] = 7.0; + dataIn[15] = 4.0; + dataIn[16] = -5.0; + dataIn[17] = 3.0; + dataIn[18] = -3.0; + dataIn[19] = 0.0; + dataIn[20] = -7.0; + dataIn[21] = 7.0; + dataIn[22] = 8.0; + dataIn[23] = 9.0; + dataIn[24] = 4.0; + dataIn[25] = 2.0; + dataIn[26] = 3.0; + dataIn[27] = 1.0; + mitk::PixelType pixelType = mitk::MakeScalarPixelType(); + m_Input->Initialize(pixelType, 3, dimensions); + m_Input->SetImportVolume(dataIn, mitk::Image::ImportMemoryManagementType::CopyMemory); + + delete dataIn; + + // expectedOutputImage initialisieren + unsigned int* dimOfOutputImage = new unsigned int[3]; + dimOfOutputImage[0] = 1; + dimOfOutputImage[1] = 7; + dimOfOutputImage[2] = 4; + + double* dataOut = new double[28]; + dataOut[0] = 0.430212; + dataOut[1] = 0.520996; + dataOut[2] = 0.458767; + dataOut[3] = 0.556835; + dataOut[4] = 0.0676849; + dataOut[5] = -0.105548; + dataOut[6] = -0.0850128; + dataOut[7] = 0.481268; + dataOut[8] = 0.0782124; + dataOut[9] = -0.0592848; + dataOut[10] = -0.18708; + dataOut[11] = -0.0795045; + dataOut[12] = 0.503429; + dataOut[13] = 0.681167; + dataOut[14] = 0.455248; // ab dem dritten Eigenvector kommen vorzeicheninvertierte Werte raus + dataOut[15] = 0.190691; // damit der Test nicht fehlschläg sind die Vorzeichen 14-27 invertiert + dataOut[16] = -0.496905; + dataOut[17] = -0.144413; + dataOut[18] = -0.233733; + dataOut[19] = 0.224215; + dataOut[20] = -0.619446; + dataOut[21] = 0.130432; + dataOut[22] = -0.172334; + dataOut[23] = -0.565834; + dataOut[24] = 0.523608; + dataOut[25] = -0.261511; + dataOut[26] = -0.437446; + dataOut[27] = 0.314973; + m_ExpectedOutput->Initialize(pixelType, 3, dimOfOutputImage); + m_ExpectedOutput->SetImportVolume(dataOut, mitk::Image::ImportMemoryManagementType::CopyMemory); + + delete dataOut; + + MITK_INFO << "StartTest: 7x4 ComplexSVD"; m_Filter->SetInput(m_Input); m_Filter->Update(); mitk::Image::Pointer outputFilter = m_Filter->GetOutput(); CPPUNIT_ASSERT(mitk::Equal(*outputFilter, *m_ExpectedOutput, mitk::eps, true)); } void tearDown() { m_Filter = nullptr; m_Input = nullptr; m_ExpectedOutput = nullptr; } }; MITK_TEST_SUITE_REGISTRATION(mitkPhotoacousticArtifacts)