diff --git a/Modules/PhotoacousticsAlgorithms/CMakeLists.txt b/Modules/PhotoacousticsAlgorithms/CMakeLists.txt index b54cc87313..b107f07b58 100644 --- a/Modules/PhotoacousticsAlgorithms/CMakeLists.txt +++ b/Modules/PhotoacousticsAlgorithms/CMakeLists.txt @@ -1,16 +1,16 @@ -set(dependencies_list MitkCore MitkAlgorithmsExt) +set(dependencies_list MitkCore MitkAlgorithmsExt MitkOpenCVVideoSupport) IF(MITK_USE_OpenCL) add_definitions(-DPHOTOACOUSTICS_USE_GPU) set(dependencies_list ${dependencies_list} MitkOpenCL) message("Using OpenCL in PhotoacousticAlgorithms") ENDIF(MITK_USE_OpenCL) MITK_CREATE_MODULE( SUBPROJECTS DEPENDS ${dependencies_list} #AUTOLOAD_WITH MitkCore INCLUDE_DIRS PUBLIC Algorithms/ITKUltrasound Algorithms Algorithms/OCL INTERNAL_INCLUDE_DIRS ${INCLUDE_DIRS_INTERNAL} - PACKAGE_DEPENDS ITK|ITKFFT+ITKImageCompose+ITKImageIntensity OpenCV + PACKAGE_DEPENDS ITK|ITKFFT+ITKImageCompose+ITKImageIntensity OpenCV ) diff --git a/Modules/PhotoacousticsAlgorithms/include/mitkPhotoacousticMotionCorrectionFilter.h b/Modules/PhotoacousticsAlgorithms/include/mitkPhotoacousticMotionCorrectionFilter.h index 056b5c19a4..0caebebb08 100644 --- a/Modules/PhotoacousticsAlgorithms/include/mitkPhotoacousticMotionCorrectionFilter.h +++ b/Modules/PhotoacousticsAlgorithms/include/mitkPhotoacousticMotionCorrectionFilter.h @@ -1,80 +1,83 @@ /*=================================================================== 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 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 protected: PhotoacousticMotionCorrectionFilter(); ~PhotoacousticMotionCorrectionFilter() override; void GenerateData() override; //##Description //## @brief Time when Header was last initialized /* itk::TimeStamp m_TimeOfHeaderInitialization; */ private: // InputData - mitk::Image::Pointer m_paImage, m_usImage; + mitk::Image::Pointer m_paImage, m_usImage, m_paCompensated, m_usCompensated; // Parameters double m_pyr_scale, m_poly_sigma; - int m_levels, m_winsize, m_iterations, m_poly_n, m_flags, m_batch; + unsigned int m_levels, m_winsize, m_iterations, m_poly_n, m_flags, m_batch; /* // Stuff that OpenCV needs */ - cv::Mat m_startImage, m_stopImage, m_flow; - /* cv::UMat m_startImageG, m_StopImageG, m_uflow; */ + 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(); }; } #endif diff --git a/Modules/PhotoacousticsAlgorithms/source/mitkPhotoacousticMotionCorrectionFilter.cpp b/Modules/PhotoacousticsAlgorithms/source/mitkPhotoacousticMotionCorrectionFilter.cpp index 960c9d3a06..0f6564521c 100644 --- a/Modules/PhotoacousticsAlgorithms/source/mitkPhotoacousticMotionCorrectionFilter.cpp +++ b/Modules/PhotoacousticsAlgorithms/source/mitkPhotoacousticMotionCorrectionFilter.cpp @@ -1,62 +1,164 @@ /*=================================================================== 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() -{ +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 + // 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; SetNumberOfIndexedInputs(2); - SetNumberOfIndexedOutputs(1); + SetNumberOfIndexedOutputs(2); } -mitk::PhotoacousticMotionCorrectionFilter::~PhotoacousticMotionCorrectionFilter() -{ -} +mitk::PhotoacousticMotionCorrectionFilter:: + ~PhotoacousticMotionCorrectionFilter() {} -// void mitk::PhotoacousticMotionCorrectionFilter::SetInput(Image::Pointer paImage, -// Image::Pointer usImage) +// void mitk::PhotoacousticMotionCorrectionFilter::SetInput(Image::Pointer +// paImage, +// Image::Pointer +// usImage) // { // m_paImage = paImage; // m_usImage = usImage; // } -void mitk::PhotoacousticMotionCorrectionFilter::GenerateData() -{ - auto input1 = this->GetInput(0); - auto input2 = this->GetInput(1); +void mitk::PhotoacousticMotionCorrectionFilter::GenerateData() { + MITK_INFO << "Start motion compensation."; + + m_paImage = this->GetInput(0); + m_usImage = this->GetInput(1); - if(input1 && input2) - { - m_paImage = input1; - m_usImage = input2; + // Check that we actually got some images + if (!m_paImage || !m_usImage) { + // TODO: Throw some error here + } + // Check that the image dimensions are the same + if (m_paImage->GetDimension() != m_usImage->GetDimension() && + m_usImage->GetDimension() == 3) { + MITK_INFO << "Mismatching image dimensions detected in the motion " + "compensation filter."; + // TODO: Throw some error here + } + for (unsigned int i = 0; i < m_paImage->GetDimension(); i++) { + if (m_paImage->GetDimensions()[i] != m_usImage->GetDimensions()[i]) { + MITK_INFO << "Mismatching image dimensions detected in the motion " + "compensation filter."; + // TODO: Throw some error here + } + } - mitk::CastToItkImage(m_paImage, m_itkPaImage); - // TODO: Make sure that only 2d images enter here - m_startImage = itk::OpenCVImageBridge::ITKImageToCVMat< itk::Image >(m_itkPaImage); + // Initialize output images + if (!m_paCompensated) { + m_paCompensated = mitk::Image::New(); + } + if (!m_usCompensated) { + m_usCompensated = mitk::Image::New(); } + m_paCompensated->Initialize(m_paImage->GetPixelType(), + m_paImage->GetDimension(), + m_paImage->GetDimensions()); + m_usCompensated->Initialize(m_usImage->GetPixelType(), + m_usImage->GetDimension(), + m_usImage->GetDimensions()); + + // TODO: remove debug messages + + // Initialize the slices + mitk::Image::Pointer pa_slice = mitk::Image::New(); + mitk::Image::Pointer us_slice = mitk::Image::New(); + pa_slice->Initialize(m_paImage->GetPixelType(), 2, + m_paImage->GetDimensions()); + us_slice->Initialize(m_usImage->GetPixelType(), 2, + m_usImage->GetDimensions()); + + MITK_INFO << "Start iteration."; + // Iterate over all the slices + for (unsigned int i = 0; i < m_paImage->GetDimensions()[2]; i++) { + + // Get a read accessor for each slice + mitk::ImageReadAccessor pa_accessor(m_paImage, m_paImage->GetSliceData(i)); + mitk::ImageReadAccessor us_accessor(m_paImage, m_usImage->GetSliceData(i)); + + // Write the correct image data into the slice + pa_slice->SetImportVolume(pa_accessor.GetData()); + us_slice->SetImportVolume(us_accessor.GetData()); + + // Convert them first to an itk::Image and then to a cv::Mat + mitk::CastToItkImage(pa_slice, m_itkPaImage); + mitk::CastToItkImage(us_slice, m_itkUsImage); + MITK_INFO << "Generate Matrix."; + + m_PaMatC = itk::OpenCVImageBridge::ITKImageToCVMat>( + m_itkPaImage); + m_UsMatC = itk::OpenCVImageBridge::ITKImageToCVMat>( + m_itkUsImage); + + m_PaMat = m_PaMatC.getUMat( cv::ACCESS_READ ); + m_UsMat = m_UsMatC.getUMat( cv::ACCESS_READ ); + + MITK_INFO << "Matrix generated."; + // At the beginning of a batch we set new references and the compensation + // can be skipped. + // TODO: handle m_batch == 0 + // if (i % m_batch == 0) { + // m_UsRef = m_UsMat; + // m_UsRes = m_UsMatC; + // m_PaRes = m_PaMatC; + // continue; + // } else { + // // Calculate the flow using the Farneback algorithm + // // TODO: flags hard coded to 0, rethink. + // 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 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); + // } + + // TODO: Actually do something, not just retransform + m_PaRes = m_PaMatC; + m_UsRes = m_UsMatC; + + m_OpenCVToImageFilter->SetOpenCVMat(m_PaRes); + m_OpenCVToImageFilter->Update(); + pa_slice = m_OpenCVToImageFilter->GetOutput(); + mitk::ImageReadAccessor pa_slice_accessor(pa_slice); + m_paCompensated->SetSlice(pa_slice_accessor.GetData(), i); + m_OpenCVToImageFilter->SetOpenCVMat(m_UsRes); + m_OpenCVToImageFilter->Update(); + us_slice = m_OpenCVToImageFilter->GetOutput(); + mitk::ImageReadAccessor us_slice_accessor(us_slice); + m_paCompensated->SetSlice(us_slice_accessor.GetData(), i); +} + + this->SetNthOutput(1, m_usCompensated); + this->SetNthOutput(0, m_paCompensated); + MITK_INFO << "We succeeded in running through the whole thing!"; }