diff --git a/Modules/PhotoacousticArtifacts/src/Algorithms/mitkPhotoacousticArtifact.cpp b/Modules/PhotoacousticArtifacts/src/Algorithms/mitkPhotoacousticArtifact.cpp index e6862b3be1..8316022b69 100644 --- a/Modules/PhotoacousticArtifacts/src/Algorithms/mitkPhotoacousticArtifact.cpp +++ b/Modules/PhotoacousticArtifacts/src/Algorithms/mitkPhotoacousticArtifact.cpp @@ -1,118 +1,224 @@ /*=================================================================== 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 #include -// step 1 in eigene klasse extrahieren +// Methode soll von einem eingeladenen 2D mitkImage die Hauptkomponenten berechnen, und die gewünschte +// Anzahl (1..n) zum Entrauschen abziehen. Mathematisch wird dabei die angegebene Anzahl der ersten +// Singulärwerte Null gesetzt und dann wierder zurückgerechnet. +template +static void SVDFilter2D(mitk::Image::Pointer inputImage, int svdNumber, mitk::Image::Pointer outputImage) +{ + //Step 1: inputImage in matrix m überführen, danach transponieren + unsigned int* dimensions = inputImage->GetDimensions(); //0=x, 1=y + unsigned int numberOfColums = dimensions[0]; + unsigned int numberOfRows = dimensions[1]; + //unsigned int zDimImmer1 = dimensions[2]; + + Eigen::MatrixXd m(numberOfRows, numberOfColums); + mitk::ImageReadAccessor readAccessor(inputImage, inputImage->GetVolumeData()); + const TPixel* data = (const TPixel*)readAccessor.GetData(); + for (unsigned int i = 0; i < numberOfColums; ++i) + { + for (unsigned int j = 0; j < numberOfRows; ++j) + { + m(j, i) = (double)data[i*numberOfRows + j]; + } + } + + //m.transpose(); + + //Step 2: eigenvectoren der thin Matrix U berechnen & die singularValues + // Compute Thin U, singularValues S and V + // Set desired number of components, which should be removed/substracted to zero -> (weight are zeroed for denoising) + // multiply all members (A=U*S*V(transpose)) back together again + Eigen::JacobiSVD svdSolver(m.transpose(), Eigen::ComputeThinU | Eigen::ComputeThinV); + Eigen::MatrixXd matrixU = svdSolver.matrixU(); + //MITK_INFO << endl << "Matrix U: " << endl << matrixU << endl; + Eigen::MatrixXd singularValues = svdSolver.singularValues(); + + // singularValues to matrix (m x m), if m < n + + // TODO: Fallunterscheidung: if (m > n) : SValues n x n; if (m < n) : SValues m x m !!!!!!!!!!!!!!!!!!!!!!!!!!! + // TODO: optimierung mit sparse matrix eigen + Eigen::MatrixXd singularValuesInMatrix = Eigen::MatrixXd::Zero(numberOfColums, numberOfColums); + auto SVdata = singularValues.data(); + for (unsigned int i = 0; i < numberOfColums; ++i) + { + singularValuesInMatrix(i, i) = (double)SVdata[i]; + } + + //MITK_INFO << endl << "Singulaerwerte: " << endl << singularValuesInMatrix << endl; + Eigen::MatrixXd matrixV = svdSolver.matrixV(); + //MITK_INFO << endl << "Matrix V: " << endl << matrixV << endl; + + //singularValuesInMatrix(0, 0) = 0.0; + //MITK_INFO << endl << "Reduzierte Singulaerwerte: " << endl << singularValuesInMatrix << endl; + + Eigen::MatrixXd denoisedMatrix = matrixU*singularValuesInMatrix*matrixV.transpose(); + + denoisedMatrix.transposeInPlace(); + + //MITK_INFO << endl << "Denoised Matrix: " << endl << denoisedMatrix << endl; -// fertige matrix übergeben, die zuvor aus einem bild überführt wurde + //Step 3: denoisedMatrix in outputImage überführen + unsigned int* dimOutputImage = new unsigned int[3]; + dimOutputImage[0] = dimensions[0]; + dimOutputImage[1] = dimensions[1]; + dimOutputImage[2] = dimensions[2]; + + mitk::PixelType pixelType = mitk::MakeScalarPixelType(); + outputImage->Initialize(pixelType, 3, dimOutputImage); + outputImage->SetSpacing(inputImage->GetGeometry()->GetSpacing()); + if (numberOfColums > numberOfRows) // not sure if this checks are necessary + { + for (unsigned int i = 0; i < numberOfRows; ++i) + { + outputImage->SetImportVolume(denoisedMatrix.col(i).data(), i); + } + } + else { + for (unsigned int i = 0; i < numberOfColums; ++i) + { + outputImage->SetImportVolume(denoisedMatrix.col(i).data(), i); + } + } +} + +// TODO: step 1 in eigene klasse extrahieren + +// TODO: fertige matrix übergeben, die zuvor aus einem bild überführt wurde +// Methode bekommt zuvor angepasste Matrix als input, es wird die ThinU Matrix berechnet und +// die Hauptkomponenten/Eigenvektoren werden zurück zum Bild konvertiert template -static void ComputeSVD(mitk::Image::Pointer inputImage, int svdNumber, mitk::Image::Pointer outputImage) +static void ComputeThinU(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 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; - if (numberOfSlices > vectorLength) - { - 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; + //MITK_INFO << endl << "Matrix U: " << endl << matrixU << endl; + //if (numberOfSlices > vectorLength) + //{ + // 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]; if (numberOfSlices > vectorLength) { dimOutputImage[2] = vectorLength; } else { dimOutputImage[2] = numberOfSlices; } mitk::PixelType pixelType = mitk::MakeScalarPixelType(); outputImage->Initialize(pixelType, 3, dimOutputImage); + outputImage->SetSpacing(inputImage->GetGeometry()->GetSpacing()); if (numberOfSlices > vectorLength) { 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); + //convert input from short/float to double type + //typedef itk::Image ShortImageType; + //typedef itk::Image DoubleImageType; + + //typedef itk::CastImageFilter CastFilterType; + //CastFilterType::Pointer castFilter = CastFilterType::New(); + //mitk::ImageToItk::Pointer filter = mitk::ImageToItk::New(); + //filter->SetInput(inputImage); + //filter->Update(); + //castFilter->SetInput(filter->GetOutput()); + //castFilter->Update(); + + //mitk::Image::Pointer newInputImage = mitk::Image::New(); + //auto castImage = castFilter->GetOutput(); + //newInputImage->InitializeByItk(castImage); + //newInputImage->SetVolume(castImage->GetBufferPointer()); + + // 3D thing + //ComputeThinU(newInputImage, 1, outputImage); + + // 2D thing + SVDFilter2D(inputImage, 1, outputImage); } diff --git a/Modules/PhotoacousticArtifacts/test/mitkPhotoacousticArtifactsTest.cpp b/Modules/PhotoacousticArtifacts/test/mitkPhotoacousticArtifactsTest.cpp index 8895496637..14b2565c7a 100644 --- a/Modules/PhotoacousticArtifacts/test/mitkPhotoacousticArtifactsTest.cpp +++ b/Modules/PhotoacousticArtifacts/test/mitkPhotoacousticArtifactsTest.cpp @@ -1,305 +1,412 @@ /*=================================================================== 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(Test2x2SimpleSVD); - MITK_TEST(Test2x3ComplexSVD); - MITK_TEST(Test7x3ComplexSVD); - MITK_TEST(Test7x4ComplexSVD); + //MITK_TEST(Test2x2SimpleSVD); + //MITK_TEST(Test2x3ComplexSVD); + //MITK_TEST(Test7x3ComplexSVD); + //MITK_TEST(Test7x4ComplexSVD); + MITK_TEST(Test3x4SVDFilter2D); + MITK_TEST(Test3x2SVDFilter2D); 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 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: 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 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 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; 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: 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[2] = -0.419598; //Fehler eingebaut; richtig: -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)); + CPPUNIT_ASSERT(mitk::Equal(*outputFilter, *m_ExpectedOutput, 0.000001, 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)); + CPPUNIT_ASSERT(mitk::Equal(*outputFilter, *m_ExpectedOutput, 0.000001, true)); + } + + // InputImage 3x4, wird danach transponiert -> 4x3 + void Test3x4SVDFilter2D() + { + // inputImage initialisieren + unsigned int* dimensions = new unsigned int[3]; + dimensions[0] = 4; // dim[0]=x + dimensions[1] = 3; // dim[1]=y + dimensions[2] = 1; // dim[2]=z + + double* dataIn = new double[12]; + dataIn[0] = 1.00446; + dataIn[1] = 5.00249; + dataIn[2] = 9.00723; + dataIn[3] = 2.00421; + dataIn[4] = 6.00007; + dataIn[5] = 10.0048; + dataIn[6] = 3.00397; + dataIn[7] = 6.99766; + dataIn[8] = 11.0024; + dataIn[9] = 4.00372; + dataIn[10] = 7.99524; + dataIn[11] = 11.9999; + 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] = 4; + dimOfOutputImage[1] = 3; + dimOfOutputImage[2] = 1; + + double* dataOut = new double[12]; + dataOut[0] = -1.12277; + dataOut[1] = -0.320719; + dataOut[2] = 0.479816; + dataOut[3] = -0.444206; + dataOut[4] = -0.126887; + dataOut[5] = 0.189832; + dataOut[6] = 0.234357; + dataOut[7] = 0.0669441; + dataOut[8] = -0.100153; + dataOut[9] = 0.91292; + dataOut[10] = 0.260776; + dataOut[11] = -0.390137; + m_ExpectedOutput->Initialize(pixelType, 3, dimOfOutputImage); + m_ExpectedOutput->SetImportVolume(dataOut, mitk::Image::ImportMemoryManagementType::CopyMemory); + + delete dataOut; + + MITK_INFO << "StartTest: 3x4 SVDFilter 2D"; + m_Filter->SetInput(m_Input); + m_Filter->Update(); + mitk::Image::Pointer outputFilter = m_Filter->GetOutput(); + CPPUNIT_ASSERT(mitk::Equal(*outputFilter, *m_ExpectedOutput, 0.001, true)); + } + + void Test3x2SVDFilter2D() + { + // inputImage initialisieren + unsigned int* dimensions = new unsigned int[3]; + dimensions[0] = 2; // dim[0]=x + dimensions[1] = 3; // dim[1]=y + dimensions[2] = 1; // dim[2]=z + + double* dataIn = new double[6]; + dataIn[0] = 7.0; + dataIn[1] = 17.0; + dataIn[2] = 7.0; + dataIn[3] = 7.0; + dataIn[4] = 17.0; + dataIn[5] = 77.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] = 2; + dimOfOutputImage[1] = 3; + dimOfOutputImage[2] = 1; + + double* dataOut = new double[6]; + dataOut[0] = 5.85346; + dataOut[1] = 14.2155; + dataOut[2] = -4.15134; + dataOut[3] = -0.85444; + dataOut[4] = -2.07507; + dataOut[5] = 0.605979; + m_ExpectedOutput->Initialize(pixelType, 3, dimOfOutputImage); + m_ExpectedOutput->SetImportVolume(dataOut, mitk::Image::ImportMemoryManagementType::CopyMemory); + + delete dataOut; + + MITK_INFO << "StartTest: 3x2 SVDFilter 2D"; + m_Filter->SetInput(m_Input); + m_Filter->Update(); + mitk::Image::Pointer outputFilter = m_Filter->GetOutput(); + CPPUNIT_ASSERT(mitk::Equal(*outputFilter, *m_ExpectedOutput, 0.0001, true)); } void tearDown() { m_Filter = nullptr; m_Input = nullptr; m_ExpectedOutput = nullptr; } }; MITK_TEST_SUITE_REGISTRATION(mitkPhotoacousticArtifacts) diff --git a/Plugins/org.mitk.gui.qt.photoacousticartifacts/src/internal/PhotoacousticArtifact.cpp b/Plugins/org.mitk.gui.qt.photoacousticartifacts/src/internal/PhotoacousticArtifact.cpp index 3920834c23..22d562547d 100644 --- a/Plugins/org.mitk.gui.qt.photoacousticartifacts/src/internal/PhotoacousticArtifact.cpp +++ b/Plugins/org.mitk.gui.qt.photoacousticartifacts/src/internal/PhotoacousticArtifact.cpp @@ -1,122 +1,122 @@ /*=================================================================== 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 "PhotoacousticArtifact.h" // QT includes (GUI) #include // Berry includes (selection service) #include #include // MITK includes (general) //#include #include #include // Includes for image casting between ITK and MITK #include #include const std::string PhotoacousticArtifact::VIEW_ID = "photoacousticproject.views.photoacoustic.artifacts"; void PhotoacousticArtifact::CreateQtPartControl(QWidget* parent) { m_Controls.setupUi(parent); - connect(m_Controls.applySVDPushButton, SIGNAL(clicked()), this, SLOT(ProcessSelectedImage())); + connect(m_Controls.computeNoisePushButton, SIGNAL(clicked()), this, SLOT(ProcessSelectedImage())); } void PhotoacousticArtifact::SetFocus() { - m_Controls.applySVDPushButton->setFocus(); + m_Controls.computeNoisePushButton->setFocus(); } void PhotoacousticArtifact::OnSelectionChanged(berry::IWorkbenchPart::Pointer, - const QList& dataNodes) + const QList& dataNodes) // TODO: dasselbe noch für "now select a raw 3D image" { for (const auto& dataNode : dataNodes) { if (dataNode.IsNotNull() && dynamic_cast(dataNode->GetData()) != nullptr) { m_Controls.selectImageLabel->setVisible(false); return; } } m_Controls.selectImageLabel->setVisible(true); } void PhotoacousticArtifact::ProcessSelectedImage() { auto selectedDataNodes = this->GetDataManagerSelection(); if (selectedDataNodes.empty()) return; auto firstSelectedDataNode = selectedDataNodes.front(); if (firstSelectedDataNode.IsNull()) { QMessageBox::information(nullptr, "Photoacoustic Artifact", "Please load and select an image before starting image processing."); return; } auto data = firstSelectedDataNode->GetData(); if (data != nullptr) { mitk::Image::Pointer image = dynamic_cast(data); // Something is selected and it contains data, but is it an image? if (image.IsNotNull()) { auto imageName = firstSelectedDataNode->GetName(); auto numberSVC = m_Controls.numberSVCSpinBox->value(); MITK_INFO << "Process image \"" << imageName << "\" ..."; - auto svc = mitk::PhotoacousticArtifact::New(); - svc->SetInput(image); + auto computeNoise = mitk::PhotoacousticArtifact::New(); + computeNoise->SetInput(image); //svc->SetNumberSVC(svc); - svc->Update(); + computeNoise->Update(); - mitk::Image::Pointer processedImage = svc->GetOutput(); + mitk::Image::Pointer noiseImage = computeNoise->GetOutput(); - if (processedImage.IsNull() || !processedImage->IsInitialized()) + if (noiseImage.IsNull() || !noiseImage->IsInitialized()) return; MITK_INFO << " done"; // Stuff the resulting image into a data node, set some properties, // and add it to the data storage, which will eventually display the // image in the application. - auto processedImageDataNode = mitk::DataNode::New(); - processedImageDataNode->SetData(processedImage); + auto noiseImageDataNode = mitk::DataNode::New(); + noiseImageDataNode->SetData(noiseImage); //QString name = QString("%1 (SVC: %2)").arg(imageName.c_str()).arg(svc); //processedImageDataNode->SetName(name.toStdString()); mitk::LevelWindow levelWindow; if (firstSelectedDataNode->GetLevelWindow(levelWindow)) - processedImageDataNode->SetLevelWindow(levelWindow); + noiseImageDataNode->SetLevelWindow(levelWindow); - this->GetDataStorage()->Add(processedImageDataNode); + this->GetDataStorage()->Add(noiseImageDataNode); } } } diff --git a/Plugins/org.mitk.gui.qt.photoacousticartifacts/src/internal/PhotoacousticArtifactViewControls.ui b/Plugins/org.mitk.gui.qt.photoacousticartifacts/src/internal/PhotoacousticArtifactViewControls.ui index 9497b098ea..366e73d01e 100644 --- a/Plugins/org.mitk.gui.qt.photoacousticartifacts/src/internal/PhotoacousticArtifactViewControls.ui +++ b/Plugins/org.mitk.gui.qt.photoacousticartifacts/src/internal/PhotoacousticArtifactViewControls.ui @@ -1,83 +1,517 @@ PhotoacousticArtifactViewControls 0 0 403 427 Photoacoustic Artifacts QLabel { color: rgb(255, 0, 0) } - Please select a RAW data image. + Please select an empty raw 3D data image. - - - - - 0 - - - 500 - - - - - - - Choose a number of SVC(s): - - - - + - + - Apply SVD + Compute noise + + + + + + + Qt::Vertical + + + + 20 + 60 + + + + + + + + + + + + + 255 + 0 + 0 + + + + + + + 255 + 0 + 0 + + + + + + + 255 + 127 + 127 + + + + + + + 255 + 63 + 63 + + + + + + + 127 + 0 + 0 + + + + + + + 170 + 0 + 0 + + + + + + + 255 + 0 + 0 + + + + + + + 255 + 255 + 255 + + + + + + + 0 + 0 + 0 + + + + + + + 255 + 255 + 255 + + + + + + + 255 + 0 + 0 + + + + + + + 0 + 0 + 0 + + + + + + + 255 + 127 + 127 + + + + + + + 255 + 255 + 220 + + + + + + + 0 + 0 + 0 + + + + + + + + + 255 + 0 + 0 + + + + + + + 255 + 0 + 0 + + + + + + + 255 + 127 + 127 + + + + + + + 255 + 63 + 63 + + + + + + + 127 + 0 + 0 + + + + + + + 170 + 0 + 0 + + + + + + + 255 + 0 + 0 + + + + + + + 255 + 255 + 255 + + + + + + + 0 + 0 + 0 + + + + + + + 255 + 255 + 255 + + + + + + + 255 + 0 + 0 + + + + + + + 0 + 0 + 0 + + + + + + + 255 + 127 + 127 + + + + + + + 255 + 255 + 220 + + + + + + + 0 + 0 + 0 + + + + + + + + + 127 + 0 + 0 + + + + + + + 255 + 0 + 0 + + + + + + + 255 + 127 + 127 + + + + + + + 255 + 63 + 63 + + + + + + + 127 + 0 + 0 + + + + + + + 170 + 0 + 0 + + + + + + + 120 + 120 + 120 + + + + + + + 255 + 255 + 255 + + + + + + + 127 + 0 + 0 + + + + + + + 255 + 0 + 0 + + + + + + + 255 + 0 + 0 + + + + + + + 0 + 0 + 0 + + + + + + + 255 + 0 + 0 + + + + + + + 255 + 255 + 220 + + + + + + + 0 + 0 + 0 + + + + + + + + Now select a raw 3D data image which should be denoised. + + + + + + + Choose a number of SVC(s): (how strong the image will be denoised) + + + + + + + 0 + + + 500 - Reconstruct image + Denoise image Qt::Vertical QSizePolicy::Expanding 20 220