diff --git a/Modules/PhotoacousticsAlgorithms/include/mitkPhotoacousticMotionCorrectionFilter.h b/Modules/PhotoacousticsAlgorithms/include/mitkPhotoacousticMotionCorrectionFilter.h index 83bf01b5cb..2a98bbed0e 100644 --- a/Modules/PhotoacousticsAlgorithms/include/mitkPhotoacousticMotionCorrectionFilter.h +++ b/Modules/PhotoacousticsAlgorithms/include/mitkPhotoacousticMotionCorrectionFilter.h @@ -1,90 +1,113 @@ /*=================================================================== 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 _MITKPHOTOACOUSTICSMOTIONCORRECTIONFILTER_H_ #define _MITKPHOTOACOUSTICSMOTIONCORRECTIONFILTER_H_ #include "mitkImageToImageFilter.h" #include #include "opencv2/imgproc.hpp" // TODO: Find out why build fails with this option or replace with something else /* #include "opencv2/opencv.hpp" */ #include "opencv2/video/tracking.hpp" #include "itkOpenCVImageBridge.h" #include #include "mitkImageCast.h" #include #include namespace mitk { /*! * \brief Class implementing a mitk::ImageToImageFilter for PAUS motion correction. * * TODO: Write down all the parameters needed. * The filter takes a stack of PA and US images. It then computes the optical flow * within the US image and compensates the PA images for the flow. Afterwards it * returns the stack of PA images. */ class MITKPHOTOACOUSTICSALGORITHMS_EXPORT PhotoacousticMotionCorrectionFilter : public ImageToImageFilter { public: mitkClassMacro(PhotoacousticMotionCorrectionFilter, ImageToImageFilter); itkFactorylessNewMacro(Self); /* itkCloneMacro(Self); */ /* void SetInput(Image::Pointer paImage, Image::Pointer usImage); */ - // TODO: implement setters and getters for all the parameters + void SetBatchSize(unsigned int); + unsigned int GetBatchSize(); + void SetPyrScale(double); + double GetPyrScale(); + void SetLevels(unsigned int); + unsigned int GetLevels(); + void SetWindowSize(unsigned int); + unsigned int GetWindowSize(); + void SetIterations(unsigned int); + unsigned int GetIterations(); + void SetPolyN(unsigned int); + unsigned int GetPolyN(); + void SetPolySigma(double); + double GetPolySigma(); + void SetFlags(unsigned int); + unsigned int GetFlags(); + + // Wrapper for SetInput, GetInput and GetOutput + void SetPaInput(const mitk::Image::Pointer); + mitk::Image::Pointer GetPaInput(); + void SetUsInput(const mitk::Image::Pointer); + mitk::Image::Pointer GetUsInput(); + mitk::Image::Pointer GetPaOutput(); + mitk::Image::Pointer GetUsOutput(); protected: PhotoacousticMotionCorrectionFilter(); ~PhotoacousticMotionCorrectionFilter() override; void GenerateData() override; void CheckInput(mitk::Image::Pointer paImage, mitk::Image::Pointer usImage); void InitializeOutput(mitk::Image::Pointer paInput, mitk::Image::Pointer usInput, mitk::Image::Pointer paOutput, mitk::Image::Pointer usOutput); void SetOutputData(mitk::Image::Pointer input, mitk::Image::Pointer output); void PerformCorrection(mitk::Image::Pointer paInput, mitk::Image::Pointer usInput, mitk::Image::Pointer paOutput, mitk::Image::Pointer usOutput); cv::Mat GetMatrix(const mitk::Image::Pointer input, unsigned int i); void EnterMatrixInPosition(cv::Mat mat, mitk::Image::Pointer output, unsigned int i); //##Description //## @brief Time when Header was last initialized /* itk::TimeStamp m_TimeOfHeaderInitialization; */ private: // InputData /* mitk::Image::Pointer m_paImage, m_usImage, m_paCompensated, m_usCompensated; */ // Parameters double m_pyr_scale, m_poly_sigma; unsigned int m_levels, m_winsize, m_iterations, m_poly_n, m_flags, m_batch; /* // Stuff that OpenCV needs */ cv::UMat m_UsRef, m_PaMat, m_UsMat, m_Flow; cv::Mat m_PaRes, m_UsRes, m_PaMatC, m_UsMatC; // middle step conversion from MITK::image to cv::Mat // TODO: Note that there is always a float conversion inbetween itk::Image::Pointer m_itkPaImage, m_itkUsImage; mitk::OpenCVToMitkImageFilter::Pointer m_OpenCVToImageFilter = mitk::OpenCVToMitkImageFilter::New(); mitk::ImageToOpenCVImageFilter::Pointer m_ImageToOpenCVFilter = mitk::ImageToOpenCVImageFilter::New(); }; } #endif diff --git a/Modules/PhotoacousticsAlgorithms/source/mitkPhotoacousticMotionCorrectionFilter.cpp b/Modules/PhotoacousticsAlgorithms/source/mitkPhotoacousticMotionCorrectionFilter.cpp index 912ac00f8e..6cad15d997 100644 --- a/Modules/PhotoacousticsAlgorithms/source/mitkPhotoacousticMotionCorrectionFilter.cpp +++ b/Modules/PhotoacousticsAlgorithms/source/mitkPhotoacousticMotionCorrectionFilter.cpp @@ -1,186 +1,272 @@ /*=================================================================== 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 "./mitkPhotoacousticMotionCorrectionFilter.h" #include mitk::PhotoacousticMotionCorrectionFilter:: PhotoacousticMotionCorrectionFilter() { // Set the defaults for the OpticalFlowFarneback algorithm // The values are taken directly out of Thomas's // US-CV-Based-Optical-Flow-Carotis.ipyn m_batch = 5; m_pyr_scale = 0.5; m_levels = 1; m_winsize = 40; m_iterations = 2; m_poly_n = 7; m_poly_sigma = 1.5; m_flags = 0; this->SetNumberOfIndexedInputs(2); this->SetNumberOfIndexedOutputs(2); this->SetNthOutput(0, mitk::Image::New()); this->SetNthOutput(1, mitk::Image::New()); } mitk::PhotoacousticMotionCorrectionFilter:: ~PhotoacousticMotionCorrectionFilter() {} -// TODO: Find out how to throw the right errors +// Setters and Getters +// TODO: refactor member variables +void mitk::PhotoacousticMotionCorrectionFilter::SetBatchSize(unsigned int batch) { + m_batch = batch; +} + +unsigned int mitk::PhotoacousticMotionCorrectionFilter::GetBatchSize() { + return m_batch; +} + +void mitk::PhotoacousticMotionCorrectionFilter::SetPyrScale(double pyr_scale) { + m_pyr_scale = pyr_scale; +} + +double mitk::PhotoacousticMotionCorrectionFilter::GetPyrScale() { + return m_pyr_scale; +} + +void mitk::PhotoacousticMotionCorrectionFilter::SetLevels(unsigned int levels) { + m_levels = levels; +} + +unsigned int mitk::PhotoacousticMotionCorrectionFilter::GetLevels() { + return m_levels; +} + +void mitk::PhotoacousticMotionCorrectionFilter::SetWindowSize(unsigned int winsize) { + m_winsize = winsize; +} + +unsigned int mitk::PhotoacousticMotionCorrectionFilter::GetWindowSize() { + return m_winsize; +} + +void mitk::PhotoacousticMotionCorrectionFilter::SetIterations(unsigned int iterations) { + m_iterations = iterations; +} + +unsigned int mitk::PhotoacousticMotionCorrectionFilter::GetIterations() { + return m_iterations; +} + +void mitk::PhotoacousticMotionCorrectionFilter::SetPolyN(unsigned int poly_n) { + m_poly_n = poly_n; +} + +unsigned int mitk::PhotoacousticMotionCorrectionFilter::GetPolyN() { + return m_poly_n; +} + +void mitk::PhotoacousticMotionCorrectionFilter::SetPolySigma(double poly_sigma) { + m_poly_sigma = poly_sigma; +} + +double mitk::PhotoacousticMotionCorrectionFilter::GetPolySigma() { + return m_poly_sigma; +} + +void mitk::PhotoacousticMotionCorrectionFilter::SetFlags(unsigned int flags) { + m_flags = flags; +} + +unsigned int mitk::PhotoacousticMotionCorrectionFilter::GetFlags() { + return m_flags; +} + +void mitk::PhotoacousticMotionCorrectionFilter::SetPaInput(const mitk::Image::Pointer input) { + this->SetInput(0, input); +} + +mitk::Image::Pointer mitk::PhotoacousticMotionCorrectionFilter::GetPaInput() { + return this->GetInput(0); +} + +void mitk::PhotoacousticMotionCorrectionFilter::SetUsInput(const mitk::Image::Pointer input) { + this->SetInput(1, input); +} + +mitk::Image::Pointer mitk::PhotoacousticMotionCorrectionFilter::GetUsInput() { + return this->GetInput(1); +} + +mitk::Image::Pointer mitk::PhotoacousticMotionCorrectionFilter::GetPaOutput() { + return this->GetOutput(0); +} + +mitk::Image::Pointer mitk::PhotoacousticMotionCorrectionFilter::GetUsOutput() { + return this->GetOutput(1); +} + void mitk::PhotoacousticMotionCorrectionFilter::CheckInput(mitk::Image::Pointer paImage, mitk::Image::Pointer usImage) { // Check that we actually got some images if (!paImage || !usImage) { - // TODO: Throw some error here MITK_INFO << "We did not get two images!"; throw std::invalid_argument("One of the images was NULL."); } // Check that the image dimensions are the same if (paImage->GetDimension() != 3 || usImage->GetDimension() != 3) { MITK_INFO << "Mismatching image dimensions detected in the motion " "compensation filter."; - // TODO: Throw some error here throw std::invalid_argument("Both images must have dimension 3."); } // Check that each dimension has the same size for (unsigned int i = 0; i < paImage->GetDimension(); i++) { if (paImage->GetDimensions()[i] != usImage->GetDimensions()[i]) { MITK_INFO << "Mismatching image dimensions detected in the motion " "compensation filter."; - // TODO: Throw some error here throw std::invalid_argument("Both images must have the same length in each dimension."); } } } void mitk::PhotoacousticMotionCorrectionFilter::InitializeOutput(mitk::Image::Pointer paInput, mitk::Image::Pointer usInput, mitk::Image::Pointer paOutput, mitk::Image::Pointer usOutput) { if (paOutput->GetDimension() != 3) { MITK_INFO << "I jump in here."; this->SetOutputData(paInput, paOutput); this->SetOutputData(usInput, usOutput); } for (unsigned int i = 0; i < usOutput->GetDimension(); i++) { if (usOutput->GetDimensions()[i] != usInput->GetDimensions()[i]) { this->SetOutputData(paInput, paOutput); this->SetOutputData(usInput, usOutput); break; } } } void mitk::PhotoacousticMotionCorrectionFilter::SetOutputData(mitk::Image::Pointer input, mitk::Image::Pointer output) { output->Initialize(input); mitk::ImageReadAccessor accessor(input); output->SetImportVolume(accessor.GetData()); } void mitk::PhotoacousticMotionCorrectionFilter::PerformCorrection(mitk::Image::Pointer paInput, mitk::Image::Pointer usInput, mitk::Image::Pointer paOutput, mitk::Image::Pointer usOutput) { for(unsigned int i = 0; i < paInput->GetDimensions()[2]; i++) { // Get the 2d matrix from slice at i m_PaMatC = this->GetMatrix(paInput, i); m_UsMatC = this->GetMatrix(usInput, i); // Transform to cv::UMat m_PaMat = m_PaMatC.getUMat( cv::ACCESS_READ ).clone(); m_UsMat = m_UsMatC.getUMat( cv::ACCESS_READ ).clone(); // If batch size was set to 0, use one single batch for the whole data set. unsigned int batch; if (m_batch == 0) { batch = paInput->GetDimensions()[2]; } else { batch = m_batch; } // At the beginning of a batch we set the new reference image and directly write the result images if (i % batch == 0) { m_UsRef = m_UsMat.clone(); m_UsRes = m_UsMatC.clone(); m_PaRes = m_PaMatC.clone(); } else { cv::calcOpticalFlowFarneback(m_UsRef, m_UsMat, m_Flow, m_pyr_scale, m_levels, m_winsize, m_iterations, m_poly_n, m_poly_sigma, 0); // Apply the flow to the matrices cv::remap(m_PaMatC, m_PaRes, m_Flow, cv::noArray(), cv::INTER_LINEAR); cv::remap(m_UsMatC, m_UsRes, m_Flow, cv::noArray(), cv::INTER_LINEAR); } // Enter the matrix as a slice at position i into output this->EnterMatrixInPosition(m_PaRes, paOutput, i); this->EnterMatrixInPosition(m_UsRes, usOutput, i); } } cv::Mat mitk::PhotoacousticMotionCorrectionFilter::GetMatrix(const mitk::Image::Pointer input, unsigned int i) { // Access slice i of image input mitk::ImageReadAccessor accessor(input, input->GetSliceData(i)); mitk::Image::Pointer slice = mitk::Image::New(); slice->Initialize(input->GetPixelType(), 2, input->GetDimensions()); slice->SetVolume(accessor.GetData()); // Transform the slice to matrix m_ImageToOpenCVFilter->SetImage(slice); return m_ImageToOpenCVFilter->GetOpenCVMat(); } void mitk::PhotoacousticMotionCorrectionFilter::EnterMatrixInPosition(cv::Mat mat, mitk::Image::Pointer output, unsigned int i) { m_OpenCVToImageFilter->SetOpenCVMat(mat); m_OpenCVToImageFilter->Update(); mitk::Image::Pointer slice = m_OpenCVToImageFilter->GetOutput(); mitk::ImageReadAccessor accessor(slice); output->SetSlice(accessor.GetData(), i); } // TODO: remove debug messages void mitk::PhotoacousticMotionCorrectionFilter::GenerateData() { MITK_INFO << "Start motion compensation."; mitk::Image::Pointer paInput = this->GetInput(0); mitk::Image::Pointer usInput = this->GetInput(1); mitk::Image::Pointer paOutput = this->GetOutput(0); mitk::Image::Pointer usOutput = this->GetOutput(1); // Check that we have two input images with agreeing dimensions try { this->CheckInput(paInput, usInput); } catch(std::invalid_argument& e) { throw e; } // Check the output images and (re-)initialize, if necessary. this->InitializeOutput(paInput, usInput, paOutput, usOutput); // m_ImageToOpenCVFilter->SetImage(paInput); this->PerformCorrection(paInput, usInput, paOutput, usOutput); // MITK_INFO << "Input: " << usInput; // MITK_INFO << "Output: " << usOutput; // MITK_INFO << "Output2: " << paOutput; MITK_INFO << "We succeeded in running through the whole thing!"; }