OTB  7.4.0
Orfeo Toolbox
otbDisparityTranslateFilter.hxx
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2005-2020 Centre National d'Etudes Spatiales (CNES)
3  *
4  * This file is part of Orfeo Toolbox
5  *
6  * https://www.orfeo-toolbox.org/
7  *
8  * Licensed under the Apache License, Version 2.0 (the "License");
9  * you may not use this file except in compliance with the License.
10  * You may obtain a copy of the License at
11  *
12  * http://www.apache.org/licenses/LICENSE-2.0
13  *
14  * Unless required by applicable law or agreed to in writing, software
15  * distributed under the License is distributed on an "AS IS" BASIS,
16  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
17  * See the License for the specific language governing permissions and
18  * limitations under the License.
19  */
20 
21 #ifndef otbDisparityTranslateFilter_hxx
22 #define otbDisparityTranslateFilter_hxx
23 
25 
26 #include "itkImageRegionIteratorWithIndex.h"
27 #include "itkImageRegionIterator.h"
28 
29 namespace otb
30 {
31 
32 template <class TDisparityImage, class TGridImage, class TSensorImage, class TMaskImage>
34  : m_NoDataValue(-32768)
35 {
36  // Set the number of inputs (1 moving image by default -> 3 inputs)
37  this->SetNumberOfRequiredInputs(6);
38  this->SetNumberOfRequiredInputs(1);
39 
40  // Set the outputs
41  this->SetNumberOfRequiredOutputs(2);
42  this->SetNthOutput(0, TDisparityImage::New());
43  this->SetNthOutput(1, TDisparityImage::New());
44 }
45 
46 template <class TDisparityImage, class TGridImage, class TSensorImage, class TMaskImage>
48 {
49  // Process object is not const-correct so the const casting is required.
50  this->SetNthInput(0, const_cast<TDisparityImage*>(hmap));
51 }
52 
53 template <class TDisparityImage, class TGridImage, class TSensorImage, class TMaskImage>
55 {
56  // Process object is not const-correct so the const casting is required.
57  this->SetNthInput(1, const_cast<TDisparityImage*>(vmap));
58 }
59 
60 template <class TDisparityImage, class TGridImage, class TSensorImage, class TMaskImage>
62 {
63  // Process object is not const-correct so the const casting is required.
64  this->SetNthInput(2, const_cast<TGridImage*>(lgrid));
65 }
66 
67 template <class TDisparityImage, class TGridImage, class TSensorImage, class TMaskImage>
69 {
70  // Process object is not const-correct so the const casting is required.
71  this->SetNthInput(3, const_cast<TGridImage*>(rgrid));
72 }
73 
74 template <class TDisparityImage, class TGridImage, class TSensorImage, class TMaskImage>
76 {
77  // Process object is not const-correct so the const casting is required.
78  this->SetNthInput(4, const_cast<TMaskImage*>(mask));
79 }
80 
81 template <class TDisparityImage, class TGridImage, class TSensorImage, class TMaskImage>
83 {
84  // Process object is not const-correct so the const casting is required.
85  this->SetNthInput(5, const_cast<TSensorImage*>(left));
86 }
87 
88 template <class TDisparityImage, class TGridImage, class TSensorImage, class TMaskImage>
90 {
91  if (this->GetNumberOfInputs() < 1)
92  {
93  return nullptr;
94  }
95  return static_cast<const TDisparityImage*>(this->itk::ProcessObject::GetInput(0));
96 }
97 
98 template <class TDisparityImage, class TGridImage, class TSensorImage, class TMaskImage>
100 {
101  if (this->GetNumberOfInputs() < 2)
102  {
103  return nullptr;
104  }
105  return static_cast<const TDisparityImage*>(this->itk::ProcessObject::GetInput(1));
106 }
107 
108 template <class TDisparityImage, class TGridImage, class TSensorImage, class TMaskImage>
110 {
111  if (this->GetNumberOfInputs() < 3)
112  {
113  return nullptr;
114  }
115  return static_cast<const TGridImage*>(this->itk::ProcessObject::GetInput(2));
116 }
117 
118 template <class TDisparityImage, class TGridImage, class TSensorImage, class TMaskImage>
120 {
121  if (this->GetNumberOfInputs() < 4)
122  {
123  return nullptr;
124  }
125  return static_cast<const TGridImage*>(this->itk::ProcessObject::GetInput(3));
126 }
127 
128 template <class TDisparityImage, class TGridImage, class TSensorImage, class TMaskImage>
130 {
131  if (this->GetNumberOfInputs() < 5)
132  {
133  return nullptr;
134  }
135  return static_cast<const TMaskImage*>(this->itk::ProcessObject::GetInput(4));
136 }
137 
138 template <class TDisparityImage, class TGridImage, class TSensorImage, class TMaskImage>
140 {
141  if (this->GetNumberOfInputs() < 6)
142  {
143  return nullptr;
144  }
145  return static_cast<const TSensorImage*>(this->itk::ProcessObject::GetInput(5));
146 }
147 
148 template <class TDisparityImage, class TGridImage, class TSensorImage, class TMaskImage>
150 {
151  if (this->GetNumberOfOutputs() < 1)
152  {
153  return nullptr;
154  }
155  return static_cast<TDisparityImage*>(this->itk::ProcessObject::GetOutput(0));
156 }
157 
158 template <class TDisparityImage, class TGridImage, class TSensorImage, class TMaskImage>
160 {
161  if (this->GetNumberOfOutputs() < 2)
162  {
163  return nullptr;
164  }
165  return static_cast<TDisparityImage*>(this->itk::ProcessObject::GetOutput(1));
166 }
167 
168 template <class TDisparityImage, class TGridImage, class TSensorImage, class TMaskImage>
170 {
171  // this->Superclass::GenerateOutputInformation();
172 
173  TDisparityImage* horizOut = this->GetHorizontalDisparityMapOutput();
174  TDisparityImage* vertiOut = this->GetVerticalDisparityMapOutput();
175 
176  const TSensorImage* leftIn = this->GetLeftSensorImageInput();
177 
178  horizOut->CopyInformation(leftIn);
179  vertiOut->CopyInformation(leftIn);
180 
181  // Set the NoData value
182  std::vector<bool> noDataValueAvailable;
183  noDataValueAvailable.push_back(true);
184  std::vector<double> noDataValue;
185  noDataValue.push_back(m_NoDataValue);
186  itk::MetaDataDictionary& dict = horizOut->GetMetaDataDictionary();
187  itk::EncapsulateMetaData<std::vector<bool>>(dict, MetaDataKey::NoDataValueAvailable, noDataValueAvailable);
188  itk::EncapsulateMetaData<std::vector<double>>(dict, MetaDataKey::NoDataValue, noDataValue);
189  dict = vertiOut->GetMetaDataDictionary();
190  itk::EncapsulateMetaData<std::vector<bool>>(dict, MetaDataKey::NoDataValueAvailable, noDataValueAvailable);
191  itk::EncapsulateMetaData<std::vector<double>>(dict, MetaDataKey::NoDataValue, noDataValue);
192 }
193 
194 template <class TDisparityImage, class TGridImage, class TSensorImage, class TMaskImage>
196 {
197  this->Superclass::GenerateInputRequestedRegion();
198 
199  TGridImage* leftGrid = const_cast<TGridImage*>(this->GetInverseEpipolarLeftGrid());
200  TGridImage* rightGrid = const_cast<TGridImage*>(this->GetDirectEpipolarRightGrid());
201 
202  leftGrid->SetRequestedRegionToLargestPossibleRegion();
203  rightGrid->SetRequestedRegionToLargestPossibleRegion();
204 
205  leftGrid->UpdateOutputData();
206  rightGrid->UpdateOutputData();
207 
208  TSensorImage* leftIn = const_cast<TSensorImage*>(this->GetLeftSensorImageInput());
209  RegionType emptyRegion = leftIn->GetLargestPossibleRegion();
210  emptyRegion.SetSize(0, 0);
211  emptyRegion.SetSize(1, 0);
212  leftIn->SetRequestedRegion(emptyRegion);
213 
214  TDisparityImage* horizOut = this->GetHorizontalDisparityMapOutput();
215 
216  TDisparityImage* horizIn = const_cast<TDisparityImage*>(this->GetHorizontalDisparityMapInput());
217  TDisparityImage* vertiIn = const_cast<TDisparityImage*>(this->GetVerticalDisparityMapInput());
218 
219  TMaskImage* maskIn = const_cast<TMaskImage*>(this->GetDisparityMaskInput());
220 
221  RegionType requested = this->GetHorizontalDisparityMapOutput()->GetRequestedRegion();
222  RegionType inputlargest = this->GetHorizontalDisparityMapInput()->GetLargestPossibleRegion();
223  RegionType inputRequested;
224  GridRegionType gridLargest = leftGrid->GetLargestPossibleRegion();
225 
226  IndexType minIndex, maxIndex;
227  // -Wmaybe-uninitialized
228  minIndex.Fill(itk::NumericTraits<IndexValueType>::Zero);
229  maxIndex.Fill(itk::NumericTraits<IndexValueType>::Zero);
230 
231  IndexType corners[4];
232  for (int i = 0; i < 4; i++)
233  corners[i].Fill(itk::NumericTraits<IndexValueType>::Zero);
234 
235  corners[0] = requested.GetIndex();
236  corners[1] = requested.GetIndex();
237  corners[1][0] += static_cast<IndexValueType>(requested.GetSize()[0]) - 1;
238  corners[2] = requested.GetIndex();
239  corners[2][1] += static_cast<IndexValueType>(requested.GetSize()[1]) - 1;
240  corners[3] = requested.GetIndex();
241  corners[3][0] += static_cast<IndexValueType>(requested.GetSize()[0]) - 1;
242  corners[3][1] += static_cast<IndexValueType>(requested.GetSize()[1]) - 1;
243  for (unsigned int k = 0; k < 4; ++k)
244  {
245  PointType pointSensor;
246  horizOut->TransformIndexToPhysicalPoint(corners[k], pointSensor);
247  itk::ContinuousIndex<double, 2> indexGrid;
248  leftGrid->TransformPhysicalPointToContinuousIndex(pointSensor, indexGrid);
249  IndexType ul;
250  ul[0] = static_cast<long>(std::floor(indexGrid[0]));
251  ul[1] = static_cast<long>(std::floor(indexGrid[1]));
252  if (ul[0] < gridLargest.GetIndex()[0])
253  ul[0] = gridLargest.GetIndex()[0];
254  if (ul[1] < gridLargest.GetIndex()[1])
255  ul[1] = gridLargest.GetIndex()[1];
256  if (ul[0] > static_cast<IndexValueType>(gridLargest.GetIndex()[0] + gridLargest.GetSize()[0] - 2))
257  ul[0] = (gridLargest.GetIndex()[0] + gridLargest.GetSize()[0] - 2);
258  if (ul[1] > static_cast<IndexValueType>(gridLargest.GetIndex()[1] + gridLargest.GetSize()[1] - 2))
259  ul[1] = (gridLargest.GetIndex()[1] + gridLargest.GetSize()[1] - 2);
260 
261  IndexType ur = ul;
262  ur[0] += 1;
263  IndexType ll = ul;
264  ll[1] += 1;
265  IndexType lr = ul;
266  lr[0] += 1;
267  lr[1] += 1;
268 
269  double rx = indexGrid[0] - static_cast<double>(ul[0]);
270  double ry = indexGrid[1] - static_cast<double>(ul[1]);
271  PointType pointEpi = pointSensor;
272 
273  pointEpi[0] += (1. - ry) * ((1. - rx) * leftGrid->GetPixel(ul)[0] + rx * leftGrid->GetPixel(ur)[0]) +
274  ry * ((1. - rx) * leftGrid->GetPixel(ll)[0] + rx * leftGrid->GetPixel(lr)[0]);
275  pointEpi[1] += (1. - ry) * ((1. - rx) * leftGrid->GetPixel(ul)[1] + rx * leftGrid->GetPixel(ur)[1]) +
276  ry * ((1. - rx) * leftGrid->GetPixel(ll)[1] + rx * leftGrid->GetPixel(lr)[1]);
277  itk::ContinuousIndex<double, 2> indexEpi;
278 
279  horizIn->TransformPhysicalPointToContinuousIndex(pointEpi, indexEpi);
280  if (k == 0)
281  {
282  minIndex[0] = static_cast<long>(std::floor(indexEpi[0]));
283  minIndex[1] = static_cast<long>(std::floor(indexEpi[1]));
284  maxIndex[0] = static_cast<long>(std::ceil(indexEpi[0]));
285  maxIndex[1] = static_cast<long>(std::ceil(indexEpi[1]));
286  }
287  else
288  {
289  if (minIndex[0] > static_cast<long>(std::floor(indexEpi[0])))
290  minIndex[0] = static_cast<long>(std::floor(indexEpi[0]));
291  if (minIndex[1] > static_cast<long>(std::floor(indexEpi[1])))
292  minIndex[1] = static_cast<long>(std::floor(indexEpi[1]));
293  if (maxIndex[0] < static_cast<long>(std::ceil(indexEpi[0])))
294  maxIndex[0] = static_cast<long>(std::ceil(indexEpi[0]));
295  if (maxIndex[1] < static_cast<long>(std::ceil(indexEpi[1])))
296  maxIndex[1] = static_cast<long>(std::ceil(indexEpi[1]));
297  }
298  }
299 
300  inputRequested.SetIndex(minIndex);
301  inputRequested.SetSize(0, static_cast<unsigned long>(maxIndex[0] - minIndex[0]));
302  inputRequested.SetSize(1, static_cast<unsigned long>(maxIndex[1] - minIndex[1]));
303 
304  if (!inputRequested.Crop(inputlargest))
305  {
306  inputRequested.SetSize(0, 0);
307  inputRequested.SetSize(1, 0);
308  inputRequested.SetIndex(inputlargest.GetIndex());
309  }
310 
311  horizIn->SetRequestedRegion(inputRequested);
312  if (vertiIn)
313  vertiIn->SetRequestedRegion(inputRequested);
314  if (maskIn)
315  maskIn->SetRequestedRegion(inputRequested);
316 }
317 
318 template <class TDisparityImage, class TGridImage, class TSensorImage, class TMaskImage>
320  itk::ThreadIdType itkNotUsed(threadId))
321 {
322  const TGridImage* leftGrid = this->GetInverseEpipolarLeftGrid();
323  const TGridImage* rightGrid = this->GetDirectEpipolarRightGrid();
324 
325  TDisparityImage* horizOut = this->GetHorizontalDisparityMapOutput();
326  TDisparityImage* vertiOut = this->GetVerticalDisparityMapOutput();
327 
328  const TDisparityImage* horizIn = this->GetHorizontalDisparityMapInput();
329  const TDisparityImage* vertiIn = this->GetVerticalDisparityMapInput();
330 
331  const TMaskImage* maskIn = this->GetDisparityMaskInput();
332 
333  GridRegionType leftLargest = leftGrid->GetLargestPossibleRegion();
334  GridRegionType rightLargest = rightGrid->GetLargestPossibleRegion();
335  RegionType buffered = horizIn->GetBufferedRegion();
336 
337  const bool emptyInputRegion = buffered.GetNumberOfPixels() == 0;
338 
339  typedef itk::ImageRegionIteratorWithIndex<TDisparityImage> DispIterator;
340  DispIterator horizIter(horizOut, outputRegionForThread);
341  DispIterator vertiIter(vertiOut, outputRegionForThread);
342 
343  horizIter.GoToBegin();
344  vertiIter.GoToBegin();
345 
346  while (!horizIter.IsAtEnd() && !vertiIter.IsAtEnd())
347  {
348 
349  if (!emptyInputRegion)
350  {
351  PointType pointSensor;
352  horizOut->TransformIndexToPhysicalPoint(horizIter.GetIndex(), pointSensor);
353 
354  itk::ContinuousIndex<double, 2> indexGrid;
355  leftGrid->TransformPhysicalPointToContinuousIndex(pointSensor, indexGrid);
356 
357  // Interpolate in left grid
358  IndexType ul;
359  ul[0] = static_cast<long>(std::floor(indexGrid[0]));
360  ul[1] = static_cast<long>(std::floor(indexGrid[1]));
361  if (ul[0] < leftLargest.GetIndex()[0])
362  ul[0] = leftLargest.GetIndex()[0];
363  if (ul[1] < leftLargest.GetIndex()[1])
364  ul[1] = leftLargest.GetIndex()[1];
365  if (ul[0] > static_cast<IndexValueType>(leftLargest.GetIndex()[0] + leftLargest.GetSize()[0] - 2))
366  ul[0] = (leftLargest.GetIndex()[0] + leftLargest.GetSize()[0] - 2);
367  if (ul[1] > static_cast<IndexValueType>(leftLargest.GetIndex()[1] + leftLargest.GetSize()[1] - 2))
368  ul[1] = (leftLargest.GetIndex()[1] + leftLargest.GetSize()[1] - 2);
369 
370  IndexType ur = ul;
371  ur[0] += 1;
372  IndexType ll = ul;
373  ll[1] += 1;
374  IndexType lr = ul;
375  lr[0] += 1;
376  lr[1] += 1;
377 
378  double rx = indexGrid[0] - static_cast<double>(ul[0]);
379  double ry = indexGrid[1] - static_cast<double>(ul[1]);
380  PointType pointEpi = pointSensor;
381 
382  pointEpi[0] += (1. - ry) * ((1. - rx) * leftGrid->GetPixel(ul)[0] + rx * leftGrid->GetPixel(ur)[0]) +
383  ry * ((1. - rx) * leftGrid->GetPixel(ll)[0] + rx * leftGrid->GetPixel(lr)[0]);
384  pointEpi[1] += (1. - ry) * ((1. - rx) * leftGrid->GetPixel(ul)[1] + rx * leftGrid->GetPixel(ur)[1]) +
385  ry * ((1. - rx) * leftGrid->GetPixel(ll)[1] + rx * leftGrid->GetPixel(lr)[1]);
386 
387  itk::ContinuousIndex<double, 2> indexEpi;
388  horizIn->TransformPhysicalPointToContinuousIndex(pointEpi, indexEpi);
389 
390  // Interpolate in disparity map
391  ul[0] = static_cast<long>(std::floor(indexEpi[0]));
392  ul[1] = static_cast<long>(std::floor(indexEpi[1]));
393  if (ul[0] < buffered.GetIndex()[0])
394  ul[0] = buffered.GetIndex()[0];
395  if (ul[1] < buffered.GetIndex()[1])
396  ul[1] = buffered.GetIndex()[1];
397  if (ul[0] > static_cast<IndexValueType>(buffered.GetIndex()[0] + buffered.GetSize()[0] - 2))
398  ul[0] = (buffered.GetIndex()[0] + buffered.GetSize()[0] - 2);
399  if (ul[1] > static_cast<IndexValueType>(buffered.GetIndex()[1] + buffered.GetSize()[1] - 2))
400  ul[1] = (buffered.GetIndex()[1] + buffered.GetSize()[1] - 2);
401 
402  ur = ul;
403  ur[0] += 1;
404  ll = ul;
405  ll[1] += 1;
406  lr = ul;
407  lr[0] += 1;
408  lr[1] += 1;
409 
410  // check if all corners are valid
411  if (!maskIn || (maskIn && maskIn->GetPixel(ul) > 0 && maskIn->GetPixel(ur) > 0 && maskIn->GetPixel(ll) > 0 && maskIn->GetPixel(lr) > 0))
412  {
413  rx = indexEpi[0] - static_cast<double>(ul[0]);
414  ry = indexEpi[1] - static_cast<double>(ul[1]);
415 
416  itk::ContinuousIndex<double, 2> indexRight(indexEpi);
417 
418  indexRight[0] += (1. - ry) * ((1. - rx) * horizIn->GetPixel(ul) + rx * horizIn->GetPixel(ur)) +
419  ry * ((1. - rx) * horizIn->GetPixel(ll) + rx * horizIn->GetPixel(lr));
420  if (vertiIn)
421  {
422  indexRight[1] += (1. - ry) * ((1. - rx) * vertiIn->GetPixel(ul) + rx * vertiIn->GetPixel(ur)) +
423  ry * ((1. - rx) * vertiIn->GetPixel(ll) + rx * vertiIn->GetPixel(lr));
424  }
425 
426  PointType pointRight;
427  horizIn->TransformContinuousIndexToPhysicalPoint(indexRight, pointRight);
428 
429  itk::ContinuousIndex<double, 2> indexGridRight;
430  rightGrid->TransformPhysicalPointToContinuousIndex(pointRight, indexGridRight);
431 
432  // Interpolate in right grid
433  ul[0] = static_cast<long>(std::floor(indexGridRight[0]));
434  ul[1] = static_cast<long>(std::floor(indexGridRight[1]));
435  if (ul[0] < rightLargest.GetIndex()[0])
436  ul[0] = rightLargest.GetIndex()[0];
437  if (ul[1] < rightLargest.GetIndex()[1])
438  ul[1] = rightLargest.GetIndex()[1];
439  if (ul[0] > static_cast<IndexValueType>(rightLargest.GetIndex()[0] + rightLargest.GetSize()[0] - 2))
440  ul[0] = (rightLargest.GetIndex()[0] + rightLargest.GetSize()[0] - 2);
441  if (ul[1] > static_cast<IndexValueType>(rightLargest.GetIndex()[1] + rightLargest.GetSize()[1] - 2))
442  ul[1] = (rightLargest.GetIndex()[1] + rightLargest.GetSize()[1] - 2);
443 
444  ur = ul;
445  ur[0] += 1;
446  ll = ul;
447  ll[1] += 1;
448  lr = ul;
449  lr[0] += 1;
450  lr[1] += 1;
451 
452  rx = indexGridRight[0] - static_cast<double>(ul[0]);
453  ry = indexGridRight[1] - static_cast<double>(ul[1]);
454 
455  PointType pointSensorRight = pointRight;
456 
457  pointSensorRight[0] += (1. - ry) * ((1. - rx) * rightGrid->GetPixel(ul)[0] + rx * rightGrid->GetPixel(ur)[0]) +
458  ry * ((1. - rx) * rightGrid->GetPixel(ll)[0] + rx * rightGrid->GetPixel(lr)[0]);
459  pointSensorRight[1] += (1. - ry) * ((1. - rx) * rightGrid->GetPixel(ul)[1] + rx * rightGrid->GetPixel(ur)[1]) +
460  ry * ((1. - rx) * rightGrid->GetPixel(ll)[1] + rx * rightGrid->GetPixel(lr)[1]);
461 
462  horizIter.Set(pointSensorRight[0] - pointSensor[0]);
463  vertiIter.Set(pointSensorRight[1] - pointSensor[1]);
464  }
465  else
466  {
467  horizIter.Set(m_NoDataValue);
468  vertiIter.Set(m_NoDataValue);
469  }
470  }
471  else
472  {
473  horizIter.Set(m_NoDataValue);
474  vertiIter.Set(m_NoDataValue);
475  }
476  ++horizIter;
477  ++vertiIter;
478  }
479 }
480 }
481 
482 #endif
void SetLeftSensorImageInput(const TSensorImage *left)
void SetInverseEpipolarLeftGrid(const TGridImage *lgrid)
OTBMetadata_EXPORT char const * NoDataValueAvailable
void ThreadedGenerateData(const RegionType &outputRegionForThread, itk::ThreadIdType threadId) override
const TGridImage * GetDirectEpipolarRightGrid() const
void SetHorizontalDisparityMapInput(const TDisparityImage *hmap)
const TDisparityImage * GetVerticalDisparityMapInput() const
const TDisparityImage * GetHorizontalDisparityMapInput() const
void SetDirectEpipolarRightGrid(const TGridImage *rgrid)
const TMaskImage * GetDisparityMaskInput() const
const TGridImage * GetInverseEpipolarLeftGrid() const
The "otb" namespace contains all Orfeo Toolbox (OTB) classes.
void SetDisparityMaskInput(const TMaskImage *mask)
OTBMetadata_EXPORT char const * NoDataValue
const TSensorImage * GetLeftSensorImageInput() const
DispMapType::IndexValueType IndexValueType
void SetVerticalDisparityMapInput(const TDisparityImage *vmap)