diff --git a/Modules/PhotoacousticsAlgorithms/include/mitkPhotoacousticMotionCorrectionFilter.h b/Modules/PhotoacousticsAlgorithms/include/mitkPhotoacousticMotionCorrectionFilter.h index f6613ca47f..54454c80e1 100644 --- a/Modules/PhotoacousticsAlgorithms/include/mitkPhotoacousticMotionCorrectionFilter.h +++ b/Modules/PhotoacousticsAlgorithms/include/mitkPhotoacousticMotionCorrectionFilter.h @@ -1,119 +1,288 @@ /*=================================================================== 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 +// 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 +#include + +#define IMAGE_DIMENSION 3 /*!< All images need to have dimension 3*/ + +namespace mitk { +/*! +* \brief Class implementing a mitk::ImageToImageFilter for PAUS motion +* correction. +* +* The filter takes a stack of PA and US images. It then computes the optical +* flow +* within the US image and compensates the PA and US images for the flow. +* Afterwards it +* returns the stack of PA and US images. +* +* @see +* https://docs.opencv.org/3.0-beta/modules/video/doc/motion_analysis_and_object_tracking.html#calcopticalflowfarneback +*/ +class MITKPHOTOACOUSTICSALGORITHMS_EXPORT PhotoacousticMotionCorrectionFilter + : public ImageToImageFilter { +public: + mitkClassMacro(PhotoacousticMotionCorrectionFilter, ImageToImageFilter); + + itkFactorylessNewMacro(Self); + + // Setters and Getters for the class variables + itkSetMacro(BatchSize, unsigned int); + itkSetMacro(PyrScale, double); + itkSetMacro(Levels, unsigned int); + itkSetMacro(WinSize, unsigned int); + itkSetMacro(Iterations, unsigned int); + itkSetMacro(PolyN, unsigned int); + itkSetMacro(PolySigma, double); + itkSetMacro(Flags, unsigned int); + itkGetConstMacro(BatchSize, unsigned int); + itkGetConstMacro(PyrScale, double); + itkGetConstMacro(Levels, unsigned int); + itkGetConstMacro(WinSize, unsigned int); + itkGetConstMacro(Iterations, unsigned int); + itkGetConstMacro(PolyN, unsigned int); + itkGetConstMacro(PolySigma, double); + itkGetConstMacro(Flags, unsigned int); + + // Wrapper for SetInput, GetInput and GetOutput + /*! + * \brief Wrapper which sets the photoacoustic image as the correct input + * + * This method is a wrapper around the @c SetInput method. It is implemented + * for convenience such that you do not have to remember which input is for + * which image. + * + * @param input The photoacoustic image + */ + void SetPaInput(mitk::Image::Pointer); + /*! + * \brief Wrapper which gets the photoacoustic image out of the correct input + * + * This method is a wrapper around the @c GetInput method. It is implemented + * for convenience such that you do not have to remember which input is for + * which image. + * + * @return The photoacoustic image + */ + mitk::Image::Pointer GetPaInput(); + /*! + * \brief Wrapper which sets the ultrasonic image as the correct input + * + * This method is a wrapper around the @c SetInput method. It is implemented + * for convenience such that you do not have to remember which input is for + * which image. + * + * @param input The ultrasonic image + */ + void SetUsInput(mitk::Image::Pointer); + /*! + * \brief Wrapper which gets the ultrasonic image out of the correct input + * + * This method is a wrapper around the @c GetInput method. It is implemented + * for convenience such that you do not have to remember which input is for + * which image. + * + * @return The ultrasonic image + */ + mitk::Image::Pointer GetUsInput(); + /*! + * \brief Wrapper which gets the photoacoustic image out of the correct output + * + * This method is a wrapper around the @c GetOutput method. It is implemented + * for convenience such that you do not have to remember which output is for + * which image. + * + * @return The photoacoustic image + */ + mitk::Image::Pointer GetPaOutput(); + /*! + * \brief Wrapper which gets the ultrasonic image out of the correct output + * + * This method is a wrapper around the @c GetOutput method. It is implemented + * for convenience such that you do not have to remember which output is for + * which image. + * + * @return The ultrasonic image + */ + mitk::Image::Pointer GetUsOutput(); + +protected: + PhotoacousticMotionCorrectionFilter(); + + ~PhotoacousticMotionCorrectionFilter() override; -namespace mitk -{ /*! - * \brief Class implementing a mitk::ImageToImageFilter for PAUS motion correction. - * - * The filter takes a stack of PA and US images. It then computes the optical flow - * within the US image and compensates the PA and US images for the flow. Afterwards it - * returns the stack of PA and US images. - * - * @param m_Batch Determines how many slices belong together and will be motion compensated with regard to the first image in the batch. If the variable is set to 0, the whole time series will be processed as one batch. - * @param m_PyrScale See @c pyr_scale in @c cv::calcOpticalFlowFarneback - * @param m_Levels See @c levels in @c cv::calcOpticalFlowFarneback - * @param m_WinSize See @c winsize in @c cv::calcOpticalFlowFarneback - * @param m_Iterations See @c iterations in @c cv::calcOpticalFlowFarneback - * @param m_PolyN See @c poly_n in @c cv::calcOpticalFlowFarneback - * @param m_PolySigma See @c poly_sigma in @c cv::calcOpticalFlowFarneback - * @param m_Flags See @c flags in @c cv::calcOpticalFlowFarneback - * @see https://docs.opencv.org/3.0-beta/modules/video/doc/motion_analysis_and_object_tracking.html#calcopticalflowfarneback - */ - class MITKPHOTOACOUSTICSALGORITHMS_EXPORT PhotoacousticMotionCorrectionFilter : public ImageToImageFilter - { - public: - mitkClassMacro(PhotoacousticMotionCorrectionFilter, ImageToImageFilter); - - itkFactorylessNewMacro(Self); - - // Setters and Getters for the class variables - itkSetMacro(Batch, unsigned int); - itkSetMacro(PyrScale, double); - itkSetMacro(Levels, unsigned int); - itkSetMacro(WinSize, unsigned int); - itkSetMacro(Iterations, unsigned int); - itkSetMacro(PolyN, unsigned int); - itkSetMacro(PolySigma, double); - itkSetMacro(Flags, unsigned int); - itkGetMacro(Batch, unsigned int); - itkGetMacro(PyrScale, double); - itkGetMacro(Levels, unsigned int); - itkGetMacro(WinSize, unsigned int); - itkGetMacro(Iterations, unsigned int); - itkGetMacro(PolyN, unsigned int); - itkGetMacro(PolySigma, double); - itkGetMacro(Flags, unsigned int); - - // Wrapper for SetInput, GetInput and GetOutput - void SetPaInput(mitk::Image::Pointer); - mitk::Image::Pointer GetPaInput(); - void SetUsInput(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_PyrScale, m_PolySigma; - unsigned int m_Levels, m_WinSize, m_Iterations, m_PolyN, 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 - mitk::OpenCVToMitkImageFilter::Pointer m_OpenCVToImageFilter = mitk::OpenCVToMitkImageFilter::New(); - mitk::ImageToOpenCVImageFilter::Pointer m_ImageToOpenCVFilter = mitk::ImageToOpenCVImageFilter::New(); - }; + * \brief Apply OpenCV algorithm to compensate motion in a 2d image time + * series + * + * This method uses two 3d mitk images. Both will be interpreted as time + * series of 2d images. @c GetInput(0) should be a photoacoustic image whereas + * @c GetInput(1) should be an ultrasound image. The input will be validated + * and then converted to OpenCV matrices. In the end the Farneback algorithm + * will be used to compute the optical flow in consecutive images and + * compensate for this flow. The Output will be two 3d mitk images of the same + * dimensions as the input containing the compensated data. + * + * @warning The input images need to be 3-dimensional (with the same size in + * each dimension). Otherwise, an @c invalid_argument exception will be + * thrown. + * @throws invalid_argument + */ + void GenerateData() override; + /*! + * \brief Validate the input images + * + * The input images have to be non-empty, 3d and have to coincide in the + * length in each dimension. If any of these conditions are violated, the + * method will throw an @c invalid_argument exception. + * + * @param paImage A mitk image + * @param usImage A mitk image + * @warning If the two images are not 3d and do not coincide in the length in + * each dimension, this method will throw an @c invalid_argument exception. + * @throws invalid_argument + */ + void CheckInput(mitk::Image::Pointer paImage, mitk::Image::Pointer usImage); + /*! + * \brief Assure that the output images have the same dimensions as the input + * images. + * + * The output images need to have the same dimensions as the input images. + * This will be checked here. If the dimensions do not match, the output will + * be reinitialized and the image data from the input images will be copied to + * the output images (in order to make sure that they have a valid data + * pointer). + * + * @param paInput Pointer to the photoacoustic input image + * @param usInput Pointer to the ultrasonic input image + * @param paOutput Pointer to the photoacoustic output image + * @param usOutput Pointer to the ultrasonic output image + */ + void InitializeOutputIfNecessary(mitk::Image::Pointer paInput, + mitk::Image::Pointer usInput, + mitk::Image::Pointer paOutput, + mitk::Image::Pointer usOutput); + /*! + * \brief Copy the image data from the input image to the output image + * + * This method copys the image data from @p input to @p output. This method + * assumes that the dimensions of the two images match and will not test this. + * + * @param input A mitk image + * @param output A mitk image + */ + void InitializeOutput(mitk::Image::Pointer input, mitk::Image::Pointer output); + /*! + * \brief This method performs the actual motion compensation. + * + * This method uses the ultrasonic input image @p usInput to compute the + * optical flow in the time series of 2d images. Then it compensates both the + * @p usInput and @p paInput for it and saves the result in @p usOutput and @p + * paOutput respectively. In the background the OpenCV Farneback algorithm is + * used for the flow determination. + * + * @param paInput The photoacoustic input image + * @param usInput The ultrasonic input image + * @param paOutput The photoacoustic output image + * @param usOutput The ultrasonic output image + */ + void PerformCorrection(mitk::Image::Pointer paInput, + mitk::Image::Pointer usInput, + mitk::Image::Pointer paOutput, + mitk::Image::Pointer usOutput); + /*! + * \brief Extract a 2d slice as OpenCV matrix. + * + * This method extracts slice @p i from the 3-dimensional image @p input and + * converts it to a OpenCV matrix. Internally, the + * mitkImageToOpenCVImageFilter is used. + * + * @param input A 3d image from which a slice is extracted as a 2d OpenCV + * matrix. + * @param i Determines the slice to be extracted. + * @return returns a OpenCV matrix containing the 2d slice. + */ + cv::Mat GetMatrix(const mitk::Image::Pointer input, unsigned int i); + /*! + * \brief Insert a OpenCV matrix as a slice into an image + * + * This method converts the 2d OpenCV matrix @p mat into an mitk image using + * the mitkOpenCVToMitkImageFilter. Afterwards it inserts the image as slice + * @p i into the 3d mitk image @p output. + * + * @param mat The matrix to be inserted as a slice + * @param output The 3d image the matrix is inserted into + * @param i The index of the slice to be replaced. + */ + void InsertMatrixAsSlice(cv::Mat mat, mitk::Image::Pointer output, + unsigned int i); + +private: + // Parameters + double m_PyrScale; /*!< See @c pyr_scale in @c cv::calcOpticalFlowFarneback +*/ + double m_PolySigma; /*!< See @c poly_sigma in @c cv::calcOpticalFlowFarneback +*/ + unsigned int + m_Levels; /*!< See @c levels in @c cv::calcOpticalFlowFarneback */ + unsigned int + m_WinSize; /*!< See @c winsize in @c cv::calcOpticalFlowFarneback */ + unsigned int + m_Iterations; /*!< See @c iterations in @c cv::calcOpticalFlowFarneback */ + unsigned int m_PolyN; /*!< See @c poly_n in @c cv::calcOpticalFlowFarneback */ + unsigned int m_Flags; /*!< See @c flags in @c cv::calcOpticalFlowFarneback */ + unsigned int + m_BatchSize; /*!< Determines how many slices belong together and will be + motion compensated with regard to the first image in the + batch. If the variable is set to 0, the whole time series will + be processed as one batch. */ + + // Stuff that OpenCV needs + cv::Mat m_UsRef; /*!< Contains the reference ultrasonic image to which the + motion compensation is compared to.*/ + cv::Mat + m_Flow; /*!< Contains the optical flow between @c m_UsRef and @c m_UsMat*/ + cv::Mat m_PaRes; /*!< Contains the motion compensated photoacoustic image*/ + cv::Mat m_UsRes; /*!< Contains the motion compensated ultrasonic image*/ + cv::Mat m_PaMat; /*!< Contains the latest photoacoustic image to be motion compensated*/ + cv::Mat m_UsMat; /*!< Contains the latest ultrasonic image on which the + optical flow is to be computed */ + + mitk::OpenCVToMitkImageFilter::Pointer m_OpenCVToImageFilter = + mitk::OpenCVToMitkImageFilter::New(); /*!< Filter which converts an OpenCV matrix to a mitk image */ + mitk::ImageToOpenCVImageFilter::Pointer m_ImageToOpenCVFilter = + mitk::ImageToOpenCVImageFilter::New(); /*!< Filter which converts a mitk image to an OpenCV matrix */ +}; } #endif diff --git a/Modules/PhotoacousticsAlgorithms/source/mitkPhotoacousticMotionCorrectionFilter.cpp b/Modules/PhotoacousticsAlgorithms/source/mitkPhotoacousticMotionCorrectionFilter.cpp index f11f914af2..bf17f312d0 100644 --- a/Modules/PhotoacousticsAlgorithms/source/mitkPhotoacousticMotionCorrectionFilter.cpp +++ b/Modules/PhotoacousticsAlgorithms/source/mitkPhotoacousticMotionCorrectionFilter.cpp @@ -1,310 +1,204 @@ /*=================================================================== 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_BatchSize = 5; m_PyrScale = 0.5; m_Levels = 1; m_WinSize = 40; m_Iterations = 2; m_PolyN = 7; m_PolySigma = 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() {} -/*! - * \brief Wrapper which sets the photoacoustic image as the correct input - * - * This method is a wrapper around the @c SetInput method. It is implemented for convenience such that you do not have to remember which input is for which image. - * - * @param input The photoacoustic image - */ -void mitk::PhotoacousticMotionCorrectionFilter::SetPaInput(mitk::Image::Pointer input) { +void mitk::PhotoacousticMotionCorrectionFilter::SetPaInput( + mitk::Image::Pointer input) { this->SetInput(0, input); } -/*! - * \brief Wrapper which gets the photoacoustic image out of the correct input - * - * This method is a wrapper around the @c GetInput method. It is implemented for convenience such that you do not have to remember which input is for which image. - * - * @return The photoacoustic image - */ mitk::Image::Pointer mitk::PhotoacousticMotionCorrectionFilter::GetPaInput() { return this->GetInput(0); } -/*! - * \brief Wrapper which sets the ultrasonic image as the correct input - * - * This method is a wrapper around the @c SetInput method. It is implemented for convenience such that you do not have to remember which input is for which image. - * - * @param input The ultrasonic image - */ -void mitk::PhotoacousticMotionCorrectionFilter::SetUsInput(mitk::Image::Pointer input) { +void mitk::PhotoacousticMotionCorrectionFilter::SetUsInput( + mitk::Image::Pointer input) { this->SetInput(1, input); } -/*! - * \brief Wrapper which gets the ultrasonic image out of the correct input - * - * This method is a wrapper around the @c GetInput method. It is implemented for convenience such that you do not have to remember which input is for which image. - * - * @return The ultrasonic image - */ mitk::Image::Pointer mitk::PhotoacousticMotionCorrectionFilter::GetUsInput() { return this->GetInput(1); } -/*! - * \brief Wrapper which gets the photoacoustic image out of the correct output - * - * This method is a wrapper around the @c GetOutput method. It is implemented for convenience such that you do not have to remember which output is for which image. - * - * @return The photoacoustic image - */ mitk::Image::Pointer mitk::PhotoacousticMotionCorrectionFilter::GetPaOutput() { return this->GetOutput(0); } -/*! - * \brief Wrapper which gets the ultrasonic image out of the correct output - * - * This method is a wrapper around the @c GetOutput method. It is implemented for convenience such that you do not have to remember which output is for which image. - * - * @return The ultrasonic image - */ mitk::Image::Pointer mitk::PhotoacousticMotionCorrectionFilter::GetUsOutput() { return this->GetOutput(1); } -/*! - * \brief Validate the input images - * - * The input images have to be non-empty, 3d and have to coincide in the length in each dimension. If any of these conditions are violated, the method will throw an @c invalid_argument exception. - * - * @param paImage A mitk image - * @param usImage A mitk image - * @warning If the two images are not 3d and do not coincide in the length in each dimension, this method will throw an @c invalid_argument exception. - */ -void mitk::PhotoacousticMotionCorrectionFilter::CheckInput(mitk::Image::Pointer paImage, mitk::Image::Pointer usImage) { +void mitk::PhotoacousticMotionCorrectionFilter::CheckInput( + mitk::Image::Pointer paImage, mitk::Image::Pointer usImage) { // Check that we actually got some images if (!paImage || !usImage) { - MITK_INFO << "We did not get two images!"; + MITK_ERROR << "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."; + if (paImage->GetDimension() != IMAGE_DIMENSION || usImage->GetDimension() != IMAGE_DIMENSION) { + MITK_ERROR << "Mismatching image dimensions detected in the motion " + "compensation filter."; 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++) { + 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."; - throw std::invalid_argument("Both images must have the same length in each dimension."); + MITK_ERROR << "Mismatching image dimensions detected in the motion " + "compensation filter."; + throw std::invalid_argument( + "Both images must have the same length in each dimension."); } } - } -/*! - * \brief Assure that the output images have the same dimensions as the input images. - * - * The output images need to have the same dimensions as the input images. This will be checked here. If the dimensions do not match, the output will be reinitialized and the image data from the input images will be copied to the output images (in order to make sure that they have a valid data pointer). - * - * @param paInput Pointer to the photoacoustic input image - * @param usInput Pointer to the ultrasonic input image - * @param paOutput Pointer to the photoacoustic output image - * @param usOutput Pointer to the ultrasonic output image - */ -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); +void mitk::PhotoacousticMotionCorrectionFilter::InitializeOutputIfNecessary( + mitk::Image::Pointer paInput, mitk::Image::Pointer usInput, + mitk::Image::Pointer paOutput, mitk::Image::Pointer usOutput) { + if (paOutput->GetDimension() != IMAGE_DIMENSION) { + this->InitializeOutput(paInput, paOutput); + this->InitializeOutput(usInput, usOutput); } - for (unsigned int i = 0; i < usOutput->GetDimension(); i++) { + for (unsigned int i = 0; i < usOutput->GetDimension(); ++i) { if (usOutput->GetDimensions()[i] != usInput->GetDimensions()[i]) { - this->SetOutputData(paInput, paOutput); - this->SetOutputData(usInput, usOutput); + this->InitializeOutput(paInput, paOutput); + this->InitializeOutput(usInput, usOutput); break; - } + } } } -/*! - * \brief Copy the image data from the input image to the output image - * - * This method copys the image data from @p input to @p output. This method assumes that the dimensions of the two images match and will not test this. - * - * @param input A mitk image - * @param output A mitk image - */ -void mitk::PhotoacousticMotionCorrectionFilter::SetOutputData(mitk::Image::Pointer input, mitk::Image::Pointer output) { +void mitk::PhotoacousticMotionCorrectionFilter::InitializeOutput( + mitk::Image::Pointer input, mitk::Image::Pointer output) { output->Initialize(input); mitk::ImageReadAccessor accessor(input); output->SetImportVolume(accessor.GetData()); } -/*! - * \brief This method performs the actual motion compensation. - * - * This method uses the ultrasonic input image @p usInput to compute the optical flow in the time series of 2d images. Then it compensates both the @p usInput and @p paInput for it and saves the result in @p usOutput and @p paOutput respectively. In the background the OpenCV Farneback algorithm is used for the flow determination. - * - * @param paInput The photoacoustic input image - * @param usInput The ultrasonic input image - * @param paOutput The photoacoustic output image - * @param usOutput The ultrasonic output image - */ -void mitk::PhotoacousticMotionCorrectionFilter::PerformCorrection(mitk::Image::Pointer paInput, mitk::Image::Pointer usInput, mitk::Image::Pointer paOutput, mitk::Image::Pointer usOutput) { +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); + // If batch size was set to 0, use one single batch for the whole data set. + unsigned int batch; + if (m_BatchSize == 0) { + batch = paInput->GetDimensions()[IMAGE_DIMENSION - 1]; + } else { + batch = m_BatchSize; + } - // Transform to cv::UMat - m_PaMat = m_PaMatC.getUMat( cv::ACCESS_READ ).clone(); - m_UsMat = m_UsMatC.getUMat( cv::ACCESS_READ ).clone(); + for (unsigned int i = 0; i < paInput->GetDimensions()[IMAGE_DIMENSION - 1]; ++i) { - // 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; - } + // Get the 2d matrix from slice at i + m_PaMat = this->GetMatrix(paInput, i); + m_UsMat = this->GetMatrix(usInput, i); - // At the beginning of a batch we set the new reference image and directly write the result images - + // 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(); + m_UsRes = m_UsMat.clone(); + m_PaRes = m_PaMat.clone(); } else { - cv::calcOpticalFlowFarneback(m_UsRef, m_UsMat, m_Flow, m_PyrScale, m_Levels, - m_WinSize, m_Iterations, m_PolyN, + cv::calcOpticalFlowFarneback(m_UsRef, m_UsMat, m_Flow, m_PyrScale, + m_Levels, m_WinSize, m_Iterations, m_PolyN, m_PolySigma, m_Flags); + // TODO: reshape flow // 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); + cv::remap(m_PaMat, m_PaRes, m_Flow, cv::noArray(), cv::INTER_LINEAR); + cv::remap(m_UsMat, 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); + this->InsertMatrixAsSlice(m_PaRes, paOutput, i); + this->InsertMatrixAsSlice(m_UsRes, usOutput, i); } } -/*! - * \brief Extract a 2d slice as OpenCV matrix. - * - * This method extracts slice @p i from the 3-dimensional image @p input and converts it to a OpenCV matrix. Internally, the mitkImageToOpenCVImageFilter is used. - * - * @param input A 3d image from which a slice is extracted as a 2d OpenCV matrix. - * @param i Determines the slice to be extracted. - * @return returns a OpenCV matrix containing the 2d slice. - */ - -cv::Mat mitk::PhotoacousticMotionCorrectionFilter::GetMatrix(const mitk::Image::Pointer input, unsigned int 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->Initialize(input->GetPixelType(), IMAGE_DIMENSION - 1, input->GetDimensions()); slice->SetVolume(accessor.GetData()); // Transform the slice to matrix m_ImageToOpenCVFilter->SetImage(slice); return m_ImageToOpenCVFilter->GetOpenCVMat(); } -/*! - * \brief Insert a OpenCV matrix as a slice into an image - * - * This method converts the 2d OpenCV matrix @p mat into an mitk image using the mitkOpenCVToMitkImageFilter. Afterwards it inserts the image as slice @p i into the 3d mitk image @p output. - * - * @param mat The matrix to be inserted as a slice - * @param output The 3d image the matrix is inserted into - * @param i The index of the slice to be replaced. - */ - -void mitk::PhotoacousticMotionCorrectionFilter::EnterMatrixInPosition(cv::Mat mat, mitk::Image::Pointer output, unsigned int i) { +void mitk::PhotoacousticMotionCorrectionFilter::InsertMatrixAsSlice( + 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); } - -/*! - * \brief Apply OpenCV algorithm to compensate motion in a 2d image time series - * - * This method uses two 3d mitk images. Both will be interpreted as time series of 2d images. @c GetInput(0) should be a photoacoustic image whereas @c GetInput(1) should be an ultrasound image. The input will be validated and then converted to OpenCV matrices. In the end the Farneback algorithm will be used to compute the optical flow in consecutive images and compensate for this flow. The Output will be two 3d mitk images of the same dimensions as the input containing the compensated data. - * - * @warning The input images need to be 3-dimensional (with the same size in each dimension). Otherwise, an @c invalid_argument exception will be thrown. - */ // 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; - } + this->CheckInput(paInput, usInput); // Check the output images and (re-)initialize, if necessary. - this->InitializeOutput(paInput, usInput, paOutput, usOutput); + this->InitializeOutputIfNecessary(paInput, usInput, paOutput, usOutput); // m_ImageToOpenCVFilter->SetImage(paInput); this->PerformCorrection(paInput, usInput, paOutput, usOutput); MITK_INFO << "Motion compensation accomplished."; } diff --git a/Modules/PhotoacousticsAlgorithms/test/mitkPhotoacousticMotionCorrectionFilterTest.cpp b/Modules/PhotoacousticsAlgorithms/test/mitkPhotoacousticMotionCorrectionFilterTest.cpp index d44c3ebdd3..8ac6ea91c6 100644 --- a/Modules/PhotoacousticsAlgorithms/test/mitkPhotoacousticMotionCorrectionFilterTest.cpp +++ b/Modules/PhotoacousticsAlgorithms/test/mitkPhotoacousticMotionCorrectionFilterTest.cpp @@ -1,182 +1,182 @@ /*=================================================================== 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 class mitkPhotoacousticMotionCorrectionFilterTestSuite : public mitk::TestFixture { CPPUNIT_TEST_SUITE(mitkPhotoacousticMotionCorrectionFilterTestSuite); MITK_TEST(testSettingFirstInput); MITK_TEST(testSettingSecondInput); // MITK_TEST(testNoThirdInput); MITK_TEST(testGetFirstEmptyOutput); MITK_TEST(testGetSecondEmptyOutput); MITK_TEST(testSameInOutDimensions); MITK_TEST(testNo3DError); MITK_TEST(testSameInputDimension1); MITK_TEST(testSameInputDimension2); MITK_TEST(testNullPtr1); MITK_TEST(testNullPtr2); MITK_TEST(testNullPtr3); MITK_TEST(testSameInputDimensions); MITK_TEST(testStaticSliceCorrection); CPPUNIT_TEST_SUITE_END(); private: mitk::PhotoacousticMotionCorrectionFilter::Pointer filter; mitk::Image::Pointer image, image4d, image2d, image3d1slice; float * data2d; public: void setUp() override { // get the filter I need filter = mitk::PhotoacousticMotionCorrectionFilter::New(); // get a 3d mitk image image = mitk::Image::New(); image2d = mitk::Image::New(); image4d = mitk::Image::New(); image3d1slice = mitk::Image::New(); unsigned int * dimensions = new unsigned int[3] {2, 2, 2}; unsigned int * dimensions2 = new unsigned int[3] {2, 2, 1}; unsigned int * dim2 = new unsigned int[4] {2, 2, 2, 2}; mitk::PixelType pt = mitk::MakeScalarPixelType(); float * data = new float[8] {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; float * data4d = new float[16]; data2d = new float[4] {1.0, 0.0, 0.0, 0.0}; image->Initialize(pt, 3, dimensions); image->SetVolume(data); image4d->Initialize(pt, 4, dim2); image4d->SetVolume(data4d); image2d->Initialize(pt, 2, dimensions); image2d->SetVolume(data2d); image3d1slice->Initialize(pt, 3, dimensions2); image3d1slice->SetVolume(data); delete[] dimensions; delete[] dimensions2; delete[] data; delete[] data4d; delete[] dim2; } void tearDown() override { delete[] data2d; } void testSettingFirstInput() { filter->SetInput(0, image); mitk::Image::Pointer out = filter->GetInput(0); CPPUNIT_ASSERT_EQUAL(image, out); } void testSettingSecondInput() { filter->SetInput(1, image); mitk::Image::Pointer out = filter->GetInput(1); CPPUNIT_ASSERT_EQUAL(image, out); } // void testNoThirdInput() { // CPPUNIT_ASSERT_NO_THROW(filter->SetInput(2, image)); // // CPPUNIT_ASSERT(true); // } void testGetFirstEmptyOutput() { mitk::Image::Pointer out = filter->GetOutput(0); unsigned int dim = 0; CPPUNIT_ASSERT_EQUAL(dim, out->GetDimension()); } void testGetSecondEmptyOutput() { mitk::Image::Pointer out = filter->GetOutput(1); unsigned int dim = 0; CPPUNIT_ASSERT_EQUAL(dim, out->GetDimension()); } void testSameInOutDimensions() { filter->SetInput(0, image); filter->SetInput(1, image); filter->Update(); mitk::Image::Pointer out = filter->GetOutput(0); CPPUNIT_ASSERT_EQUAL(image->GetDimension(), out->GetDimension()); out = filter->GetOutput(1); CPPUNIT_ASSERT_EQUAL(image->GetDimension(), out->GetDimension()); } void testNo3DError() { filter->SetInput(0, image4d); filter->SetInput(1, image4d); CPPUNIT_ASSERT_THROW(filter->Update(), std::invalid_argument); } // I only test for dim input0 <= dim input1, because otherwise // itk will throw an error before the program even reaches the // parts I wrote. // Same for testSameInputDimension2 void testSameInputDimension1() { filter->SetInput(0, image); filter->SetInput(1, image4d); CPPUNIT_ASSERT_THROW(filter->Update(), std::invalid_argument); } // See previous comment void testSameInputDimension2() { filter->SetInput(0, image2d); filter->SetInput(1, image); CPPUNIT_ASSERT_THROW(filter->Update(), std::invalid_argument); } // I tried to catch the error myself in the filter, but itk does some magic beforehand void testNullPtr1() { - filter->SetInput(0, NULL); + filter->SetInput(0, nullptr); filter->SetInput(1, image); CPPUNIT_ASSERT_THROW(filter->Update(), itk::ExceptionObject); } // Now I am allowed to catch it myself, because the first input is fine -.- void testNullPtr2() { filter->SetInput(0, image); - filter->SetInput(1, NULL); + filter->SetInput(1, nullptr); CPPUNIT_ASSERT_THROW(filter->Update(), std::invalid_argument); } void testNullPtr3() { - filter->SetInput(0, NULL); - filter->SetInput(1, NULL); + filter->SetInput(0, nullptr); + filter->SetInput(1, nullptr); CPPUNIT_ASSERT_THROW(filter->Update(), itk::ExceptionObject); } void testSameInputDimensions() { filter->SetInput(0, image3d1slice); filter->SetInput(1, image); CPPUNIT_ASSERT_THROW(filter->Update(), std::invalid_argument); } void testStaticSliceCorrection() { // image->SetSlice(data2d, 0); // image->SetSlice(data2d, 1); filter->SetInput(0, image); filter->SetInput(1, image); filter->Update(); mitk::Image::Pointer out1 = filter->GetOutput(0); mitk::Image::Pointer out2 = filter->GetOutput(1); MITK_ASSERT_EQUAL(image, out1, "Check that static image does not get changed."); } }; MITK_TEST_SUITE_REGISTRATION(mitkPhotoacousticMotionCorrectionFilter)