Medical Imaging Interaction Toolkit  2018.4.99-6a3ea89d
Medical Imaging Interaction Toolkit
mitkPhotoacousticMotionCorrectionFilter.cpp
Go to the documentation of this file.
1 /*============================================================================
2 
3 The Medical Imaging Interaction Toolkit (MITK)
4 
5 Copyright (c) German Cancer Research Center (DKFZ)
6 All rights reserved.
7 
8 Use of this source code is governed by a 3-clause BSD license that can be
9 found in the LICENSE file.
10 
11 ============================================================================*/
12 
14 #include <mitkImageReadAccessor.h>
15 
18  // Set the defaults for the OpticalFlowFarneback algorithm
19  // The values are taken directly out of Thomas's
20  // US-CV-Based-Optical-Flow-Carotis.ipyn
21  m_BatchSize = 5;
22  m_PyrScale = 0.5;
23  m_Levels = 1;
24  m_WinSize = 40;
25  m_Iterations = 2;
26  m_PolyN = 7;
27  m_PolySigma = 1.5;
28  m_Flags = 0;
29 
30  m_MaxValue = 255.0;
31  m_MinValue = 0.0;
32 
33  this->SetNumberOfIndexedInputs(2);
34  this->SetNumberOfIndexedOutputs(2);
35  this->SetNthOutput(0, mitk::Image::New());
36  this->SetNthOutput(1, mitk::Image::New());
37 }
38 
41 
43  mitk::Image::Pointer input) {
44  this->SetInput(0, input);
45 }
46 
48  return this->GetInput(0);
49 }
50 
52  mitk::Image::Pointer input) {
53  this->SetInput(1, input);
54 }
55 
57  return this->GetInput(1);
58 }
59 
61  return this->GetOutput(0);
62 }
63 
65  return this->GetOutput(1);
66 }
67 
69  mitk::Image::Pointer paImage, mitk::Image::Pointer usImage) {
70  // Check that we actually got some images
71  if (!paImage || !usImage) {
72  MITK_ERROR << "We did not get two images!";
73  throw std::invalid_argument("One of the images was NULL.");
74  }
75 
76  // Check that the image dimensions are the same
77  if (paImage->GetDimension() != IMAGE_DIMENSION || usImage->GetDimension() != IMAGE_DIMENSION) {
78  MITK_ERROR << "Mismatching image dimensions detected in the motion "
79  "compensation filter.";
80  throw std::invalid_argument("Both images must have dimension 3.");
81  }
82 
83  // Check that each dimension has the same size
84  for (unsigned int i = 0; i < paImage->GetDimension(); ++i) {
85  if (paImage->GetDimensions()[i] != usImage->GetDimensions()[i]) {
86  MITK_ERROR << "Mismatching image dimensions detected in the motion "
87  "compensation filter.";
88  throw std::invalid_argument(
89  "Both images must have the same length in each dimension.");
90  }
91  }
92 }
93 
96  mitk::Image::Pointer paOutput, mitk::Image::Pointer usOutput) {
97  if (paOutput->GetDimension() != IMAGE_DIMENSION) {
98  this->InitializeOutput(paInput, paOutput);
99  this->InitializeOutput(usInput, usOutput);
100  }
101 
102  for (unsigned int i = 0; i < usOutput->GetDimension(); ++i) {
103  if (usOutput->GetDimensions()[i] != usInput->GetDimensions()[i]) {
104  this->InitializeOutput(paInput, paOutput);
105  this->InitializeOutput(usInput, usOutput);
106  break;
107  }
108  }
109 }
110 
113  output->Initialize(input);
114  mitk::ImageReadAccessor accessor(input);
115  output->SetImportVolume(accessor.GetData());
116 }
117 
120  mitk::Image::Pointer paOutput, mitk::Image::Pointer usOutput) {
121 
122  // If batch size was set to 0, use one single batch for the whole data set.
123  unsigned int batch;
124  if (m_BatchSize == 0) {
125  batch = paInput->GetDimensions()[IMAGE_DIMENSION - 1];
126  } else {
127  batch = m_BatchSize;
128  }
129 
130  for (unsigned int i = 0; i < paInput->GetDimensions()[IMAGE_DIMENSION - 1]; ++i) {
131 
132  // Get the 2d matrix from slice at i
133  m_PaMat = this->GetMatrix(paInput, i);
134  m_UsMat = this->GetMatrix(usInput, i);
135 
136  // At the beginning of a batch we set the new reference image and directly
137  // write the result images
138  if (i % batch == 0) {
139  // Rescale reference to full char resolution
140  m_UsRef = this->FitMatrixToChar(m_UsMat);
141  m_UsRes = m_UsMat.clone();
142  m_PaRes = m_PaMat.clone();
143  } else {
144  cv::Mat UsMatRescaled = this->FitMatrixToChar(m_UsMat);
145  cv::calcOpticalFlowFarneback(m_UsRef, UsMatRescaled, m_Flow, m_PyrScale,
146  m_Levels, m_WinSize, m_Iterations, m_PolyN,
147  m_PolySigma, m_Flags);
148 
149  m_Map = this->ComputeFlowMap(m_Flow);
150 
151  // Apply the flow to the matrices
152  cv::remap(m_PaMat, m_PaRes, m_Map, cv::noArray(), cv::INTER_LINEAR);
153  cv::remap(m_UsMat, m_UsRes, m_Map, cv::noArray(), cv::INTER_LINEAR);
154  }
155 
156  // Enter the matrix as a slice at position i into output
157  this->InsertMatrixAsSlice(m_PaRes, paOutput, i);
158  this->InsertMatrixAsSlice(m_UsRes, usOutput, i);
159  }
160 }
161 
162 // Copied from https://stackoverflow.com/questions/17459584/opencv-warping-image-based-on-calcopticalflowfarneback
164  cv::Mat map(flow.size(), flow.type());
165 
166  for (int y = 0; y < map.rows; ++y) {
167  for(int x = 0; x < map.cols; ++x) {
168  cv::Point2f f = flow.at<cv::Point2f>(y,x);
169  map.at<cv::Point2f>(y,x) = cv::Point2f(x + f.x, y + f.y);
170  }
171  }
172 
173  return map;
174 }
175 
177 
178  if (m_MaxValue == m_MinValue) {
179 
180  return mat.clone();
181  }
182 
183  return MAX_MATRIX*(mat.clone() - m_MinValue) / (m_MaxValue - m_MinValue);
184 }
185 
187  const mitk::Image::Pointer input, unsigned int i) {
188 
189  // Access slice i of image input
190  mitk::ImageReadAccessor accessor(input, input->GetSliceData(i));
192  slice->Initialize(input->GetPixelType(), IMAGE_DIMENSION - 1, input->GetDimensions());
193  slice->SetVolume(accessor.GetData());
194 
195  // Transform the slice to matrix
196  m_ImageToOpenCVFilter->SetImage(slice);
197  return m_ImageToOpenCVFilter->GetOpenCVMat();
198 }
199 
201  cv::Mat mat, mitk::Image::Pointer output, unsigned int i) {
202 
203  m_OpenCVToImageFilter->SetOpenCVMat(mat);
204  m_OpenCVToImageFilter->Update();
205  mitk::Image::Pointer slice = m_OpenCVToImageFilter->GetOutput();
206 
207  mitk::ImageReadAccessor accessor(slice);
208  output->SetSlice(accessor.GetData(), i);
209 }
210 
211 // TODO: remove debug messages
213  MITK_INFO << "Start motion compensation.";
214 
215  mitk::Image::Pointer paInput = this->GetInput(0);
216  mitk::Image::Pointer usInput = this->GetInput(1);
217  mitk::Image::Pointer paOutput = this->GetOutput(0);
218  mitk::Image::Pointer usOutput = this->GetOutput(1);
219 
220  // Check that we have two input images with agreeing dimensions
221  this->CheckInput(paInput, usInput);
222 
223  // Check the output images and (re-)initialize, if necessary.
224  this->InitializeOutputIfNecessary(paInput, usInput, paOutput, usOutput);
225 
226  // Set Max and Min of ultrasonic image
227  this->m_MaxValue = usInput->GetStatistics()->GetScalarValueMax();
228  this->m_MinValue = usInput->GetStatistics()->GetScalarValueMin();
229 
230  // m_ImageToOpenCVFilter->SetImage(paInput);
231  this->PerformCorrection(paInput, usInput, paOutput, usOutput);
232 
233  MITK_INFO << "Motion compensation accomplished.";
234 }
mitk::Image::Pointer GetPaOutput()
Wrapper which gets the photoacoustic image out of the correct output.
cv::Mat GetMatrix(const mitk::Image::Pointer input, unsigned int i)
Extract a 2d slice as OpenCV matrix.
void PerformCorrection(mitk::Image::Pointer paInput, mitk::Image::Pointer usInput, mitk::Image::Pointer paOutput, mitk::Image::Pointer usOutput)
This method performs the actual motion compensation.
mitk::Image::Pointer GetPaInput()
Wrapper which gets the photoacoustic image out of the correct input.
#define MITK_INFO
Definition: mitkLogMacros.h:18
#define MITK_ERROR
Definition: mitkLogMacros.h:20
mitk::Image::Pointer GetUsInput()
Wrapper which gets the ultrasonic image out of the correct input.
virtual void SetInput(const InputImageType *image)
cv::Mat FitMatrixToChar(cv::Mat mat)
Rescale matrix such that the values lie between 0 and 255.
void SetPaInput(mitk::Image::Pointer)
Wrapper which sets the photoacoustic image as the correct input.
mitk::Image::Pointer GetUsOutput()
Wrapper which gets the ultrasonic image out of the correct output.
void InitializeOutputIfNecessary(mitk::Image::Pointer paInput, mitk::Image::Pointer usInput, mitk::Image::Pointer paOutput, mitk::Image::Pointer usOutput)
Assure that the output images have the same dimensions as the input images.
cv::Mat ComputeFlowMap(cv::Mat)
Compute the remapping map from an optical flow.
void SetUsInput(mitk::Image::Pointer)
Wrapper which sets the ultrasonic image as the correct input.
void InsertMatrixAsSlice(cv::Mat mat, mitk::Image::Pointer output, unsigned int i)
Insert a OpenCV matrix as a slice into an image.
void InitializeOutput(mitk::Image::Pointer input, mitk::Image::Pointer output)
Copy the image data from the input image to the output image.
void CheckInput(mitk::Image::Pointer paImage, mitk::Image::Pointer usImage)
Validate the input images.
static Pointer New()
MITKMATCHPOINTREGISTRATION_EXPORT ResultImageType::Pointer map(const InputImageType *input, const RegistrationType *registration, bool throwOnOutOfInputAreaError=false, const double &paddingValue=0, const ResultImageGeometryType *resultGeometry=nullptr, bool throwOnMappingError=true, const double &errorValue=0, mitk::ImageMappingInterpolator::Type interpolatorType=mitk::ImageMappingInterpolator::Linear)
InputImageType * GetInput(void)
OutputType * GetOutput()
Get the output data of this image source object.
void GenerateData() override
Apply OpenCV algorithm to compensate motion in a 2d image time series.
ImageReadAccessor class to get locked read access for a particular image part.
const void * GetData() const
Gives const access to the data.