diff --git a/Modules/PhotoacousticArtifacts/src/Algorithms/mitkPhotoacousticArtifact.cpp b/Modules/PhotoacousticArtifacts/src/Algorithms/mitkPhotoacousticArtifact.cpp index 5482f657b5..13ffd5b13a 100644 --- a/Modules/PhotoacousticArtifacts/src/Algorithms/mitkPhotoacousticArtifact.cpp +++ b/Modules/PhotoacousticArtifacts/src/Algorithms/mitkPhotoacousticArtifact.cpp @@ -1,221 +1,263 @@ /*=================================================================== 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 // 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 numberOfZeroSVCs, mitk::Image::Pointer outputImage) { //Step 1: inputImage in matrix m überführen, danach transponieren --> wird bei Step 2 gemacht 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]; } } //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 // transponieren Eigen::MatrixXd singularValuesInMatrix = Eigen::MatrixXd::Zero(numberOfColums, numberOfColums); auto SVdata = singularValues.data(); for (unsigned int i = numberOfZeroSVCs; 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; Eigen::MatrixXd denoisedMatrix = matrixU*singularValuesInMatrix*matrixV.transpose(); denoisedMatrix.transposeInPlace(); //MITK_INFO << endl << "Denoised Matrix: " << endl << denoisedMatrix << endl; //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 ComputeThinU(mitk::Image::Pointer inputImage, int svdNumber, mitk::Image::Pointer outputImage) +static void ComputeSVD3D(mitk::Image::Pointer inputImage, unsigned int numberOfSVCs, 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::JacobiSVD svdSolver(m, Eigen::ComputeThinU | Eigen::ComputeThinV); Eigen::MatrixXd matrixU = svdSolver.matrixU(); - + Eigen::MatrixXd singularValues = svdSolver.singularValues(); + Eigen::MatrixXd singularValuesInMatrix = Eigen::MatrixXd::Zero(numberOfSlices, numberOfSlices); + auto SVdata = singularValues.data(); + for (unsigned int i = 0; i < numberOfSVCs; ++i) + { + //unsigned int i = 2; + singularValuesInMatrix(i, i) = (double)SVdata[i]; + } + Eigen::MatrixXd matrixV = svdSolver.matrixV(); + Eigen::MatrixXd noiseMatrix3D = matrixU *singularValuesInMatrix *matrixV.transpose(); + //___________________________________________________________________________________________________________________________ //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; - } + dimOutputImage[2] = 1; // average bild + //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); - } - } + + // average bilden + Eigen::VectorXd averageImage = noiseMatrix3D.rowwise().mean(); + + outputImage->SetImportSlice(averageImage.data(), 0); + + //if (numberOfSlices > vectorLength) + //{ + // for (unsigned int i = 0; i < vectorLength; ++i) + // { + // outputImage->SetImportSlice(noiseMatrix3D.col(i).data(), i); + // } + //} + //else { + // for (unsigned int i = 0; i < numberOfSlices; ++i) + // { + // outputImage->SetImportSlice(noiseMatrix3D.col(i).data(), i); + // } + //} } +// methode zum denoisen 3DFall +//template +//static void SVDFilter3D(mitk::Image::Pointer inputImage, mitk::Image::Pointer outputImage, mitk::Image::Pointer noiseImage) +//{ +// mitk::Image::Pointer noiseImage = ComputeSVD3D->GetOutput(); +// //subtrahieren +// +// unsigned int* dimensions = inputImage->GetDimensions(); //0=x, 1=y, 2=z +// //unsigned int numberOfSlices = dimensions[2]; +// //unsigned int vectorLength = dimensions[0] * dimensions[1]; +// +// mitk::Image::Pointer sliceOfInputImage = inputImage->GetSliceData(0); +// mitk::Image::Pointer denoisedImage = sliceOfInputImage - noiseImage; +// +// unsigned int* dimOutputImage = new unsigned int[3]; +// dimOutputImage[0] = dimensions[0]; +// dimOutputImage[1] = dimensions[1]; +// dimOutputImage[2] = 1; +// mitk::PixelType pixelType = mitk::MakeScalarPixelType(); +// outputImage->Initialize(pixelType, 3, dimOutputImage); +// outputImage->SetSpacing(inputImage->GetGeometry()->GetSpacing()); +// +// outputImage->SetImportSlice(denoisedImage.data(), 0); +//} + mitk::PhotoacousticArtifact::PhotoacousticArtifact() : m_NumberOfSVC(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(); // TODO: switch case für jeden type //convert input from short/float to double type typedef itk::Image WhatEverImageType; 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); + ComputeSVD3D(newInputImage, m_NumberOfSVC, outputImage); // 2D thing - SVDFilter2D(newInputImage, m_NumberOfSVC, outputImage); + //SVDFilter2D(newInputImage, m_NumberOfSVC, outputImage); } 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 7c99f2f554..de7f1b482f 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,139 @@ /*=================================================================== 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 +#include +#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.computeNoise3DPushButton, SIGNAL(clicked()), this, SLOT(ProcessSelectedImage())); - connect(m_Controls.denoise2DPushButton, SIGNAL(clicked()), this, SLOT(ProcessSelectedImage())); + + m_Controls.selectComputedNoiseComboBox->SetDataStorage(this->GetDataStorage()); + m_Controls.selectComputedNoiseComboBox->SetPredicate(mitk::NodePredicateAnd::New( + mitk::TNodePredicateDataType::New(), + mitk::NodePredicateNot::New(mitk::NodePredicateProperty::New("helper object")))); + m_Controls.selectSignalImageComboBox->SetDataStorage(this->GetDataStorage()); + m_Controls.selectSignalImageComboBox->SetPredicate(mitk::NodePredicateAnd::New( + mitk::TNodePredicateDataType::New(), + mitk::NodePredicateNot::New(mitk::NodePredicateProperty::New("helper object")))); + + //connect(m_Controls.denoise2DPushButton, SIGNAL(clicked()), this, SLOT(ProcessSelectedImage())); + connect(m_Controls.computeNoise3DPushButton, SIGNAL(clicked()), this, SLOT(ProcessSelectedImage())); } void PhotoacousticArtifact::SetFocus() { m_Controls.computeNoise3DPushButton->setFocus(); } void PhotoacousticArtifact::OnSelectionChanged(berry::IWorkbenchPart::Pointer, 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.selectSignal2DImageLabel->setVisible(false); + m_Controls.selectEmpty3DImageLabel->setVisible(false); // selectSignal2DImageLabel return; } } - m_Controls.selectSignal2DImageLabel->setVisible(true); + m_Controls.selectEmpty3DImageLabel->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 numberOfSVC = m_Controls.numberOfSVC2DSpinBox->value(); + auto numberOfSVC3D = m_Controls.numberOfSVC3DSpinBox->value(); + auto numberOfSVC2D = m_Controls.numberOfSVC2DSpinBox->value(); MITK_INFO << "Process image \"" << imageName << "\" ..."; auto computeNoise = mitk::PhotoacousticArtifact::New(); computeNoise->SetInput(image); - computeNoise->SetNumberOfSVC(numberOfSVC); + computeNoise->SetNumberOfSVC(numberOfSVC3D); computeNoise->Update(); mitk::Image::Pointer noiseImage = computeNoise->GetOutput(); 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 noiseImageDataNode = mitk::DataNode::New(); noiseImageDataNode->SetData(noiseImage); - QString name = QString("%1_%2SVC").arg(imageName.c_str()).arg(numberOfSVC); + //QString name = QString("%1_%2SVC").arg(imageName.c_str()).arg(numberOfSVC2D); + //noiseImageDataNode->SetName(name.toStdString()); + QString name = QString("%1_averageNoise_%2SVC").arg(imageName.c_str()).arg(numberOfSVC3D); noiseImageDataNode->SetName(name.toStdString()); mitk::LevelWindow levelWindow; if (firstSelectedDataNode->GetLevelWindow(levelWindow)) noiseImageDataNode->SetLevelWindow(levelWindow); 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 9a219822e4..6e420dbca5 100644 --- a/Plugins/org.mitk.gui.qt.photoacousticartifacts/src/internal/PhotoacousticArtifactViewControls.ui +++ b/Plugins/org.mitk.gui.qt.photoacousticartifacts/src/internal/PhotoacousticArtifactViewControls.ui @@ -1,560 +1,181 @@ PhotoacousticArtifactViewControls 0 0 403 427 Photoacoustic Artifacts QLabel { color: rgb(255, 0, 0) } - [3D] (1): Please select an empty raw 3D data image. + [3D] (1): Please select an empty noise 3D data image. + + + + Choose a number of SVC(s): (how strong the signal image will be denoised later) + + + + + + + 0 + + + 500 + + + Compute noise Qt::Vertical 20 10 - - - - - - - 255 - 0 - 0 - - - - - - - 255 - 0 - 0 - - - - - - - 255 - 127 - 127 - - - - - - - 255 - 63 - 63 - - - - - - - 127 - 0 - 0 - - - - - - - 170 - 0 - 0 - - - - - - - 255 - 0 - 0 - - - - - - - 255 - 255 - 255 - - - - - - - 255 - 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 - - - - - - - 255 - 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 - - - - - - - 255 - 0 - 0 - - - - - - - 255 - 0 - 0 - - - - - - - 255 - 0 - 0 - - - - - - - 0 - 0 - 0 - - - - - - - 255 - 0 - 0 - - - - - - - 255 - 255 - 220 - - - - - - - 0 - 0 - 0 - - - - - - QLabel { color: rgb(255, 0, 0) } - [3D] (2): Now select a raw PA signal 3D data image which should be denoised. + [3D] (2): Now select the computed averageNoise image and a PA signal 3D data image which should be denoised. + + + true - + + + + + + + - Choose a number of SVC(s): (how strong the image will be denoised) + Denoise stack of images - - - 0 + + + Qt::Vertical - - 500 + + + 20 + 20 + - + - - - Denoise stack of images + + + Qt::Horizontal - + Qt::Vertical + + QSizePolicy::Expanding + 20 20 QLabel { color: rgb(255, 0, 0) } [2D]: Please select a raw PA signal 2D data image which should be denoised. Choose a number of SVC(s): (how strong the image will be denoised) Denoise image - + Qt::Vertical - - QSizePolicy::Expanding - 20 - 20 + 40 + + + QmitkDataStorageComboBox + QComboBox +
qmitkdatastoragecombobox.h
+ +