diff --git a/Modules/PhotoacousticArtifacts/include/mitkPhotoacousticArtifact.h b/Modules/PhotoacousticArtifacts/include/mitkPhotoacousticArtifact.h index ce20b14908..9838a7d0c7 100644 --- a/Modules/PhotoacousticArtifacts/include/mitkPhotoacousticArtifact.h +++ b/Modules/PhotoacousticArtifacts/include/mitkPhotoacousticArtifact.h @@ -1,53 +1,53 @@ /*=================================================================== 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 mitkPhotoacousticArtifact_h #define mitkPhotoacousticArtifact_h #include #include #include #include #include namespace mitk { class MITKPHOTOACOUSTICARTIFACTS_EXPORT PhotoacousticArtifact final : public mitk::ImageToImageFilter { public: // All classes that derive from an ITK-based MITK class need at least the // following two macros. Make sure you don't declare the constructor public // to force clients of your class to follow the ITK convention for // instantiating classes via the static New() method. mitkClassMacro(PhotoacousticArtifact, mitk::ImageToImageFilter) itkFactorylessNewMacro(Self) - itkSetMacro(SVD, int) - itkGetMacro(SVD, int) + itkSetMacro(NumberOfSVC, int) + itkGetMacro(NumberOfSVC, int) protected: private: PhotoacousticArtifact(); ~PhotoacousticArtifact(); void GenerateData() override; - int m_SVD; + int m_NumberOfSVC; }; } #endif mitkPhotoacousticArtifact_h diff --git a/Modules/PhotoacousticArtifacts/src/Algorithms/mitkPhotoacousticArtifact.cpp b/Modules/PhotoacousticArtifacts/src/Algorithms/mitkPhotoacousticArtifact.cpp index 7e3024dd9a..5482f657b5 100644 --- a/Modules/PhotoacousticArtifacts/src/Algorithms/mitkPhotoacousticArtifact.cpp +++ b/Modules/PhotoacousticArtifacts/src/Algorithms/mitkPhotoacousticArtifact.cpp @@ -1,223 +1,221 @@ /*=================================================================== 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 svdNumber, mitk::Image::Pointer outputImage) // svdNumber = numberOfZeroSVCs +static void SVDFilter2D(mitk::Image::Pointer inputImage, int numberOfZeroSVCs, mitk::Image::Pointer outputImage) { - //Step 1: inputImage in matrix m überführen, danach transponieren + //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]; } } - //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: 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 = svdNumber; i < numberOfColums; ++i) // i = numberOfZeroSVCs + 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) { //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; //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) + : 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 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); // 2D thing - SVDFilter2D(newInputImage, 20, 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 22d562547d..7c99f2f554 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.computeNoisePushButton, SIGNAL(clicked()), this, SLOT(ProcessSelectedImage())); + //connect(m_Controls.computeNoise3DPushButton, SIGNAL(clicked()), this, SLOT(ProcessSelectedImage())); + connect(m_Controls.denoise2DPushButton, SIGNAL(clicked()), this, SLOT(ProcessSelectedImage())); } void PhotoacousticArtifact::SetFocus() { - m_Controls.computeNoisePushButton->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.selectImageLabel->setVisible(false); + m_Controls.selectSignal2DImageLabel->setVisible(false); return; } } - m_Controls.selectImageLabel->setVisible(true); + m_Controls.selectSignal2DImageLabel->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(); + auto numberOfSVC = m_Controls.numberOfSVC2DSpinBox->value(); MITK_INFO << "Process image \"" << imageName << "\" ..."; auto computeNoise = mitk::PhotoacousticArtifact::New(); computeNoise->SetInput(image); - //svc->SetNumberSVC(svc); - + computeNoise->SetNumberOfSVC(numberOfSVC); 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 (SVC: %2)").arg(imageName.c_str()).arg(svc); - //processedImageDataNode->SetName(name.toStdString()); + QString name = QString("%1_%2SVC").arg(imageName.c_str()).arg(numberOfSVC); + 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 366e73d01e..9a219822e4 100644 --- a/Plugins/org.mitk.gui.qt.photoacousticartifacts/src/internal/PhotoacousticArtifactViewControls.ui +++ b/Plugins/org.mitk.gui.qt.photoacousticartifacts/src/internal/PhotoacousticArtifactViewControls.ui @@ -1,517 +1,560 @@ PhotoacousticArtifactViewControls 0 0 403 427 Photoacoustic Artifacts - + QLabel { color: rgb(255, 0, 0) } - Please select an empty raw 3D data image. + [3D] (1): Please select an empty raw 3D data image. - + Compute noise - + Qt::Vertical 20 - 60 + 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 - 0 + 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 - 0 + 255 0 0 255 255 255 255 0 0 0 0 0 255 127 127 255 255 220 0 0 0 - 127 + 255 0 0 255 0 0 255 127 127 255 63 63 127 0 0 170 0 0 - 120 - 120 - 120 + 255 + 0 + 0 255 255 255 - 127 + 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) } + - Now select a raw 3D data image which should be denoised. + [3D] (2): Now select a raw PA signal 3D data image which should be denoised. - + Choose a number of SVC(s): (how strong the image will be denoised) - + 0 500 - + + + Denoise stack of images + + + + + + + Qt::Vertical + + + + 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 - 220 + 20