OTB  10.0.0
Orfeo Toolbox
otbDisparityMapToDEMFilter.hxx
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2005-2024 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 otbDisparityMapToDEMFilter_hxx
22 #define otbDisparityMapToDEMFilter_hxx
23 
24 #include "otbDisparityMapToDEMFilter.h"
25 #include "itkImageRegionConstIteratorWithIndex.h"
26 #include "itkImageRegionIterator.h"
27 
28 namespace otb
29 {
30 
31 template <class TDisparityImage, class TInputImage, class TOutputDEMImage, class TEpipolarGridImage, class TMaskImage>
32 DisparityMapToDEMFilter<TDisparityImage, TInputImage, TOutputDEMImage, TEpipolarGridImage, TMaskImage>::DisparityMapToDEMFilter()
33 {
34  this->DynamicMultiThreadingOff();
35  // Set the number of inputs
36  this->SetNumberOfRequiredInputs(7);
37  this->SetNumberOfRequiredInputs(1);
38 
39  // Set the outputs
40  this->SetNumberOfRequiredOutputs(1);
41  this->SetNthOutput(0, TOutputDEMImage::New());
42 
43  // Default DEM reconstruction parameters
44  m_ElevationMin = 0.0;
45  m_ElevationMax = 100.0;
46  m_DEMGridStep = 10.0;
47 
48  m_InputSplitter = SplitterType::New();
49  m_UsedInputSplits = 1;
50 }
51 
52 template <class TDisparityImage, class TInputImage, class TOutputDEMImage, class TEpipolarGridImage, class TMaskImage>
53 DisparityMapToDEMFilter<TDisparityImage, TInputImage, TOutputDEMImage, TEpipolarGridImage, TMaskImage>::~DisparityMapToDEMFilter()
54 {
55 }
56 
57 template <class TDisparityImage, class TInputImage, class TOutputDEMImage, class TEpipolarGridImage, class TMaskImage>
58 void DisparityMapToDEMFilter<TDisparityImage, TInputImage, TOutputDEMImage, TEpipolarGridImage, TMaskImage>::SetHorizontalDisparityMapInput(
59  const TDisparityImage* hmap)
60 {
61  // Process object is not const-correct so the const casting is required.
62  this->SetNthInput(0, const_cast<TDisparityImage*>(hmap));
63 }
64 
65 template <class TDisparityImage, class TInputImage, class TOutputDEMImage, class TEpipolarGridImage, class TMaskImage>
66 void DisparityMapToDEMFilter<TDisparityImage, TInputImage, TOutputDEMImage, TEpipolarGridImage, TMaskImage>::SetVerticalDisparityMapInput(
67  const TDisparityImage* vmap)
68 {
69  // Process object is not const-correct so the const casting is required.
70  this->SetNthInput(1, const_cast<TDisparityImage*>(vmap));
71 }
72 
73 template <class TDisparityImage, class TInputImage, class TOutputDEMImage, class TEpipolarGridImage, class TMaskImage>
74 void DisparityMapToDEMFilter<TDisparityImage, TInputImage, TOutputDEMImage, TEpipolarGridImage, TMaskImage>::SetLeftInput(const TInputImage* image)
75 {
76  // Process object is not const-correct so the const casting is required.
77  this->SetNthInput(2, const_cast<TInputImage*>(image));
78 }
79 
80 template <class TDisparityImage, class TInputImage, class TOutputDEMImage, class TEpipolarGridImage, class TMaskImage>
81 void DisparityMapToDEMFilter<TDisparityImage, TInputImage, TOutputDEMImage, TEpipolarGridImage, TMaskImage>::SetRightInput(const TInputImage* image)
82 {
83  // Process object is not const-correct so the const casting is required.
84  this->SetNthInput(3, const_cast<TInputImage*>(image));
85 }
86 
87 template <class TDisparityImage, class TInputImage, class TOutputDEMImage, class TEpipolarGridImage, class TMaskImage>
88 void DisparityMapToDEMFilter<TDisparityImage, TInputImage, TOutputDEMImage, TEpipolarGridImage, TMaskImage>::SetLeftEpipolarGridInput(
89  const TEpipolarGridImage* grid)
90 {
91  // Process object is not const-correct so the const casting is required.
92  this->SetNthInput(4, const_cast<TEpipolarGridImage*>(grid));
93 }
94 
95 template <class TDisparityImage, class TInputImage, class TOutputDEMImage, class TEpipolarGridImage, class TMaskImage>
96 void DisparityMapToDEMFilter<TDisparityImage, TInputImage, TOutputDEMImage, TEpipolarGridImage, TMaskImage>::SetRightEpipolarGridInput(
97  const TEpipolarGridImage* grid)
98 {
99  // Process object is not const-correct so the const casting is required.
100  this->SetNthInput(5, const_cast<TEpipolarGridImage*>(grid));
101 }
102 
103 template <class TDisparityImage, class TInputImage, class TOutputDEMImage, class TEpipolarGridImage, class TMaskImage>
104 void DisparityMapToDEMFilter<TDisparityImage, TInputImage, TOutputDEMImage, TEpipolarGridImage, TMaskImage>::SetDisparityMaskInput(const TMaskImage* mask)
105 {
106  // Process object is not const-correct so the const casting is required.
107  this->SetNthInput(6, const_cast<TMaskImage*>(mask));
108 }
109 
110 template <class TDisparityImage, class TInputImage, class TOutputDEMImage, class TEpipolarGridImage, class TMaskImage>
111 const TDisparityImage*
112 DisparityMapToDEMFilter<TDisparityImage, TInputImage, TOutputDEMImage, TEpipolarGridImage, TMaskImage>::GetHorizontalDisparityMapInput() const
113 {
114  if (this->GetNumberOfInputs() < 1)
115  {
116  return nullptr;
117  }
118  return static_cast<const TDisparityImage*>(this->itk::ProcessObject::GetInput(0));
119 }
120 
121 template <class TDisparityImage, class TInputImage, class TOutputDEMImage, class TEpipolarGridImage, class TMaskImage>
122 const TDisparityImage*
123 DisparityMapToDEMFilter<TDisparityImage, TInputImage, TOutputDEMImage, TEpipolarGridImage, TMaskImage>::GetVerticalDisparityMapInput() const
124 {
125  if (this->GetNumberOfInputs() < 2)
126  {
127  return nullptr;
128  }
129  return static_cast<const TDisparityImage*>(this->itk::ProcessObject::GetInput(1));
130 }
131 
132 template <class TDisparityImage, class TInputImage, class TOutputDEMImage, class TEpipolarGridImage, class TMaskImage>
133 const TInputImage* DisparityMapToDEMFilter<TDisparityImage, TInputImage, TOutputDEMImage, TEpipolarGridImage, TMaskImage>::GetLeftInput() const
134 {
135  if (this->GetNumberOfInputs() < 3)
136  {
137  return nullptr;
138  }
139  return static_cast<const TInputImage*>(this->itk::ProcessObject::GetInput(2));
140 }
141 
142 template <class TDisparityImage, class TInputImage, class TOutputDEMImage, class TEpipolarGridImage, class TMaskImage>
143 const TInputImage* DisparityMapToDEMFilter<TDisparityImage, TInputImage, TOutputDEMImage, TEpipolarGridImage, TMaskImage>::GetRightInput() const
144 {
145  if (this->GetNumberOfInputs() < 4)
146  {
147  return nullptr;
148  }
149  return static_cast<const TInputImage*>(this->itk::ProcessObject::GetInput(3));
150 }
151 
152 template <class TDisparityImage, class TInputImage, class TOutputDEMImage, class TEpipolarGridImage, class TMaskImage>
153 const TEpipolarGridImage*
154 DisparityMapToDEMFilter<TDisparityImage, TInputImage, TOutputDEMImage, TEpipolarGridImage, TMaskImage>::GetLeftEpipolarGridInput() const
155 {
156  if (this->GetNumberOfInputs() < 5)
157  {
158  return nullptr;
159  }
160  return static_cast<const TEpipolarGridImage*>(this->itk::ProcessObject::GetInput(4));
161 }
162 
163 template <class TDisparityImage, class TInputImage, class TOutputDEMImage, class TEpipolarGridImage, class TMaskImage>
164 const TEpipolarGridImage*
165 DisparityMapToDEMFilter<TDisparityImage, TInputImage, TOutputDEMImage, TEpipolarGridImage, TMaskImage>::GetRightEpipolarGridInput() const
166 {
167  if (this->GetNumberOfInputs() < 6)
168  {
169  return nullptr;
170  }
171  return static_cast<const TEpipolarGridImage*>(this->itk::ProcessObject::GetInput(5));
172 }
173 
174 template <class TDisparityImage, class TInputImage, class TOutputDEMImage, class TEpipolarGridImage, class TMaskImage>
175 const TMaskImage* DisparityMapToDEMFilter<TDisparityImage, TInputImage, TOutputDEMImage, TEpipolarGridImage, TMaskImage>::GetDisparityMaskInput() const
176 {
177  if (this->GetNumberOfInputs() < 7)
178  {
179  return nullptr;
180  }
181  return static_cast<const TMaskImage*>(this->itk::ProcessObject::GetInput(6));
182 }
183 
184 template <class TDisparityImage, class TInputImage, class TOutputDEMImage, class TEpipolarGridImage, class TMaskImage>
185 const TOutputDEMImage* DisparityMapToDEMFilter<TDisparityImage, TInputImage, TOutputDEMImage, TEpipolarGridImage, TMaskImage>::GetDEMOutput() const
186 {
187  if (this->GetNumberOfOutputs() < 1)
188  {
189  return 0;
190  }
191  return static_cast<const TOutputDEMImage*>(this->itk::ProcessObject::GetOutput(0));
192 }
193 
194 template <class TDisparityImage, class TInputImage, class TOutputDEMImage, class TEpipolarGridImage, class TMaskImage>
195 TOutputDEMImage* DisparityMapToDEMFilter<TDisparityImage, TInputImage, TOutputDEMImage, TEpipolarGridImage, TMaskImage>::GetDEMOutput()
196 {
197  if (this->GetNumberOfOutputs() < 1)
198  {
199  return nullptr;
200  }
201  return static_cast<TOutputDEMImage*>(this->itk::ProcessObject::GetOutput(0));
202 }
203 
204 template <class TDisparityImage, class TInputImage, class TOutputDEMImage, class TEpipolarGridImage, class TMaskImage>
205 void DisparityMapToDEMFilter<TDisparityImage, TInputImage, TOutputDEMImage, TEpipolarGridImage, TMaskImage>::GenerateOutputInformation()
206 {
207  // Get common area between images
208  const TInputImage* leftImgPtr = this->GetLeftInput();
209  const TInputImage* rightImgPtr = this->GetRightInput();
210  TOutputDEMImage* outputPtr = this->GetDEMOutput();
211 
212  // Set-up a transform to use the DEMHandler
213  typedef otb::GenericRSTransform<> RSTransform2DType;
214  RSTransform2DType::Pointer leftToGroundTransform = RSTransform2DType::New();
215  leftToGroundTransform->SetInputImageMetadata(&(leftImgPtr->GetImageMetadata()));
216  leftToGroundTransform->InstantiateTransform();
217 
218  RSTransform2DType::Pointer rightToGroundTransform = RSTransform2DType::New();
219  rightToGroundTransform->SetInputImageMetadata(&(rightImgPtr->GetImageMetadata()));
220  rightToGroundTransform->InstantiateTransform();
221 
222  // left image
223  typename SensorImageType::SizeType inputSize = leftImgPtr->GetLargestPossibleRegion().GetSize();
224  typename SensorImageType::PointType ulp, urp, llp, lrp;
225  itk::ContinuousIndex<double, 2> ul(leftImgPtr->GetLargestPossibleRegion().GetIndex());
226  ul[0] += -0.5;
227  ul[1] += -0.5;
228 
229  itk::ContinuousIndex<double, 2> ur(ul);
230  itk::ContinuousIndex<double, 2> lr(ul);
231  itk::ContinuousIndex<double, 2> ll(ul);
232  ur[0] += static_cast<double>(inputSize[0]);
233  lr[0] += static_cast<double>(inputSize[0]);
234  lr[1] += static_cast<double>(inputSize[1]);
235  ll[1] += static_cast<double>(inputSize[1]);
236 
237  leftImgPtr->TransformContinuousIndexToPhysicalPoint(ul, ulp);
238  leftImgPtr->TransformContinuousIndexToPhysicalPoint(ur, urp);
239  leftImgPtr->TransformContinuousIndexToPhysicalPoint(ll, llp);
240  leftImgPtr->TransformContinuousIndexToPhysicalPoint(lr, lrp);
241 
242  RSTransform2DType::OutputPointType left_ul, left_ur, left_ll, left_lr;
243  left_ul = leftToGroundTransform->TransformPoint(ulp);
244  left_ur = leftToGroundTransform->TransformPoint(urp);
245  left_ll = leftToGroundTransform->TransformPoint(llp);
246  left_lr = leftToGroundTransform->TransformPoint(lrp);
247 
248  // right image
249  inputSize = rightImgPtr->GetLargestPossibleRegion().GetSize();
250  ul = rightImgPtr->GetLargestPossibleRegion().GetIndex();
251  ul[0] += -0.5;
252  ul[1] += -0.5;
253  ur = ul;
254  lr = ul;
255  ll = ul;
256  ur[0] += static_cast<double>(inputSize[0]);
257  lr[0] += static_cast<double>(inputSize[0]);
258  lr[1] += static_cast<double>(inputSize[1]);
259  ll[1] += static_cast<double>(inputSize[1]);
260 
261  rightImgPtr->TransformContinuousIndexToPhysicalPoint(ul, ulp);
262  rightImgPtr->TransformContinuousIndexToPhysicalPoint(ur, urp);
263  rightImgPtr->TransformContinuousIndexToPhysicalPoint(ll, llp);
264  rightImgPtr->TransformContinuousIndexToPhysicalPoint(lr, lrp);
265 
266  RSTransform2DType::OutputPointType right_ul, right_ur, right_lr, right_ll;
267  right_ul = rightToGroundTransform->TransformPoint(ulp);
268  right_ur = rightToGroundTransform->TransformPoint(urp);
269  right_ll = rightToGroundTransform->TransformPoint(llp);
270  right_lr = rightToGroundTransform->TransformPoint(lrp);
271 
272  double left_xmin = std::min(std::min(std::min(left_ul[0], left_ur[0]), left_lr[0]), left_ll[0]);
273  double left_xmax = std::max(std::max(std::max(left_ul[0], left_ur[0]), left_lr[0]), left_ll[0]);
274  double left_ymin = std::min(std::min(std::min(left_ul[1], left_ur[1]), left_lr[1]), left_ll[1]);
275  double left_ymax = std::max(std::max(std::max(left_ul[1], left_ur[1]), left_lr[1]), left_ll[1]);
276 
277  double right_xmin = std::min(std::min(std::min(right_ul[0], right_ur[0]), right_lr[0]), right_ll[0]);
278  double right_xmax = std::max(std::max(std::max(right_ul[0], right_ur[0]), right_lr[0]), right_ll[0]);
279  double right_ymin = std::min(std::min(std::min(right_ul[1], right_ur[1]), right_lr[1]), right_ll[1]);
280  double right_ymax = std::max(std::max(std::max(right_ul[1], right_ur[1]), right_lr[1]), right_ll[1]);
281 
282  double box_xmin = std::max(left_xmin, right_xmin);
283  double box_xmax = std::min(left_xmax, right_xmax);
284  double box_ymin = std::max(left_ymin, right_ymin);
285  double box_ymax = std::min(left_ymax, right_ymax);
286 
287  if (box_xmin >= box_xmax || box_ymin >= box_ymax)
288  {
289  itkExceptionMacro(<< "Wrong reconstruction area, images don't overlap, check image corners");
290  }
291 
292  // Compute step :
293  // TODO : use a clean RS transform instead
294  typename TOutputDEMImage::SpacingType outSpacing;
295  outSpacing[0] = 57.295779513 * m_DEMGridStep / (6378137.0 * std::cos((box_ymin + box_ymax) * 0.5 * 0.01745329251));
296  outSpacing[1] = -57.295779513 * m_DEMGridStep / 6378137.0;
297  outputPtr->SetSignedSpacing(outSpacing);
298 
299  // Choose origin
300  typename TOutputDEMImage::PointType outOrigin;
301  outOrigin[0] = box_xmin + 0.5 * outSpacing[0];
302  outOrigin[1] = box_ymax + 0.5 * outSpacing[1];
303  outputPtr->SetOrigin(outOrigin);
304 
305  // Compute output size
306  typename DEMImageType::RegionType outRegion;
307  outRegion.SetIndex(0, 0);
308  outRegion.SetIndex(1, 0);
309  outRegion.SetSize(0, static_cast<unsigned int>((box_xmax - box_xmin) / std::abs(outSpacing[0])));
310  outRegion.SetSize(1, static_cast<unsigned int>((box_ymax - box_ymin) / std::abs(outSpacing[1])));
311 
312  outputPtr->SetLargestPossibleRegion(outRegion);
313  outputPtr->SetNumberOfComponentsPerPixel(1);
314 }
315 
316 template <class TDisparityImage, class TInputImage, class TOutputDEMImage, class TEpipolarGridImage, class TMaskImage>
317 void DisparityMapToDEMFilter<TDisparityImage, TInputImage, TOutputDEMImage, TEpipolarGridImage, TMaskImage>::GenerateInputRequestedRegion()
318 {
319  // For the epi grid : generate full buffer here !
320  TEpipolarGridImage* leftGrid = const_cast<TEpipolarGridImage*>(this->GetLeftEpipolarGridInput());
321  TEpipolarGridImage* rightGrid = const_cast<TEpipolarGridImage*>(this->GetRightEpipolarGridInput());
322 
323  leftGrid->SetRequestedRegionToLargestPossibleRegion();
324  rightGrid->SetRequestedRegionToLargestPossibleRegion();
325 
326  leftGrid->UpdateOutputData();
327  rightGrid->UpdateOutputData();
328 
329  // For the input left-right images
330  // -> No need, only use metadata and keywordlist
331  const TInputImage* leftSensor = this->GetLeftInput();
332  TOutputDEMImage* outputDEM = this->GetDEMOutput();
333 
334  typename DEMImageType::RegionType outRegion = outputDEM->GetRequestedRegion();
335  typename DEMImageType::PointType outOrigin = outputDEM->GetOrigin();
336  typename DEMImageType::SpacingType outSpacing = outputDEM->GetSignedSpacing();
337 
338  RSTransformType::Pointer groundToLeftTransform = RSTransformType::New();
339  groundToLeftTransform->SetOutputImageMetadata(&(leftSensor->GetImageMetadata()));
340  groundToLeftTransform->InstantiateTransform();
341 
342  // For the disparity maps and mask
343  // Iterate over OutputRequestedRegion corners for elevation min and max
344  TDisparityImage* horizDisp = const_cast<TDisparityImage*>(this->GetHorizontalDisparityMapInput());
345  TDisparityImage* vertiDisp = const_cast<TDisparityImage*>(this->GetVerticalDisparityMapInput());
346  TMaskImage* maskDisp = const_cast<TMaskImage*>(this->GetDisparityMaskInput());
347 
348  // We impose that both disparity map inputs have the same size
349  if (horizDisp->GetLargestPossibleRegion() != vertiDisp->GetLargestPossibleRegion())
350  {
351  itkExceptionMacro(<< "Horizontal and vertical disparity maps do not have the same size ! Horizontal largest region: "
352  << horizDisp->GetLargestPossibleRegion() << ", vertical largest region: " << vertiDisp->GetLargestPossibleRegion());
353  }
354 
355  if (maskDisp && horizDisp->GetLargestPossibleRegion() != maskDisp->GetLargestPossibleRegion())
356  {
357  itkExceptionMacro(<< "Disparity map and mask do not have the same size ! Map region : " << horizDisp->GetLargestPossibleRegion()
358  << ", mask region : " << maskDisp->GetLargestPossibleRegion());
359  }
360 
361  // up left at elevation min
362  TDPointType corners[8];
363  corners[0][0] = outOrigin[0] + outSpacing[0] * (-0.5 + static_cast<double>(outRegion.GetIndex(0)));
364  corners[0][1] = outOrigin[1] + outSpacing[1] * (-0.5 + static_cast<double>(outRegion.GetIndex(1)));
365  corners[0][2] = m_ElevationMin;
366  // up left at elevation max
367  corners[1][0] = corners[0][0];
368  corners[1][1] = corners[0][1];
369  corners[1][2] = m_ElevationMax;
370  // up right at elevation min
371  corners[2][0] = corners[0][0] + outSpacing[0] * static_cast<double>(outRegion.GetSize(0));
372  corners[2][1] = corners[0][1];
373  corners[2][2] = m_ElevationMin;
374  // up right at elevation max
375  corners[3][0] = corners[2][0];
376  corners[3][1] = corners[2][1];
377  corners[3][2] = m_ElevationMax;
378  // low right at elevation min
379  corners[4][0] = corners[0][0] + outSpacing[0] * static_cast<double>(outRegion.GetSize(0));
380  corners[4][1] = corners[0][1] + outSpacing[1] * static_cast<double>(outRegion.GetSize(1));
381  corners[4][2] = m_ElevationMin;
382  // low right at elevation max
383  corners[5][0] = corners[4][0];
384  corners[5][1] = corners[4][1];
385  corners[5][2] = m_ElevationMax;
386  // low left at elevation min
387  corners[6][0] = corners[0][0];
388  corners[6][1] = corners[0][1] + outSpacing[1] * static_cast<double>(outRegion.GetSize(1));
389  corners[6][2] = m_ElevationMin;
390  // low left at elevation max
391  corners[7][0] = corners[6][0];
392  corners[7][1] = corners[6][1];
393  corners[7][2] = m_ElevationMax;
394 
395  double x_loc, y_loc;
396  double a_grid, b_grid;
397  double u_loc[2];
398  double v_loc[2];
399  double det;
400  typename GridImageType::IndexType gridIndex;
401  typename GridImageType::IndexType gridIndexU;
402  typename GridImageType::IndexType gridIndexV;
403  typename GridImageType::PointType gridPoint;
404  typename GridImageType::PointType gridPointU;
405  typename GridImageType::PointType gridPointV;
406  typename GridImageType::PixelType origLoc(2);
407  typename GridImageType::PixelType uUnitLoc(2);
408  typename GridImageType::PixelType vUnitLoc(2);
409 
410  typename GridImageType::RegionType gridRegion = leftGrid->GetLargestPossibleRegion();
411  itk::ContinuousIndex<double, 2> gridContiIndex;
412  typename DisparityMapType::PointType epiPosition;
413  itk::ContinuousIndex<double, 2> epiContiIndex;
414  double epiIndexMin[2];
415  double epiIndexMax[2];
416  int maxGridIndex[2];
417  int minGridIndex[2];
418 
419  epiIndexMin[0] = itk::NumericTraits<double>::Zero;
420  epiIndexMin[1] = itk::NumericTraits<double>::Zero;
421 
422  epiIndexMax[0] = itk::NumericTraits<double>::Zero;
423  epiIndexMax[1] = itk::NumericTraits<double>::Zero;
424 
425  maxGridIndex[0] = static_cast<int>(gridRegion.GetIndex(0) + gridRegion.GetSize(0)) - 2;
426  maxGridIndex[1] = static_cast<int>(gridRegion.GetIndex(1) + gridRegion.GetSize(1)) - 2;
427  minGridIndex[0] = static_cast<int>(gridRegion.GetIndex(0));
428  minGridIndex[1] = static_cast<int>(gridRegion.GetIndex(1));
429 
430  for (unsigned int k = 0; k < 8; k++)
431  {
432  // compute left image coordinate
433  TDPointType tmpSensor = groundToLeftTransform->TransformPoint(corners[k]);
434 
435  // compute epipolar position
436  gridIndex[0] = minGridIndex[0];
437  gridIndex[1] = minGridIndex[1];
438 
439  // we assume 3 iterations are enough to find the 4 surrounding pixels
440  for (unsigned int s = 0; s < 3; s++)
441  {
442  gridIndexU[0] = gridIndex[0] + 1;
443  gridIndexU[1] = gridIndex[1];
444 
445  gridIndexV[0] = gridIndex[0];
446  gridIndexV[1] = gridIndex[1] + 1;
447 
448  leftGrid->TransformIndexToPhysicalPoint(gridIndex, gridPoint);
449  leftGrid->TransformIndexToPhysicalPoint(gridIndexU, gridPointU);
450  leftGrid->TransformIndexToPhysicalPoint(gridIndexV, gridPointV);
451 
452  origLoc[0] = (leftGrid->GetPixel(gridIndex))[0] + gridPoint[0];
453  origLoc[1] = (leftGrid->GetPixel(gridIndex))[1] + gridPoint[1];
454  uUnitLoc[0] = (leftGrid->GetPixel(gridIndexU))[0] + gridPointU[0];
455  uUnitLoc[1] = (leftGrid->GetPixel(gridIndexU))[1] + gridPointU[1];
456  vUnitLoc[0] = (leftGrid->GetPixel(gridIndexV))[0] + gridPointV[0];
457  vUnitLoc[1] = (leftGrid->GetPixel(gridIndexV))[1] + gridPointV[1];
458 
459  u_loc[0] = static_cast<double>(uUnitLoc[0] - origLoc[0]);
460  u_loc[1] = static_cast<double>(uUnitLoc[1] - origLoc[1]);
461 
462  v_loc[0] = static_cast<double>(vUnitLoc[0] - origLoc[0]);
463  v_loc[1] = static_cast<double>(vUnitLoc[1] - origLoc[1]);
464 
465  det = u_loc[0] * v_loc[1] - v_loc[0] * u_loc[1];
466 
467  x_loc = static_cast<double>(tmpSensor[0]) - static_cast<double>(origLoc[0]);
468  y_loc = static_cast<double>(tmpSensor[1]) - static_cast<double>(origLoc[1]);
469 
470  a_grid = (x_loc * v_loc[1] - y_loc * v_loc[0]) / det;
471  b_grid = (y_loc * u_loc[0] - x_loc * u_loc[1]) / det;
472 
473  gridContiIndex[0] = static_cast<double>(gridIndex[0]) + a_grid;
474  gridContiIndex[1] = static_cast<double>(gridIndex[1]) + b_grid;
475 
476  leftGrid->TransformContinuousIndexToPhysicalPoint(gridContiIndex, epiPosition);
477 
478  horizDisp->TransformPhysicalPointToContinuousIndex(epiPosition, epiContiIndex);
479 
480  if (0.0 < a_grid && a_grid < 1.0 && 0.0 < b_grid && b_grid < 1.0)
481  {
482  // The four nearest positions in epipolar grid have been found : stop search
483  break;
484  }
485  else
486  {
487  // Shift the gridIndex
488  int a_grid_int = static_cast<int>(gridIndex[0]) + static_cast<int>(std::floor(a_grid));
489  int b_grid_int = static_cast<int>(gridIndex[1]) + static_cast<int>(std::floor(b_grid));
490 
491  if (a_grid_int < minGridIndex[0])
492  a_grid_int = minGridIndex[0];
493  if (a_grid_int > maxGridIndex[0])
494  a_grid_int = maxGridIndex[0];
495  if (b_grid_int < minGridIndex[1])
496  b_grid_int = minGridIndex[1];
497  if (b_grid_int > maxGridIndex[1])
498  b_grid_int = maxGridIndex[1];
499 
500  gridIndex[0] = a_grid_int;
501  gridIndex[1] = b_grid_int;
502  }
503  }
504 
505  if (k == 0)
506  {
507  epiIndexMin[0] = epiContiIndex[0];
508  epiIndexMin[1] = epiContiIndex[1];
509  epiIndexMax[0] = epiContiIndex[0];
510  epiIndexMax[1] = epiContiIndex[1];
511  }
512  else
513  {
514  if (epiContiIndex[0] < epiIndexMin[0])
515  epiIndexMin[0] = epiContiIndex[0];
516  if (epiContiIndex[1] < epiIndexMin[1])
517  epiIndexMin[1] = epiContiIndex[1];
518  if (epiIndexMax[0] < epiContiIndex[0])
519  epiIndexMax[0] = epiContiIndex[0];
520  if (epiIndexMax[1] < epiContiIndex[1])
521  epiIndexMax[1] = epiContiIndex[1];
522  }
523  }
524 
525  typename DisparityMapType::RegionType inputDisparityRegion;
526  inputDisparityRegion.SetIndex(0, static_cast<int>(std::floor(epiIndexMin[0] + 0.5)));
527  inputDisparityRegion.SetIndex(1, static_cast<int>(std::floor(epiIndexMin[1] + 0.5)));
528  inputDisparityRegion.SetSize(0, 1 + static_cast<unsigned int>(std::floor(epiIndexMax[0] - epiIndexMin[0] + 0.5)));
529  inputDisparityRegion.SetSize(1, 1 + static_cast<unsigned int>(std::floor(epiIndexMax[1] - epiIndexMin[1] + 0.5)));
530 
531  // crop the disparity region at the largest possible region
532  if (inputDisparityRegion.Crop(horizDisp->GetLargestPossibleRegion()))
533  {
534  horizDisp->SetRequestedRegion(inputDisparityRegion);
535  vertiDisp->SetRequestedRegion(inputDisparityRegion);
536 
537  if (maskDisp)
538  {
539  maskDisp->SetRequestedRegion(inputDisparityRegion);
540  }
541  }
542  else
543  {
544  // use empty region
545  typename DisparityMapType::RegionType emptyRegion;
546 
547  emptyRegion.SetIndex(horizDisp->GetLargestPossibleRegion().GetIndex());
548  emptyRegion.SetSize(0, 0);
549  emptyRegion.SetSize(1, 0);
550 
551  horizDisp->SetRequestedRegion(emptyRegion);
552  vertiDisp->SetRequestedRegion(emptyRegion);
553  if (maskDisp)
554  {
555  maskDisp->SetRequestedRegion(emptyRegion);
556  }
557 
558  // no exception : it can happen in this filter
559  if (0)
560  {
561  // Couldn't crop the region (requested region is outside the largest
562  // possible region). Throw an exception.
563  // store what we tried to request (prior to trying to crop)
564  horizDisp->SetRequestedRegion(inputDisparityRegion);
565  vertiDisp->SetRequestedRegion(inputDisparityRegion);
566 
567  // build an exception
568  itk::InvalidRequestedRegionError e(__FILE__, __LINE__);
569  std::ostringstream msg;
570  msg << this->GetNameOfClass() << "::GenerateInputRequestedRegion()";
571  e.SetLocation(msg.str());
572  e.SetDescription("Requested region is (at least partially) outside the largest possible region of disparity map.");
573  e.SetDataObject(horizDisp);
574  throw e;
575  }
576  }
577 }
578 
579 template <class TDisparityImage, class TInputImage, class TOutputDEMImage, class TEpipolarGridImage, class TMaskImage>
580 void DisparityMapToDEMFilter<TDisparityImage, TInputImage, TOutputDEMImage, TEpipolarGridImage, TMaskImage>::BeforeThreadedGenerateData()
581 {
582  const TDisparityImage* horizDisp = this->GetHorizontalDisparityMapInput();
583  const TInputImage* leftSensor = this->GetLeftInput();
584  const TInputImage* rightSensor = this->GetRightInput();
585 
586  const TOutputDEMImage* outputDEM = this->GetDEMOutput();
587 
588  typename DisparityMapType::RegionType requestedRegion = horizDisp->GetRequestedRegion();
589 
590  m_UsedInputSplits = m_InputSplitter->GetNumberOfSplits(requestedRegion, this->GetNumberOfWorkUnits());
591 
592  m_LeftToGroundTransform = RSTransformType::New();
593  m_RightToGroundTransform = RSTransformType::New();
594 
595  m_LeftToGroundTransform->SetInputImageMetadata(&(leftSensor->GetImageMetadata()));
596  m_RightToGroundTransform->SetInputImageMetadata(&(rightSensor->GetImageMetadata()));
597 
598  m_LeftToGroundTransform->InstantiateTransform();
599  m_RightToGroundTransform->InstantiateTransform();
600 
601  // ensure empty regions are not processed
602  if (requestedRegion.GetSize(0) == 0 && requestedRegion.GetSize(1) == 0)
603  {
604  m_UsedInputSplits = 0;
605  }
606 
607  if (m_UsedInputSplits <= static_cast<unsigned int>(this->GetNumberOfWorkUnits()))
608  {
609  m_TempDEMRegions.clear();
610 
611  for (unsigned int i = 0; i < m_UsedInputSplits; i++)
612  {
613  typename DEMImageType::Pointer tmpImg = TOutputDEMImage::New();
614  tmpImg->SetNumberOfComponentsPerPixel(1);
615  tmpImg->SetRegions(outputDEM->GetRequestedRegion());
616  tmpImg->Allocate();
617 
618  typename DEMImageType::PixelType minElevation = static_cast<typename DEMImageType::PixelType>(m_ElevationMin);
619  tmpImg->FillBuffer(minElevation);
620 
621  m_TempDEMRegions.push_back(tmpImg);
622  }
623  }
624  else
625  {
626  itkExceptionMacro(<< "Wrong number of splits for input multithreading : " << m_UsedInputSplits);
627  }
628 }
629 
630 template <class TDisparityImage, class TInputImage, class TOutputDEMImage, class TEpipolarGridImage, class TMaskImage>
631 void DisparityMapToDEMFilter<TDisparityImage, TInputImage, TOutputDEMImage, TEpipolarGridImage, TMaskImage>::ThreadedGenerateData(
632  const RegionType& itkNotUsed(outputRegionForThread), itk::ThreadIdType threadId)
633 {
634  const TDisparityImage* horizDisp = this->GetHorizontalDisparityMapInput();
635  const TDisparityImage* vertiDisp = this->GetVerticalDisparityMapInput();
636 
637  const TMaskImage* disparityMask = this->GetDisparityMaskInput();
638 
639  const TOutputDEMImage* outputDEM = this->GetDEMOutput();
640 
641  const TEpipolarGridImage* leftGrid = this->GetLeftEpipolarGridInput();
642  const TEpipolarGridImage* rightGrid = this->GetRightEpipolarGridInput();
643 
644  typename TEpipolarGridImage::RegionType gridRegion = leftGrid->GetLargestPossibleRegion();
645 
646  TOutputDEMImage* tmpDEM = nullptr;
647  typename TOutputDEMImage::RegionType outputRequestedRegion = outputDEM->GetRequestedRegion();
648 
649  typename TDisparityImage::RegionType disparityRegion;
650  if (static_cast<unsigned int>(threadId) < m_UsedInputSplits)
651  {
652  disparityRegion = m_InputSplitter->GetSplit(threadId, m_UsedInputSplits, horizDisp->GetRequestedRegion());
653  tmpDEM = m_TempDEMRegions[threadId];
654  }
655  else
656  {
657  return;
658  }
659 
660  itk::ImageRegionConstIteratorWithIndex<DisparityMapType> horizIt(horizDisp, disparityRegion);
661  itk::ImageRegionConstIteratorWithIndex<DisparityMapType> vertiIt(vertiDisp, disparityRegion);
662 
663  horizIt.GoToBegin();
664  vertiIt.GoToBegin();
665 
666  bool useMask = false;
667  itk::ImageRegionConstIterator<MaskImageType> maskIt;
668  if (disparityMask)
669  {
670  useMask = true;
671  maskIt = itk::ImageRegionConstIterator<MaskImageType>(disparityMask, disparityRegion);
672  maskIt.GoToBegin();
673  }
674 
675  typename TDisparityImage::PointType epiPoint;
676  itk::ContinuousIndex<double, 2> gridIndexConti;
677  double subPixIndex[2];
678  typename GridImageType::IndexType ulIndex, urIndex, lrIndex, llIndex;
679  typename GridImageType::PixelType ulPixel(2);
680  typename GridImageType::PixelType urPixel(2);
681  typename GridImageType::PixelType lrPixel(2);
682  typename GridImageType::PixelType llPixel(2);
683  typename GridImageType::PixelType cPixel(2);
684 
685  typename GridImageType::PointType ulPoint;
686  typename GridImageType::PointType urPoint;
687  typename GridImageType::PointType lrPoint;
688  typename GridImageType::PointType llPoint;
689 
690  TDPointType sensorPoint;
691  TDPointType leftGroundHmin;
692  TDPointType leftGroundHmax;
693  TDPointType rightGroundHmin;
694  TDPointType rightGroundHmax;
695 
696  while (!horizIt.IsAtEnd() && !vertiIt.IsAtEnd())
697  {
698  // check mask value if any
699  if (useMask)
700  {
701  if (!(maskIt.Get() > 0))
702  {
703  ++horizIt;
704  ++vertiIt;
705  ++maskIt;
706  continue;
707  }
708  }
709 
710  // compute left ray
711  horizDisp->TransformIndexToPhysicalPoint(horizIt.GetIndex(), epiPoint);
712  leftGrid->TransformPhysicalPointToContinuousIndex(epiPoint, gridIndexConti);
713 
714  ulIndex[0] = static_cast<int>(std::floor(gridIndexConti[0]));
715  ulIndex[1] = static_cast<int>(std::floor(gridIndexConti[1]));
716  if (ulIndex[0] < gridRegion.GetIndex(0))
717  ulIndex[0] = gridRegion.GetIndex(0);
718  if (ulIndex[1] < gridRegion.GetIndex(1))
719  ulIndex[1] = gridRegion.GetIndex(1);
720  if (ulIndex[0] > (gridRegion.GetIndex(0) + static_cast<int>(gridRegion.GetSize(0)) - 2))
721  {
722  ulIndex[0] = gridRegion.GetIndex(0) + gridRegion.GetSize(0) - 2;
723  }
724  if (ulIndex[1] > (gridRegion.GetIndex(1) + static_cast<int>(gridRegion.GetSize(1)) - 2))
725  {
726  ulIndex[1] = gridRegion.GetIndex(1) + gridRegion.GetSize(1) - 2;
727  }
728  urIndex[0] = ulIndex[0] + 1;
729  urIndex[1] = ulIndex[1];
730  lrIndex[0] = ulIndex[0] + 1;
731  lrIndex[1] = ulIndex[1] + 1;
732  llIndex[0] = ulIndex[0];
733  llIndex[1] = ulIndex[1] + 1;
734  subPixIndex[0] = gridIndexConti[0] - static_cast<double>(ulIndex[0]);
735  subPixIndex[1] = gridIndexConti[1] - static_cast<double>(ulIndex[1]);
736 
737  leftGrid->TransformIndexToPhysicalPoint(ulIndex, ulPoint);
738  leftGrid->TransformIndexToPhysicalPoint(urIndex, urPoint);
739  leftGrid->TransformIndexToPhysicalPoint(lrIndex, lrPoint);
740  leftGrid->TransformIndexToPhysicalPoint(llIndex, llPoint);
741 
742  ulPixel[0] = (leftGrid->GetPixel(ulIndex))[0] + ulPoint[0];
743  ulPixel[1] = (leftGrid->GetPixel(ulIndex))[1] + ulPoint[1];
744  urPixel[0] = (leftGrid->GetPixel(urIndex))[0] + urPoint[0];
745  urPixel[1] = (leftGrid->GetPixel(urIndex))[1] + urPoint[1];
746  lrPixel[0] = (leftGrid->GetPixel(lrIndex))[0] + lrPoint[0];
747  lrPixel[1] = (leftGrid->GetPixel(lrIndex))[1] + lrPoint[1];
748  llPixel[0] = (leftGrid->GetPixel(llIndex))[0] + llPoint[0];
749  llPixel[1] = (leftGrid->GetPixel(llIndex))[1] + llPoint[1];
750  cPixel = (ulPixel * (1.0 - subPixIndex[0]) + urPixel * subPixIndex[0]) * (1.0 - subPixIndex[1]) +
751  (llPixel * (1.0 - subPixIndex[0]) + lrPixel * subPixIndex[0]) * subPixIndex[1];
752 
753  sensorPoint[0] = cPixel[0];
754  sensorPoint[1] = cPixel[1];
755  sensorPoint[2] = m_ElevationMin;
756  leftGroundHmin = m_LeftToGroundTransform->TransformPoint(sensorPoint);
757 
758  sensorPoint[2] = m_ElevationMax;
759  leftGroundHmax = m_LeftToGroundTransform->TransformPoint(sensorPoint);
760 
761  // compute right ray
762  itk::ContinuousIndex<double, 2> rightIndexEstimate;
763  rightIndexEstimate[0] = static_cast<double>((horizIt.GetIndex())[0]) + static_cast<double>(horizIt.Get());
764  rightIndexEstimate[1] = static_cast<double>((horizIt.GetIndex())[1]) + static_cast<double>(vertiIt.Get());
765 
766  horizDisp->TransformContinuousIndexToPhysicalPoint(rightIndexEstimate, epiPoint);
767  rightGrid->TransformPhysicalPointToContinuousIndex(epiPoint, gridIndexConti);
768 
769  ulIndex[0] = static_cast<int>(std::floor(gridIndexConti[0]));
770  ulIndex[1] = static_cast<int>(std::floor(gridIndexConti[1]));
771  if (ulIndex[0] < gridRegion.GetIndex(0))
772  ulIndex[0] = gridRegion.GetIndex(0);
773  if (ulIndex[1] < gridRegion.GetIndex(1))
774  ulIndex[1] = gridRegion.GetIndex(1);
775  if (ulIndex[0] > (gridRegion.GetIndex(0) + static_cast<int>(gridRegion.GetSize(0)) - 2))
776  {
777  ulIndex[0] = gridRegion.GetIndex(0) + gridRegion.GetSize(0) - 2;
778  }
779  if (ulIndex[1] > (gridRegion.GetIndex(1) + static_cast<int>(gridRegion.GetSize(1)) - 2))
780  {
781  ulIndex[1] = gridRegion.GetIndex(1) + gridRegion.GetSize(1) - 2;
782  }
783  urIndex[0] = ulIndex[0] + 1;
784  urIndex[1] = ulIndex[1];
785  lrIndex[0] = ulIndex[0] + 1;
786  lrIndex[1] = ulIndex[1] + 1;
787  llIndex[0] = ulIndex[0];
788  llIndex[1] = ulIndex[1] + 1;
789  subPixIndex[0] = gridIndexConti[0] - static_cast<double>(ulIndex[0]);
790  subPixIndex[1] = gridIndexConti[1] - static_cast<double>(ulIndex[1]);
791 
792  rightGrid->TransformIndexToPhysicalPoint(ulIndex, ulPoint);
793  rightGrid->TransformIndexToPhysicalPoint(urIndex, urPoint);
794  rightGrid->TransformIndexToPhysicalPoint(lrIndex, lrPoint);
795  rightGrid->TransformIndexToPhysicalPoint(llIndex, llPoint);
796 
797  ulPixel[0] = (rightGrid->GetPixel(ulIndex))[0] + ulPoint[0];
798  ulPixel[1] = (rightGrid->GetPixel(ulIndex))[1] + ulPoint[1];
799  urPixel[0] = (rightGrid->GetPixel(urIndex))[0] + urPoint[0];
800  urPixel[1] = (rightGrid->GetPixel(urIndex))[1] + urPoint[1];
801  lrPixel[0] = (rightGrid->GetPixel(lrIndex))[0] + lrPoint[0];
802  lrPixel[1] = (rightGrid->GetPixel(lrIndex))[1] + lrPoint[1];
803  llPixel[0] = (rightGrid->GetPixel(llIndex))[0] + llPoint[0];
804  llPixel[1] = (rightGrid->GetPixel(llIndex))[1] + llPoint[1];
805  cPixel = (ulPixel * (1.0 - subPixIndex[0]) + urPixel * subPixIndex[0]) * (1.0 - subPixIndex[1]) +
806  (llPixel * (1.0 - subPixIndex[0]) + lrPixel * subPixIndex[0]) * subPixIndex[1];
807 
808  sensorPoint[0] = cPixel[0];
809  sensorPoint[1] = cPixel[1];
810  sensorPoint[2] = m_ElevationMin;
811  rightGroundHmin = m_RightToGroundTransform->TransformPoint(sensorPoint);
812 
813  sensorPoint[2] = m_ElevationMax;
814  rightGroundHmax = m_RightToGroundTransform->TransformPoint(sensorPoint);
815 
816  // Compute ray intersection (mid-point method), TODO : implement non-iterative method from Hartley & Sturm
817  double a = (leftGroundHmax[0] - leftGroundHmin[0]) * (leftGroundHmax[0] - leftGroundHmin[0]) +
818  (leftGroundHmax[1] - leftGroundHmin[1]) * (leftGroundHmax[1] - leftGroundHmin[1]) +
819  (leftGroundHmax[2] - leftGroundHmin[2]) * (leftGroundHmax[2] - leftGroundHmin[2]);
820  double b = (rightGroundHmax[0] - rightGroundHmin[0]) * (rightGroundHmax[0] - rightGroundHmin[0]) +
821  (rightGroundHmax[1] - rightGroundHmin[1]) * (rightGroundHmax[1] - rightGroundHmin[1]) +
822  (rightGroundHmax[2] - rightGroundHmin[2]) * (rightGroundHmax[2] - rightGroundHmin[2]);
823  double c = -(leftGroundHmax[0] - leftGroundHmin[0]) * (rightGroundHmax[0] - rightGroundHmin[0]) -
824  (leftGroundHmax[1] - leftGroundHmin[1]) * (rightGroundHmax[1] - rightGroundHmin[1]) -
825  (leftGroundHmax[2] - leftGroundHmin[2]) * (rightGroundHmax[2] - rightGroundHmin[2]);
826  double g = (leftGroundHmax[0] - leftGroundHmin[0]) * (rightGroundHmin[0] - leftGroundHmin[0]) +
827  (leftGroundHmax[1] - leftGroundHmin[1]) * (rightGroundHmin[1] - leftGroundHmin[1]) +
828  (leftGroundHmax[2] - leftGroundHmin[2]) * (rightGroundHmin[2] - leftGroundHmin[2]);
829  double h = -(rightGroundHmax[0] - rightGroundHmin[0]) * (rightGroundHmin[0] - leftGroundHmin[0]) -
830  (rightGroundHmax[1] - rightGroundHmin[1]) * (rightGroundHmin[1] - leftGroundHmin[1]) -
831  (rightGroundHmax[2] - rightGroundHmin[2]) * (rightGroundHmin[2] - leftGroundHmin[2]);
832 
833  double rLeft = (b * g - c * h) / (a * b - c * c);
834  double rRight = (a * h - c * g) / (a * b - c * c);
835 
836  TDPointType leftFoot;
837  leftFoot.SetToBarycentricCombination(leftGroundHmax, leftGroundHmin, rLeft);
838 
839  TDPointType rightFoot;
840  rightFoot.SetToBarycentricCombination(rightGroundHmax, rightGroundHmin, rRight);
841 
842  TDPointType midPoint3D;
843  midPoint3D.SetToMidPoint(leftFoot, rightFoot);
844 
845  // Is point inside DEM area ?
846  typename DEMImageType::PointType midPoint2D;
847  midPoint2D[0] = midPoint3D[0];
848  midPoint2D[1] = midPoint3D[1];
849  itk::ContinuousIndex<double, 2> midIndex;
850  outputDEM->TransformPhysicalPointToContinuousIndex(midPoint2D, midIndex);
851  typename DEMImageType::IndexType cellIndex;
852 
853  // TODO JGT check if cellIndex should be calculated from the center of the pixel
854  // TransformContinuousIndexToPhysicalPoint with index [0,0] returns Origin of image
855  // TransformContinuousIndexToPhysicalPoint with index [0.5,0.5] returns a slight difference from Origin of image
856  cellIndex[0] = static_cast<int>(std::floor(midIndex[0] + 0.5));
857  cellIndex[1] = static_cast<int>(std::floor(midIndex[1] + 0.5));
858 
859  if (outputRequestedRegion.IsInside(cellIndex))
860  {
861  // Estimate local reference elevation (average, DEM or geoid) => NO NEED, ALREADY HAVE 3D RAYS
862  // double localElevation = demHandler->GetHeightAboveEllipsoid(midPoint2D);
863 
864  // Add point to its corresponding cell (keep maximum)
865  DEMPixelType cellHeight = static_cast<DEMPixelType>(midPoint3D[2]);
866  if (cellHeight > tmpDEM->GetPixel(cellIndex) && cellHeight < static_cast<DEMPixelType>(m_ElevationMax))
867  {
868  tmpDEM->SetPixel(cellIndex, cellHeight);
869  }
870  }
871 
872  ++horizIt;
873  ++vertiIt;
874 
875  if (useMask)
876  ++maskIt;
877  }
878 }
879 
880 template <class TDisparityImage, class TInputImage, class TOutputDEMImage, class TEpipolarGridImage, class TMaskImage>
881 void DisparityMapToDEMFilter<TDisparityImage, TInputImage, TOutputDEMImage, TEpipolarGridImage, TMaskImage>::AfterThreadedGenerateData()
882 {
883  TOutputDEMImage* outputDEM = this->GetDEMOutput();
884 
885  if (m_TempDEMRegions.size() < 1)
886  {
887  outputDEM->FillBuffer(m_ElevationMin);
888  return;
889  }
890 
891  itk::ImageRegionIterator<DEMImageType> outputDEMIt(outputDEM, outputDEM->GetRequestedRegion());
892  itk::ImageRegionIterator<DEMImageType> firstDEMIt(m_TempDEMRegions[0], outputDEM->GetRequestedRegion());
893 
894  outputDEMIt.GoToBegin();
895  firstDEMIt.GoToBegin();
896 
897  // Copy first DEM
898  while (!outputDEMIt.IsAtEnd() && !firstDEMIt.IsAtEnd())
899  {
900  outputDEMIt.Set(firstDEMIt.Get());
901  ++outputDEMIt;
902  ++firstDEMIt;
903  }
904 
905  // Check DEMs from other threads and keep the maximum elevation
906  for (unsigned int i = 1; i < m_TempDEMRegions.size(); i++)
907  {
908  itk::ImageRegionIterator<DEMImageType> tmpDEMIt(m_TempDEMRegions[i], outputDEM->GetRequestedRegion());
909 
910  outputDEMIt.GoToBegin();
911  tmpDEMIt.GoToBegin();
912 
913  while (!outputDEMIt.IsAtEnd() && !tmpDEMIt.IsAtEnd())
914  {
915  if (tmpDEMIt.Get() > outputDEMIt.Get())
916  {
917  outputDEMIt.Set(tmpDEMIt.Get());
918  }
919  ++outputDEMIt;
920  ++tmpDEMIt;
921  }
922  }
923 }
924 }
925 
926 #endif
otb::DisparityMapToDEMFilter::GetDisparityMaskInput
const TMaskImage * GetDisparityMaskInput() const
Definition: otbDisparityMapToDEMFilter.hxx:175
otb::DisparityMapToDEMFilter::SetLeftEpipolarGridInput
void SetLeftEpipolarGridInput(const TEpipolarGridImage *grid)
Definition: otbDisparityMapToDEMFilter.hxx:88
otb::DisparityMapToDEMFilter::GetHorizontalDisparityMapInput
const TDisparityImage * GetHorizontalDisparityMapInput() const
Definition: otbDisparityMapToDEMFilter.hxx:112
otb::DisparityMapToDEMFilter::SetRightEpipolarGridInput
void SetRightEpipolarGridInput(const TEpipolarGridImage *grid)
Definition: otbDisparityMapToDEMFilter.hxx:96
otb::DisparityMapToDEMFilter::TDPointType
RSTransformType::InputPointType TDPointType
Definition: otbDisparityMapToDEMFilter.h:97
otb::DisparityMapToDEMFilter::GetLeftInput
const TInputImage * GetLeftInput() const
Definition: otbDisparityMapToDEMFilter.hxx:133
otb::DisparityMapToDEMFilter::SetLeftInput
void SetLeftInput(const TInputImage *image)
Definition: otbDisparityMapToDEMFilter.hxx:74
otb::DisparityMapToDEMFilter::SetVerticalDisparityMapInput
void SetVerticalDisparityMapInput(const TDisparityImage *vmap)
Definition: otbDisparityMapToDEMFilter.hxx:66
otb::DisparityMapToDEMFilter::ThreadedGenerateData
void ThreadedGenerateData(const RegionType &outputRegionForThread, itk::ThreadIdType threadId) override
Definition: otbDisparityMapToDEMFilter.hxx:631
otb::DisparityMapToDEMFilter::GetLeftEpipolarGridInput
const TEpipolarGridImage * GetLeftEpipolarGridInput() const
Definition: otbDisparityMapToDEMFilter.hxx:154
otb::DisparityMapToDEMFilter::GetRightEpipolarGridInput
const TEpipolarGridImage * GetRightEpipolarGridInput() const
Definition: otbDisparityMapToDEMFilter.hxx:165
otb::DisparityMapToDEMFilter::GetRightInput
const TInputImage * GetRightInput() const
Definition: otbDisparityMapToDEMFilter.hxx:143
otb::DisparityMapToDEMFilter::SetDisparityMaskInput
void SetDisparityMaskInput(const TMaskImage *mask)
Definition: otbDisparityMapToDEMFilter.hxx:104
otb::DisparityMapToDEMFilter::GetVerticalDisparityMapInput
const TDisparityImage * GetVerticalDisparityMapInput() const
Definition: otbDisparityMapToDEMFilter.hxx:123
otb::DisparityMapToDEMFilter::GetDEMOutput
const TOutputDEMImage * GetDEMOutput() const
Definition: otbDisparityMapToDEMFilter.hxx:185
otb::DisparityMapToDEMFilter::SetHorizontalDisparityMapInput
void SetHorizontalDisparityMapInput(const TDisparityImage *hmap)
Definition: otbDisparityMapToDEMFilter.hxx:58
otb::DisparityMapToDEMFilter::SetRightInput
void SetRightInput(const TInputImage *image)
Definition: otbDisparityMapToDEMFilter.hxx:81
otb::GenericRSTransform
This is the class to handle generic remote sensing transform.
Definition: otbGenericRSTransform.h:60
otb::GenericRSTransform::Pointer
itk::SmartPointer< Self > Pointer
Definition: otbGenericRSTransform.h:66
otb
The "otb" namespace contains all Orfeo Toolbox (OTB) classes.
Definition: otbJoinContainer.h:33
Generated at Fri Mar 21 2025 16:04:13 for OTB with   doxygen 1.9.1