diff --git a/Modules/PhotoacousticArtifacts/src/Algorithms/mitkPhotoacousticArtifact.cpp b/Modules/PhotoacousticArtifacts/src/Algorithms/mitkPhotoacousticArtifact.cpp index 65b7e8d07c..dddaeb644a 100644 --- a/Modules/PhotoacousticArtifacts/src/Algorithms/mitkPhotoacousticArtifact.cpp +++ b/Modules/PhotoacousticArtifacts/src/Algorithms/mitkPhotoacousticArtifact.cpp @@ -1,84 +1,89 @@ /*=================================================================== 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 template static void ComputeSVD(mitk::Image::Pointer inputImage, int svdNumber, mitk::Image::Pointer outputImage) { //Step 1: inputImage in matrix m überführen - unsigned int* dimension = inputImage->GetDimensions(); //0=x, 1=y, 2=z - unsigned int numberOfSlices = dimension[2]; - unsigned int vectorLength = dimension[0] * dimension[1]; + 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::ComputeFullV); - Eigen::MatrixXd matrixV = svdSolver.matrixV(); - auto eigenVec1 = matrixV.col(0); - auto eigenVec2 = matrixV.col(1); - MITK_INFO << "Matrix V: " << endl << matrixV << endl; - MITK_INFO << "Eigenvector 1: " << endl << eigenVec1 << endl; - MITK_INFO << "Eigenvector 2: " << endl << eigenVec2 << endl; + Eigen::JacobiSVD svdSolver(m, Eigen::ComputeFullU); + Eigen::MatrixXd matrixU = svdSolver.matrixU(); + + MITK_INFO << endl << "Matrix U: " << endl << matrixU << endl; + for (unsigned int i = 0; i < vectorLength; ++i) + { + MITK_INFO << endl << "Eigenvector " << i + 1 << ": " << endl << matrixU.col(i) << 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; mitk::PixelType pixelType = mitk::MakeScalarPixelType(); - outputImage->Initialize(pixelType, 3, dimension); + outputImage->Initialize(pixelType, 3, dimOutputImage); for (unsigned int i = 0; i < numberOfSlices; ++i) { - outputImage->SetImportSlice(matrixV.col(i).data(), 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 1cd495a329..90074ef987 100644 --- a/Modules/PhotoacousticArtifacts/test/mitkPhotoacousticArtifactsTest.cpp +++ b/Modules/PhotoacousticArtifacts/test/mitkPhotoacousticArtifactsTest.cpp @@ -1,92 +1,136 @@ /*=================================================================== 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); 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(); + } - // inputImage initialisieren + void TestSimpleSVD() + { unsigned int* dimensions = new unsigned int[3]; - dimensions[0] = 1; - dimensions[1] = 2; - dimensions[2] = 2; + 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 - m_Input = mitk::Image::New(); - double* data = new double[4]; - data[0] = 2; - data[1] = -1; - data[2] = 2; - data[3] = 1; + // 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(data, mitk::Image::ImportMemoryManagementType::CopyMemory); + m_Input->SetImportVolume(dataIn, mitk::Image::ImportMemoryManagementType::CopyMemory); - delete data; + delete dataIn; // expectedOutputImage initialisieren - m_ExpectedOutput = mitk::Image::New(); double* dataOut = new double[4]; - dataOut[0] = 1 / sqrt(2); - dataOut[1] = 1 / sqrt(2); - dataOut[2] = -1 / sqrt(2); - dataOut[3] = 1 / sqrt(2); + 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"; + 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 TestSimpleSVD() { - MITK_INFO << "StartTest"; + void TestComplexSVD() + { + 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; + + 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"; 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_ExpectedOutput = nullptr; m_Filter = nullptr; m_Input = nullptr; + m_ExpectedOutput = nullptr; } }; MITK_TEST_SUITE_REGISTRATION(mitkPhotoacousticArtifacts)