diff --git a/Modules/PhotoacousticsAlgorithms/include/mitkPhotoacousticMotionCorrectionFilter.h b/Modules/PhotoacousticsAlgorithms/include/mitkPhotoacousticMotionCorrectionFilter.h index 54454c80e1..210d427a0d 100644 --- a/Modules/PhotoacousticsAlgorithms/include/mitkPhotoacousticMotionCorrectionFilter.h +++ b/Modules/PhotoacousticsAlgorithms/include/mitkPhotoacousticMotionCorrectionFilter.h @@ -1,288 +1,294 @@ /*=================================================================== 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 #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); +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); + 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); + // 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(); + // 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(); + protected: + PhotoacousticMotionCorrectionFilter(); - ~PhotoacousticMotionCorrectionFilter() override; + ~PhotoacousticMotionCorrectionFilter() override; - /*! - * \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); + /*! + * \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); + + /*! + * \brief Compute the remapping map from an optical flow + * + * The optical flow cannot be used directly to compensate an image. Instead we have to generate an appropriate map. + * + * @param flow The optical flow which is the base for the remapping. + * @return The remapping map. + */ + cv::Mat ComputeFlowMap(cv::Mat); -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. */ + 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 */ + // 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 */ + cv::Mat m_Map; /*!< Contains the remapping map */ - mitk::OpenCVToMitkImageFilter::Pointer m_OpenCVToImageFilter = + 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::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 bf17f312d0..f4733ed6b7 100644 --- a/Modules/PhotoacousticsAlgorithms/source/mitkPhotoacousticMotionCorrectionFilter.cpp +++ b/Modules/PhotoacousticsAlgorithms/source/mitkPhotoacousticMotionCorrectionFilter.cpp @@ -1,204 +1,219 @@ /*=================================================================== 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_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() {} void mitk::PhotoacousticMotionCorrectionFilter::SetPaInput( mitk::Image::Pointer input) { this->SetInput(0, input); } mitk::Image::Pointer mitk::PhotoacousticMotionCorrectionFilter::GetPaInput() { return this->GetInput(0); } void mitk::PhotoacousticMotionCorrectionFilter::SetUsInput( 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) { 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() != 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) { if (paImage->GetDimensions()[i] != usImage->GetDimensions()[i]) { 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."); } } } 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) { if (usOutput->GetDimensions()[i] != usInput->GetDimensions()[i]) { this->InitializeOutput(paInput, paOutput); this->InitializeOutput(usInput, usOutput); break; } } } void mitk::PhotoacousticMotionCorrectionFilter::InitializeOutput( 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) { // 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; } for (unsigned int i = 0; i < paInput->GetDimensions()[IMAGE_DIMENSION - 1]; ++i) { // 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 if (i % batch == 0) { m_UsRef = m_UsMat.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, m_PolySigma, m_Flags); - // TODO: reshape flow + m_Map = this->ComputeFlowMap(m_Flow); + // Apply the flow to the matrices - 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); + cv::remap(m_PaMat, m_PaRes, m_Map, cv::noArray(), cv::INTER_LINEAR); + cv::remap(m_UsMat, m_UsRes, m_Map, cv::noArray(), cv::INTER_LINEAR); } // Enter the matrix as a slice at position i into output this->InsertMatrixAsSlice(m_PaRes, paOutput, i); this->InsertMatrixAsSlice(m_UsRes, usOutput, i); } } +// Copied from https://stackoverflow.com/questions/17459584/opencv-warping-image-based-on-calcopticalflowfarneback +cv::Mat mitk::PhotoacousticMotionCorrectionFilter::ComputeFlowMap(cv::Mat flow) { + cv::Mat map(flow.size(), flow.type()); + + for (int y = 0; y < map.rows; ++y) { + for(int x = 0; x < map.cols; ++x) { + cv::Point2f f = flow.at(y,x); + map.at(y,x) = cv::Point2f(x + f.x, y + f.y); + } + } + + return map; +} + 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(), IMAGE_DIMENSION - 1, input->GetDimensions()); slice->SetVolume(accessor.GetData()); // Transform the slice to matrix m_ImageToOpenCVFilter->SetImage(slice); return m_ImageToOpenCVFilter->GetOpenCVMat(); } 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); } // 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 this->CheckInput(paInput, usInput); // Check the output images and (re-)initialize, if necessary. 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 8ac6ea91c6..7f87cf7316 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, 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, nullptr); CPPUNIT_ASSERT_THROW(filter->Update(), std::invalid_argument); } void testNullPtr3() { 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); + 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)