OTB  9.0.0
Orfeo Toolbox
otbMulti3DMapToDEMFilter.hxx
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2005-2022 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 otbMulti3DMapToDEMFilter_hxx
22 #define otbMulti3DMapToDEMFilter_hxx
23 
24 #include "otbMulti3DMapToDEMFilter.h"
25 #include "itkImageRegionConstIteratorWithIndex.h"
26 #include "itkImageRegionIterator.h"
27 #include "otbStreamingStatisticsVectorImageFilter.h"
28 #include "otbNoDataHelper.h"
29 
30 namespace otb
31 {
32 
33 template <class T3DImage, class TMaskImage, class TOutputDEMImage>
34 Multi3DMapToDEMFilter<T3DImage, TMaskImage, TOutputDEMImage>::Multi3DMapToDEMFilter()
35 {
36  // Set the number of inputs (1 image one optional mask)
37  this->SetNumberOfIndexedInputs(2);
38  this->SetNumberOfRequiredInputs(1);
39  // this->m_MapKeywordLists.resize(1);
40  m_DEMGridStep = 10.0;
41  // Set the outputs
42  this->SetNumberOfIndexedOutputs(1);
43  this->SetNthOutput(0, TOutputDEMImage::New());
44  // Default DEM reconstruction parameters
45  m_MapSplitterList = SplitterListType::New();
46 
47  m_NoDataValue = -32768;
48  m_ElevationMin = -100;
49  m_ElevationMax = 500;
50  m_CellFusionMode = otb::CellFusionMode::MAX;
51  m_OutputParametersFrom3DMap = -2;
52  m_IsGeographic = true;
53 
54  m_Margin[0] = 10;
55  m_Margin[1] = 10;
56 }
57 
58 template <class T3DImage, class TMaskImage, class TOutputDEMImage>
59 Multi3DMapToDEMFilter<T3DImage, TMaskImage, TOutputDEMImage>::~Multi3DMapToDEMFilter()
60 {
61 }
62 
63 template <class T3DImage, class TMaskImage, class TOutputDEMImage>
64 void Multi3DMapToDEMFilter<T3DImage, TMaskImage, TOutputDEMImage>::SetNumberOf3DMaps(unsigned int nb)
65 {
66  if (nb > 0)
67  {
68  this->SetNumberOfIndexedInputs(2 * nb);
69  }
70 }
71 
72 template <class T3DImage, class TMaskImage, class TOutputDEMImage>
73 unsigned int Multi3DMapToDEMFilter<T3DImage, TMaskImage, TOutputDEMImage>::GetNumberOf3DMaps()
74 {
75  return this->GetNumberOfInputs() / 2;
76 }
77 
78 template <class T3DImage, class TMaskImage, class TOutputDEMImage>
79 void Multi3DMapToDEMFilter<T3DImage, TMaskImage, TOutputDEMImage>::Set3DMapInput(unsigned int index, const T3DImage* map)
80 {
81  if ((2 * (index + 1)) > this->GetNumberOfInputs())
82  {
83  itkExceptionMacro(<< "Index is greater than the number of images");
84  }
85  // Process object is not const-correct so the const casting is required.
86  this->SetNthInput(2 * index, const_cast<T3DImage*>(map));
87 }
88 
89 
90 template <class T3DImage, class TMaskImage, class TOutputDEMImage>
91 void Multi3DMapToDEMFilter<T3DImage, TMaskImage, TOutputDEMImage>::SetMaskInput(unsigned int index, const TMaskImage* mask)
92 {
93  if ((2 * (index + 1)) > this->GetNumberOfInputs())
94  {
95  itkExceptionMacro(<< "Index is greater than the number of images");
96  }
97  // Process object is not const-correct so the const casting is required.
98  this->SetNthInput(2 * index + 1, const_cast<TMaskImage*>(mask));
99 }
100 
101 template <class T3DImage, class TMaskImage, class TOutputDEMImage>
102 const T3DImage* Multi3DMapToDEMFilter<T3DImage, TMaskImage, TOutputDEMImage>::Get3DMapInput(unsigned int index) const
103 {
104  if ((2 * (index + 1)) > this->GetNumberOfInputs())
105  {
106  return nullptr;
107  }
108  return static_cast<const T3DImage*>(this->itk::ProcessObject::GetInput(2 * index));
109 }
110 
111 template <class T3DImage, class TMaskImage, class TOutputDEMImage>
112 const TMaskImage* Multi3DMapToDEMFilter<T3DImage, TMaskImage, TOutputDEMImage>::GetMaskInput(unsigned int index) const
113 {
114  if ((2 * (index + 1)) > this->GetNumberOfInputs())
115  {
116  return nullptr;
117  }
118  return static_cast<const TMaskImage*>(this->itk::ProcessObject::GetInput(2 * index + 1));
119 }
120 
121 template <class T3DImage, class TMaskImage, class TOutputDEMImage>
122 const TOutputDEMImage* Multi3DMapToDEMFilter<T3DImage, TMaskImage, TOutputDEMImage>::GetDEMOutput() const
123 {
124  if (this->GetNumberOfOutputs() < 1)
125  {
126  return 0;
127  }
128  return static_cast<const TOutputDEMImage*>(this->itk::ProcessObject::GetOutput(0));
129 }
130 
131 template <class T3DImage, class TMaskImage, class TOutputDEMImage>
132 TOutputDEMImage* Multi3DMapToDEMFilter<T3DImage, TMaskImage, TOutputDEMImage>::GetDEMOutput()
133 {
134  if (this->GetNumberOfOutputs() < 1)
135  {
136  return nullptr;
137  }
138  return static_cast<TOutputDEMImage*>(this->itk::ProcessObject::GetOutput(0));
139 }
140 
141 
142 template <class T3DImage, class TMaskImage, class TOutputDEMImage>
143 void Multi3DMapToDEMFilter<T3DImage, TMaskImage, TOutputDEMImage>::SetOutputParametersFromImage()
144 {
145  int index = m_OutputParametersFrom3DMap;
146  if (static_cast<unsigned int>((2 * (index + 1))) > this->GetNumberOfInputs())
147  {
148  itkExceptionMacro(<< "input at position " << index << " is unavailable");
149  }
150 
151  unsigned int numberOfInputs = this->GetNumberOfInputs() / 2;
152  unsigned int indexStart = 0;
153  unsigned int indexEnd = numberOfInputs - 1;
154 
155  if (index != -1)
156  {
157  indexStart = index;
158  indexEnd = index;
159  }
160 
161  // compute DEM extent union of 3D map extent
162 
163  TOutputDEMImage* outputPtr = this->GetDEMOutput();
164 
165  // Set-up a transform to use the DEMHandler
166 
167  // DEM BBox
168  itk::NumericTraits<DEMPixelType>::max();
169  double box_xmin = itk::NumericTraits<DEMPixelType>::max();
170  double box_xmax = itk::NumericTraits<DEMPixelType>::NonpositiveMin();
171  double box_ymin = itk::NumericTraits<DEMPixelType>::max();
172  double box_ymax = itk::NumericTraits<DEMPixelType>::NonpositiveMin();
173 
174  for (unsigned int k = indexStart; k <= indexEnd; ++k)
175  {
176  T3DImage* imgPtr = const_cast<T3DImage*>(this->Get3DMapInput(k));
177 
178  RSTransform2DType::Pointer mapToGroundTransform = RSTransform2DType::New();
179  mapToGroundTransform->SetInputImageMetadata(&(imgPtr->GetImageMetadata()));
180 
181  /*if(!m_ProjectionRef.empty())
182  {
183  mapToGroundTransform->SetOutputProjectionRef(m_ProjectionRef);
184  }*/
185  mapToGroundTransform->InstantiateTransform();
186 
187  typename InputMapType::SizeType inputSize = imgPtr->GetLargestPossibleRegion().GetSize();
188 
189  typename InputMapType::PointType ulPt, urPt, llPt, lrPt;
190  itk::ContinuousIndex<double, 2> ulIdx(imgPtr->GetLargestPossibleRegion().GetIndex());
191  ulIdx[0] += -0.5;
192  ulIdx[1] += -0.5;
193 
194  itk::ContinuousIndex<double, 2> urIdx(ulIdx);
195  itk::ContinuousIndex<double, 2> lrIdx(ulIdx);
196  itk::ContinuousIndex<double, 2> llIdx(ulIdx);
197  urIdx[0] += static_cast<double>(inputSize[0]);
198  lrIdx[0] += static_cast<double>(inputSize[0]);
199  lrIdx[1] += static_cast<double>(inputSize[1]);
200  llIdx[1] += static_cast<double>(inputSize[1]);
201 
202  imgPtr->TransformContinuousIndexToPhysicalPoint(ulIdx, ulPt);
203  imgPtr->TransformContinuousIndexToPhysicalPoint(urIdx, urPt);
204  imgPtr->TransformContinuousIndexToPhysicalPoint(llIdx, llPt);
205  imgPtr->TransformContinuousIndexToPhysicalPoint(lrIdx, lrPt);
206 
207  RSTransform2DType::OutputPointType ul, ur, lr, ll;
208  ul = mapToGroundTransform->TransformPoint(ulPt);
209  ur = mapToGroundTransform->TransformPoint(urPt);
210  ll = mapToGroundTransform->TransformPoint(llPt);
211  lr = mapToGroundTransform->TransformPoint(lrPt);
212 
213  double xmin = std::min(std::min(std::min(ul[0], ur[0]), lr[0]), ll[0]);
214  double xmax = std::max(std::max(std::max(ul[0], ur[0]), lr[0]), ll[0]);
215  double ymin = std::min(std::min(std::min(ul[1], ur[1]), lr[1]), ll[1]);
216  double ymax = std::max(std::max(std::max(ul[1], ur[1]), lr[1]), ll[1]);
217 
218  box_xmin = std::min(box_xmin, xmin);
219  box_xmax = std::max(box_xmax, xmax);
220  box_ymin = std::min(box_ymin, ymin);
221  box_ymax = std::max(box_ymax, ymax);
222  }
223 
224  // Compute step :
225  // TODO : use a clean RS transform instead
226  typename TOutputDEMImage::SpacingType outSpacing;
227  // std::cout<<" GrisStep "<<m_DEMGridStep<<std::endl;
228  outSpacing[0] = 57.295779513 * m_DEMGridStep / (6378137.0 * std::cos((box_ymin + box_ymax) * 0.5 * 0.01745329251));
229  outSpacing[1] = -57.295779513 * m_DEMGridStep / 6378137.0;
230  outputPtr->SetSignedSpacing(outSpacing);
231 
232  // Choose origin
233  typename TOutputDEMImage::PointType outOrigin;
234  outOrigin[0] = box_xmin + 0.5 * outSpacing[0];
235  outOrigin[1] = box_ymax + 0.5 * outSpacing[1];
236  outputPtr->SetOrigin(outOrigin);
237 
238  // Compute output size
239  typename TOutputDEMImage::RegionType outRegion;
240  outRegion.SetIndex(0, 0);
241  outRegion.SetIndex(1, 0);
242  outRegion.SetSize(0, static_cast<unsigned int>(std::floor((box_xmax - box_xmin) / std::abs(outSpacing[0]) + 0.5)));
243  outRegion.SetSize(1, static_cast<unsigned int>(std::floor((box_ymax - box_ymin) / std::abs(outSpacing[1]) + 0.5)));
244  outputPtr->SetLargestPossibleRegion(outRegion);
245  outputPtr->SetNumberOfComponentsPerPixel(1);
246 
247 
248  itk::MetaDataDictionary& dictOutput = outputPtr->GetMetaDataDictionary();
249  itk::EncapsulateMetaData<std::string>(dictOutput, MetaDataKey::ProjectionRefKey, static_cast<std::string>(otb::SpatialReference::FromWGS84().ToWkt()));
250 
251  // test if WGS 84 -> true -> nothing to do
252 
253  // false project
254 
255  bool isWGS84 = !(m_ProjectionRef.compare(static_cast<std::string>(otb::SpatialReference::FromWGS84().ToWkt())));
256  if (!m_ProjectionRef.empty() && !isWGS84)
257  {
258 
259  typename OutputParametersEstimatorType::Pointer genericRSEstimator = OutputParametersEstimatorType::New();
260 
261  genericRSEstimator->SetInput(outputPtr);
262  // genericRSEstimator->SetInputProjectionRef( static_cast<std::string>(otb::SpatialReference().ToWkt()));
263  genericRSEstimator->SetOutputProjectionRef(m_ProjectionRef);
264  genericRSEstimator->Compute();
265  outputPtr->SetSignedSpacing(genericRSEstimator->GetOutputSpacing());
266  outputPtr->SetOrigin(genericRSEstimator->GetOutputOrigin());
267 
268  // Compute output size
269  typename TOutputDEMImage::RegionType outRegion2;
270  outRegion2.SetIndex(0, 0);
271  outRegion2.SetIndex(1, 0);
272  outRegion2.SetSize(0, genericRSEstimator->GetOutputSize()[0]);
273  // TODO JGT check the size
274  // outRegion.SetSize(1, static_cast<unsigned int> ((box_ymax - box_ymin) / std::abs(outSpacing[1])+1));
275  outRegion2.SetSize(1, genericRSEstimator->GetOutputSize()[1]);
276  outputPtr->SetLargestPossibleRegion(outRegion2);
277  outputPtr->SetNumberOfComponentsPerPixel(1);
278 
279 
280  itk::MetaDataDictionary& dict = outputPtr->GetMetaDataDictionary();
281 
282  itk::EncapsulateMetaData<std::string>(dict, MetaDataKey::ProjectionRefKey, m_ProjectionRef);
283  }
284  this->Modified();
285 }
286 
287 template <class T3DImage, class TMaskImage, class TOutputDEMImage>
288 void Multi3DMapToDEMFilter<T3DImage, TMaskImage, TOutputDEMImage>::GenerateOutputInformation()
289 {
290  //
291  TOutputDEMImage* outputPtr = this->GetDEMOutput();
292  //
293  if (this->m_OutputParametersFrom3DMap == -2)
294  {
295  outputPtr->SetOrigin(m_OutputOrigin);
296  outputPtr->SetSignedSpacing(m_OutputSpacing);
297 
298  typename TOutputDEMImage::RegionType outRegion;
299  outRegion.SetIndex(m_OutputStartIndex);
300  outRegion.SetSize(m_OutputSize);
301  outputPtr->SetLargestPossibleRegion(outRegion);
302  outputPtr->SetNumberOfComponentsPerPixel(1);
303 
304  if (!m_ProjectionRef.empty())
305  {
306  // fill up the metadata information for ProjectionRef
307  itk::MetaDataDictionary& dict = this->GetOutput()->GetMetaDataDictionary();
308  itk::EncapsulateMetaData<std::string>(dict, MetaDataKey::ProjectionRefKey, m_ProjectionRef);
309  }
310  }
311  else
312  {
313  this->SetOutputParametersFromImage();
314  }
315 
316  if (!m_ProjectionRef.empty())
317  {
318  OGRSpatialReference oSRS;
319 #if GDAL_VERSION_NUM >= 3000000 // importFromWkt is const-correct in GDAL 3
320  const char* wkt = m_ProjectionRef.c_str();
321 #else
322  char* wkt = const_cast<char*>(m_ProjectionRef.c_str());
323 #endif
324  oSRS.importFromWkt(&wkt);
325  m_IsGeographic = oSRS.IsGeographic(); // TODO check if this test is valid for all projection systems
326  }
327 
328  // Set the NoData value
329  std::vector<bool> noDataValueAvailable;
330  noDataValueAvailable.push_back(true);
331  std::vector<double> noDataValue;
332  noDataValue.push_back(m_NoDataValue);
333 
334  WriteNoDataFlags(noDataValueAvailable, noDataValue, outputPtr->GetImageMetadata());
335 }
336 
337 template <class T3DImage, class TMaskImage, class TOutputDEMImage>
338 void Multi3DMapToDEMFilter<T3DImage, TMaskImage, TOutputDEMImage>::GenerateInputRequestedRegion()
339 {
340  const TOutputDEMImage* outputDEM = this->GetDEMOutput();
341 
342  typename TOutputDEMImage::RegionType outRegion = outputDEM->GetRequestedRegion();
343  typename TOutputDEMImage::PointType outOrigin = outputDEM->GetOrigin();
344  typename TOutputDEMImage::SpacingType outSpacing = outputDEM->GetSignedSpacing();
345 
346  // up left at elevation min
347  TDPointType corners[8];
348  corners[0][0] = outOrigin[0] + outSpacing[0] * (-0.5 + static_cast<double>(outRegion.GetIndex(0)));
349  corners[0][1] = outOrigin[1] + outSpacing[1] * (-0.5 + static_cast<double>(outRegion.GetIndex(1)));
350  corners[0][2] = m_ElevationMin;
351  // up left at elevation max
352  corners[1][0] = corners[0][0];
353  corners[1][1] = corners[0][1];
354  corners[1][2] = m_ElevationMax;
355  // up right at elevation min
356  corners[2][0] = corners[0][0] + outSpacing[0] * static_cast<double>(outRegion.GetSize(0));
357  corners[2][1] = corners[0][1];
358  corners[2][2] = m_ElevationMin;
359  // up right at elevation max
360  corners[3][0] = corners[2][0];
361  corners[3][1] = corners[2][1];
362  corners[3][2] = m_ElevationMax;
363  // low right at elevation min
364  corners[4][0] = corners[0][0] + outSpacing[0] * static_cast<double>(outRegion.GetSize(0));
365  corners[4][1] = corners[0][1] + outSpacing[1] * static_cast<double>(outRegion.GetSize(1));
366  corners[4][2] = m_ElevationMin;
367  // low right at elevation max
368  corners[5][0] = corners[4][0];
369  corners[5][1] = corners[4][1];
370  corners[5][2] = m_ElevationMax;
371  // low left at elevation min
372  corners[6][0] = corners[0][0];
373  corners[6][1] = corners[0][1] + outSpacing[1] * static_cast<double>(outRegion.GetSize(1));
374  corners[6][2] = m_ElevationMin;
375  // low left at elevation max
376  corners[7][0] = corners[6][0];
377  corners[7][1] = corners[6][1];
378  corners[7][2] = m_ElevationMax;
379 
380  for (unsigned int k = 0; k < this->GetNumberOf3DMaps(); ++k)
381  {
382 
383  // set requested to largest and check that mask has the same size
384  T3DImage* imgPtr = const_cast<T3DImage*>(this->Get3DMapInput(k));
385 
386  RSTransformType::Pointer groundToSensorTransform = RSTransformType::New();
387  // groundToSensorTransform->SetInputOrigin(outputDEM->GetOrigin());
388  // groundToSensorTransform->SetInputSpacing(outputDEM->GetSignedSpacing());
389  groundToSensorTransform->SetInputProjectionRef(m_ProjectionRef);
390 
391  groundToSensorTransform->SetOutputImageMetadata(&(imgPtr->GetImageMetadata()));
392  groundToSensorTransform->SetOutputOrigin(imgPtr->GetOrigin());
393  groundToSensorTransform->SetOutputSpacing(imgPtr->GetSignedSpacing());
394  groundToSensorTransform->InstantiateTransform();
395 
396  typename T3DImage::RegionType mapRegion = imgPtr->GetLargestPossibleRegion();
397 
398  itk::ContinuousIndex<double, 2> mapContiIndex;
399  long int maxMapIndex[2] = {0, 0};
400  long int minMapIndex[2] = {0, 0};
401  maxMapIndex[0] = static_cast<long int>(mapRegion.GetIndex(0) + mapRegion.GetSize(0));
402  maxMapIndex[1] = static_cast<long int>(mapRegion.GetIndex(1) + mapRegion.GetSize(1));
403  minMapIndex[0] = static_cast<long int>(mapRegion.GetIndex(0));
404  minMapIndex[1] = static_cast<long int>(mapRegion.GetIndex(1));
405 
406  long int minMapRequestedIndex[2] = {0, 0};
407  minMapRequestedIndex[0] = maxMapIndex[0] + 1;
408  minMapRequestedIndex[1] = maxMapIndex[1] + 1;
409 
410  long int maxMapRequestedIndex[2] = {0, 0};
411  maxMapRequestedIndex[0] = 0;
412  maxMapRequestedIndex[1] = 0;
413 
414  for (unsigned int i = 0; i < 8; i++)
415  {
416  TDPointType tmpSensor = groundToSensorTransform->TransformPoint(corners[i]);
417 
418  minMapRequestedIndex[0] = std::min(minMapRequestedIndex[0], static_cast<long int>(tmpSensor[0] - m_Margin[0]));
419 
420  minMapRequestedIndex[1] = std::min(minMapRequestedIndex[1], static_cast<long int>(tmpSensor[1] - m_Margin[1]));
421 
422  maxMapRequestedIndex[0] = std::max(maxMapRequestedIndex[0], static_cast<long int>(tmpSensor[0] + m_Margin[0]));
423 
424  maxMapRequestedIndex[1] = std::max(maxMapRequestedIndex[1], static_cast<long int>(tmpSensor[1] + m_Margin[1]));
425 
426  minMapRequestedIndex[0] = std::max(minMapRequestedIndex[0], minMapIndex[0]);
427  minMapRequestedIndex[1] = std::max(minMapRequestedIndex[1], minMapIndex[1]);
428  maxMapRequestedIndex[0] = std::min(maxMapRequestedIndex[0], maxMapIndex[0]);
429  maxMapRequestedIndex[1] = std::min(maxMapRequestedIndex[1], maxMapIndex[1]);
430  }
431 
432  RegionType largest = imgPtr->GetLargestPossibleRegion();
433  RegionType requestedRegion = largest;
434 
435  if ((minMapRequestedIndex[0] < maxMapRequestedIndex[0]) && (minMapRequestedIndex[1] < maxMapRequestedIndex[1]))
436  {
437  requestedRegion.SetSize(0, maxMapRequestedIndex[0] - minMapRequestedIndex[0]);
438  requestedRegion.SetSize(1, maxMapRequestedIndex[1] - minMapRequestedIndex[1]);
439  requestedRegion.SetIndex(0, minMapRequestedIndex[0]);
440  requestedRegion.SetIndex(1, minMapRequestedIndex[1]);
441  }
442  else
443  {
444  requestedRegion.SetSize(0, 0);
445  requestedRegion.SetSize(1, 0);
446  requestedRegion.SetIndex(0, minMapIndex[0]);
447  requestedRegion.SetIndex(1, minMapIndex[1]);
448  }
449 
450  imgPtr->SetRequestedRegion(requestedRegion);
451 
452  TMaskImage* mskPtr = const_cast<TMaskImage*>(this->GetMaskInput(k));
453  if (mskPtr)
454  {
455 
456  if (mskPtr->GetLargestPossibleRegion() != largest)
457  {
458  itkExceptionMacro(<< "mask and map at position " << k << " have a different largest region");
459  }
460  mskPtr->SetRequestedRegion(requestedRegion);
461  }
462  }
463 }
464 
465 template <class T3DImage, class TMaskImage, class TOutputDEMImage>
466 void Multi3DMapToDEMFilter<T3DImage, TMaskImage, TOutputDEMImage>::BeforeThreadedGenerateData()
467 {
468  const TOutputDEMImage* outputDEM = this->GetDEMOutput();
469 
470  // create splits
471  // for each map we check if the input region can be split into threadNb
472  m_NumberOfSplit.resize(this->GetNumberOf3DMaps());
473 
474  unsigned int maximumRegionsNumber = 1;
475 
476  for (unsigned int k = 0; k < this->GetNumberOf3DMaps(); ++k)
477  {
478  m_MapSplitterList->PushBack(SplitterType::New());
479  T3DImage* imgPtr = const_cast<T3DImage*>(this->Get3DMapInput(k));
480 
481  typename T3DImage::RegionType requestedRegion = imgPtr->GetRequestedRegion();
482 
483  typename T3DImage::SizeType requestedSize = requestedRegion.GetSize();
484  unsigned int regionsNumber = 0;
485  if (requestedSize[0] * requestedSize[1] != 0)
486  {
487  regionsNumber = m_MapSplitterList->GetNthElement(k)->GetNumberOfSplits(requestedRegion, this->GetNumberOfThreads());
488  }
489  m_NumberOfSplit[k] = regionsNumber;
490  otbMsgDevMacro("map " << k << " will be split into " << regionsNumber << " regions");
491  if (maximumRegionsNumber < regionsNumber)
492  maximumRegionsNumber = regionsNumber;
493  }
494 
495  m_TempDEMRegions.clear();
496  m_TempDEMAccumulatorRegions.clear();
497  // m_ThreadProcessed.resize(maximumRegionsNumber);
498 
499  for (unsigned int i = 0; i < maximumRegionsNumber; i++)
500  {
501  // m_ThreadProcessed[i] = 0;
502  typename TOutputDEMImage::Pointer tmpImg = TOutputDEMImage::New();
503  tmpImg->SetNumberOfComponentsPerPixel(1); // Two components for mean calculus ?
504  tmpImg->SetRegions(outputDEM->GetRequestedRegion());
505  tmpImg->Allocate();
506 
507  tmpImg->FillBuffer(m_NoDataValue);
508  m_TempDEMRegions.push_back(tmpImg);
509 
510  typename AccumulatorImageType::Pointer tmpImg2 = AccumulatorImageType::New();
511  tmpImg2->SetNumberOfComponentsPerPixel(1); // Two components for mean calculus ?
512  tmpImg2->SetRegions(outputDEM->GetRequestedRegion());
513  tmpImg2->Allocate();
514 
515  tmpImg2->FillBuffer(0.);
516  m_TempDEMAccumulatorRegions.push_back(tmpImg2);
517  }
518 
519  if (!this->m_IsGeographic)
520  {
521  m_GroundTransform = RSTransform2DType::New();
522  m_GroundTransform->SetInputProjectionRef(static_cast<std::string>(otb::SpatialReference::FromWGS84().ToWkt()));
523  m_GroundTransform->SetOutputProjectionRef(m_ProjectionRef);
524  m_GroundTransform->InstantiateTransform();
525  }
526 }
527 
528 template <class T3DImage, class TMaskImage, class TOutputDEMImage>
529 void Multi3DMapToDEMFilter<T3DImage, TMaskImage, TOutputDEMImage>::ThreadedGenerateData(const RegionType& itkNotUsed(outputRegionForThread),
530  itk::ThreadIdType threadId)
531 {
532  TOutputDEMImage* outputPtr = this->GetOutput();
533 
534  typename OutputImageType::PointType pointRef;
535  typename OutputImageType::PointType pointRefStep;
536  typename OutputImageType::RegionType requestedRegion = outputPtr->GetRequestedRegion();
537 
538  // typename TOutputDEMImage::SpacingType step = outputPtr->GetSignedSpacing();
539 
540  // convert requested region to Long/Lat
541 
542  // typename TOutputDEMImage::SizeType size = requestedRegion.GetSize();
543 
544  typename TOutputDEMImage::IndexType index = requestedRegion.GetIndex();
545  outputPtr->TransformIndexToPhysicalPoint(index, pointRef);
546  /*
547  InputInternalPixelType regionLong1 = pointRef[0];
548  InputInternalPixelType regionLat1 = pointRef[1];
549  InputInternalPixelType regionLong2 = pointRef[0] + size[0] * step[0];
550  InputInternalPixelType regionLat2 = pointRef[1] + size[1] * step[1];
551  InputInternalPixelType minLong = std::min(regionLong1, regionLong2);
552  InputInternalPixelType minLat = std::min(regionLat1, regionLat2);
553  InputInternalPixelType maxLong = std::max(regionLong1, regionLong2);
554  InputInternalPixelType maxLat = std::max(regionLat1, regionLat2);
555  */
556 
557  TOutputDEMImage* tmpDEM = nullptr;
558  AccumulatorImageType* tmpAcc = nullptr;
559  typename TOutputDEMImage::RegionType outputRequestedRegion = outputPtr->GetRequestedRegion();
560 
561  typename T3DImage::RegionType splitRegion;
562 
563  MapPixelType position;
564 
565  itk::ImageRegionConstIterator<InputMapType> mapIt;
566  for (unsigned int k = 0; k < this->GetNumberOf3DMaps(); ++k)
567  {
568  if (m_NumberOfSplit[k] > 0)
569  {
570  T3DImage* imgPtr = const_cast<T3DImage*>(this->Get3DMapInput(k));
571  TMaskImage* mskPtr = const_cast<TMaskImage*>(this->GetMaskInput(k));
572 
573  typename InputMapType::PointType origin;
574  origin = imgPtr->GetOrigin();
575  typename InputMapType::SpacingType spacing;
576  spacing = imgPtr->GetSignedSpacing();
577 
578  if (static_cast<unsigned int>(threadId) < m_NumberOfSplit[k])
579  {
580  splitRegion = m_MapSplitterList->GetNthElement(k)->GetSplit(threadId, m_NumberOfSplit[k], imgPtr->GetRequestedRegion());
581 
582  tmpDEM = m_TempDEMRegions[threadId];
583  tmpAcc = m_TempDEMAccumulatorRegions[threadId];
584 
585  mapIt = itk::ImageRegionConstIterator<InputMapType>(imgPtr, splitRegion);
586  mapIt.GoToBegin();
587  itk::ImageRegionConstIterator<MaskImageType> maskIt;
588  bool useMask = false;
589  if (mskPtr)
590  {
591  useMask = true;
592  maskIt = itk::ImageRegionConstIterator<MaskImageType>(mskPtr, splitRegion);
593  maskIt.GoToBegin();
594  }
595 
596  while (!mapIt.IsAtEnd())
597  {
598  // check mask value if any
599  if (useMask)
600  {
601  if (!(maskIt.Get() > 0))
602  {
603  ++mapIt;
604  ++maskIt;
605  continue;
606  }
607  }
608 
609  position = mapIt.Get();
610 
611  // std::cout<<"position"<<position<<std::endl;
612  if (!this->m_IsGeographic)
613  {
614  // std::cout<<"is geographic "<<std::endl;
615  typename RSTransform2DType::InputPointType tmpPoint;
616  tmpPoint[0] = position[0];
617  tmpPoint[1] = position[1];
618  RSTransform2DType::OutputPointType groundPosition = m_GroundTransform->TransformPoint(tmpPoint);
619  position[0] = groundPosition[0];
620  position[1] = groundPosition[1];
621  }
622 
623  // test if position is in DEM BBOX
628  // Is point inside DEM area ?
629  typename OutputImageType::PointType point2D;
630  point2D[0] = position[0];
631  point2D[1] = position[1];
632  itk::ContinuousIndex<double, 2> continuousIndex;
634 
635  // std::cout<<"point2D "<<point2D<<std::endl;
636  // The DEM cell at index 'n' contains continuous indexes from 'n-0.5' to 'n+0.5'
637  outputPtr->TransformPhysicalPointToContinuousIndex(point2D, continuousIndex);
638  typename OutputImageType::IndexType cellIndex;
639  cellIndex[0] = static_cast<int>(std::floor(continuousIndex[0] + 0.5));
640  cellIndex[1] = static_cast<int>(std::floor(continuousIndex[1] + 0.5));
641  // std::cout<<"cellindex "<<cellIndex<<std::endl;
642  // index from physical
650  if (outputRequestedRegion.IsInside(cellIndex))
651  {
652  // std::cout<<"is inside "<<std::endl;
653  // Add point to its corresponding cell (keep maximum)
654  DEMPixelType cellHeight = static_cast<DEMPixelType>(position[2]);
655  // if (cellHeight > tmpDEM->GetPixel(cellIndex) && cellHeight < static_cast<DEMPixelType>(m_ElevationMax))
656  // {
657  // tmpDEM->SetPixel(cellIndex,tmpDEM->GetPixel(cellIndex)+1);
658 
659  AccumulatorPixelType accPixel = tmpAcc->GetPixel(cellIndex);
660  tmpAcc->SetPixel(cellIndex, tmpAcc->GetPixel(cellIndex) + 1);
661 
662  if (accPixel == 0)
663  {
664  tmpDEM->SetPixel(cellIndex, cellHeight);
665  }
666  else
667  {
668  DEMPixelType cellCurrentValue = tmpDEM->GetPixel(cellIndex);
669 
670  switch (this->m_CellFusionMode)
671  {
672  case otb::CellFusionMode::MIN:
673  {
674  if (cellHeight < cellCurrentValue)
675  {
676  tmpDEM->SetPixel(cellIndex, cellHeight);
677  }
678  }
679  break;
680  case otb::CellFusionMode::MAX:
681  {
682  if (cellHeight > cellCurrentValue)
683  {
684  tmpDEM->SetPixel(cellIndex, cellHeight);
685  }
686  }
687  break;
688  case otb::CellFusionMode::MEAN:
689  {
690  tmpDEM->SetPixel(cellIndex, cellCurrentValue + cellHeight);
691  }
692  break;
693  case otb::CellFusionMode::ACC:
694  {
695  }
696  break;
697  default:
698 
699  itkExceptionMacro(<< "Unexpected value cell fusion mode :" << this->m_CellFusionMode);
700  break;
701  }
702  }
703  }
704 
705  ++mapIt;
706 
707  if (useMask)
708  ++maskIt;
709  }
710  }
711  else
712  {
713  splitRegion = requestedRegion;
714  otbMsgDevMacro("map " << k << " will not be split ");
715  }
716  }
717  }
718 }
719 
720 template <class T3DImage, class TMaskImage, class TOutputDEMImage>
721 void Multi3DMapToDEMFilter<T3DImage, TMaskImage, TOutputDEMImage>::AfterThreadedGenerateData()
722 {
723 
724  TOutputDEMImage* outputDEM = this->GetOutput();
725 
726  // check is that case can occur
727  if (m_TempDEMRegions.size() < 1)
728  {
729  outputDEM->FillBuffer(m_NoDataValue);
730  return;
731  }
732 
733  itk::ImageRegionIterator<OutputImageType> outputDEMIt(outputDEM, outputDEM->GetRequestedRegion());
734  itk::ImageRegionIterator<OutputImageType> firstDEMIt(m_TempDEMRegions[0], outputDEM->GetRequestedRegion());
735  itk::ImageRegionIterator<AccumulatorImageType> firstDEMAccIt(m_TempDEMAccumulatorRegions[0], outputDEM->GetRequestedRegion());
736  // we use the first accumulator as global accumulator over tmpAcc for mean calculus
737 
738  outputDEMIt.GoToBegin();
739  firstDEMIt.GoToBegin();
740  firstDEMAccIt.GoToBegin();
741  // Copy first DEM
742 
743  while (!outputDEMIt.IsAtEnd() && !firstDEMIt.IsAtEnd() && !firstDEMAccIt.IsAtEnd())
744  {
745 
746  AccumulatorPixelType accPixel = firstDEMAccIt.Get();
747  // useless test tempDEm is initialized with NoDataValue
748  if (accPixel == 0)
749  {
750  outputDEMIt.Set(m_NoDataValue);
751  }
752  else
753  {
754 
755  DEMPixelType pixelValue = firstDEMIt.Get();
756 
757  outputDEMIt.Set(pixelValue);
758 
759  if ((this->m_CellFusionMode == otb::CellFusionMode::MEAN) && (m_TempDEMRegions.size() == 1))
760  {
761  outputDEMIt.Set(firstDEMIt.Get() / static_cast<DEMPixelType>(accPixel));
762  }
763  if (this->m_CellFusionMode == otb::CellFusionMode::ACC)
764  {
765  outputDEMIt.Set(static_cast<DEMPixelType>(accPixel));
766  }
767  }
768  ++outputDEMIt;
769  ++firstDEMIt;
770  ++firstDEMAccIt;
771  }
772 
773  // Check DEMs from other threads and keep the maximum elevation
774  for (unsigned int i = 1; i < m_TempDEMRegions.size(); i++)
775  {
776 
777  itk::ImageRegionIterator<OutputImageType> tmpDEMIt(m_TempDEMRegions[i], outputDEM->GetRequestedRegion());
778  itk::ImageRegionIterator<AccumulatorImageType> tmpDEMAccIt(m_TempDEMAccumulatorRegions[i], outputDEM->GetRequestedRegion());
779 
780  outputDEMIt.GoToBegin();
781  tmpDEMIt.GoToBegin();
782  tmpDEMAccIt.GoToBegin();
783  firstDEMAccIt.GoToBegin(); // Global Accumulator
784  while (!outputDEMIt.IsAtEnd() && !tmpDEMIt.IsAtEnd() && !tmpDEMAccIt.IsAtEnd() && !firstDEMAccIt.IsAtEnd())
785  {
786 
787  // get the accumulator value
788  AccumulatorPixelType accPixel = tmpDEMAccIt.Get();
789  if (accPixel != 0)
790  {
791 
792  DEMPixelType cellCurrentValue = outputDEMIt.Get();
793  DEMPixelType cellHeight = tmpDEMIt.Get();
794  switch (this->m_CellFusionMode)
795  {
796  case otb::CellFusionMode::MIN:
797  {
798  if ((cellHeight < cellCurrentValue) || (cellCurrentValue == m_NoDataValue))
799  {
800  outputDEMIt.Set(cellHeight);
801  }
802  }
803  break;
804  case otb::CellFusionMode::MAX:
805  {
806  if ((cellHeight > cellCurrentValue) || ((cellCurrentValue == m_NoDataValue)))
807  {
808  outputDEMIt.Set(cellHeight);
809  }
810  }
811  break;
812  case otb::CellFusionMode::MEAN:
813  {
814 
815  outputDEMIt.Set(cellCurrentValue * static_cast<DEMPixelType>(cellCurrentValue != m_NoDataValue) + cellHeight);
816  firstDEMAccIt.Set(firstDEMAccIt.Get() + accPixel);
817  }
818  break;
819  case otb::CellFusionMode::ACC:
820  {
821  firstDEMAccIt.Set(firstDEMAccIt.Get() + accPixel);
822  }
823  break;
824 
825  default:
826  itkExceptionMacro(<< "Unexpected value cell fusion mode :" << this->m_CellFusionMode);
827  break;
828  }
829  }
830 
831  if (i == (m_TempDEMRegions.size() - 1))
832  {
833  if (this->m_CellFusionMode == otb::CellFusionMode::MEAN)
834  {
835  if (static_cast<DEMPixelType>(firstDEMAccIt.Get()) != 0)
836  {
837  outputDEMIt.Set(outputDEMIt.Get() / static_cast<DEMPixelType>(firstDEMAccIt.Get()));
838  }
839  else
840  {
841  outputDEMIt.Set(m_NoDataValue);
842  }
843  }
844  else if (this->m_CellFusionMode == otb::CellFusionMode::ACC)
845  {
846  outputDEMIt.Set(static_cast<DEMPixelType>(firstDEMAccIt.Get()));
847  }
848  }
849  ++outputDEMIt;
850  ++tmpDEMIt;
851  ++tmpDEMAccIt;
852  ++firstDEMAccIt;
853  }
854  }
855 }
856 }
857 
858 
859 #endif
otb::Multi3DMapToDEMFilter::BeforeThreadedGenerateData
void BeforeThreadedGenerateData() override
Definition: otbMulti3DMapToDEMFilter.hxx:466
otb::GenericRSTransform::Pointer
itk::SmartPointer< Self > Pointer
Definition: otbGenericRSTransform.h:66
otb::Multi3DMapToDEMFilter::GenerateOutputInformation
void GenerateOutputInformation() override
Definition: otbMulti3DMapToDEMFilter.hxx:288
otb::Multi3DMapToDEMFilter::SetNumberOf3DMaps
void SetNumberOf3DMaps(unsigned int nb)
Definition: otbMulti3DMapToDEMFilter.hxx:64
otbStreamingStatisticsVectorImageFilter.h
otb::CellFusionMode::MIN
@ MIN
Definition: otbMulti3DMapToDEMFilter.h:41
otb::Multi3DMapToDEMFilter::GetMaskInput
const TMaskImage * GetMaskInput(unsigned int index) const
Definition: otbMulti3DMapToDEMFilter.hxx:112
otb::Multi3DMapToDEMFilter::Multi3DMapToDEMFilter
Multi3DMapToDEMFilter()
Definition: otbMulti3DMapToDEMFilter.hxx:34
otb::WriteNoDataFlags
void OTBMetadata_EXPORT WriteNoDataFlags(const std::vector< bool > &flags, const std::vector< double > &values, ImageMetadata &imd)
otb::CellFusionMode::ACC
@ ACC
Definition: otbMulti3DMapToDEMFilter.h:44
otb::CellFusionMode::MEAN
@ MEAN
Definition: otbMulti3DMapToDEMFilter.h:43
otb::Multi3DMapToDEMFilter::~Multi3DMapToDEMFilter
~Multi3DMapToDEMFilter() override
Definition: otbMulti3DMapToDEMFilter.hxx:59
otb
The "otb" namespace contains all Orfeo Toolbox (OTB) classes.
Definition: otbJoinContainer.h:32
otb::Image
Creation of an "otb" image which contains metadata.
Definition: otbImage.h:89
otb::Multi3DMapToDEMFilter::AccumulatorPixelType
AccumulatorImageType::PixelType AccumulatorPixelType
Definition: otbMulti3DMapToDEMFilter.h:111
otb::ImageToGenericRSOutputParameters::Pointer
itk::SmartPointer< Self > Pointer
Definition: otbImageToGenericRSOutputParameters.h:63
otb::GenericRSTransform::OutputPointType
itk::Point< ScalarType, NOutputDimensions > OutputPointType
Definition: otbGenericRSTransform.h:73
otbMulti3DMapToDEMFilter.h
otb::Multi3DMapToDEMFilter::DEMPixelType
OutputImageType::PixelType DEMPixelType
Definition: otbMulti3DMapToDEMFilter.h:105
otb::Multi3DMapToDEMFilter::RegionType
OutputImageType::RegionType RegionType
Definition: otbMulti3DMapToDEMFilter.h:104
otbNoDataHelper.h
otb::Multi3DMapToDEMFilter::TDPointType
RSTransformType::InputPointType TDPointType
Definition: otbMulti3DMapToDEMFilter.h:123
otb::Multi3DMapToDEMFilter::Set3DMapInput
void Set3DMapInput(unsigned int index, const T3DImage *hmap)
Definition: otbMulti3DMapToDEMFilter.hxx:79
otb::Image::GetSignedSpacing
SpacingType GetSignedSpacing() const
Definition: otbImage.hxx:38
otb::CellFusionMode::MAX
@ MAX
Definition: otbMulti3DMapToDEMFilter.h:42
otb::GenericRSTransform::InputPointType
itk::Point< ScalarType, NInputDimensions > InputPointType
Definition: otbGenericRSTransform.h:72
otb::Multi3DMapToDEMFilter::GetNumberOf3DMaps
unsigned int GetNumberOf3DMaps()
Definition: otbMulti3DMapToDEMFilter.hxx:73
otb::Multi3DMapToDEMFilter::Get3DMapInput
const T3DImage * Get3DMapInput(unsigned int index) const
Definition: otbMulti3DMapToDEMFilter.hxx:102
otb::Multi3DMapToDEMFilter::GenerateInputRequestedRegion
void GenerateInputRequestedRegion() override
Definition: otbMulti3DMapToDEMFilter.hxx:338
otbMsgDevMacro
#define otbMsgDevMacro(x)
Definition: otbMacro.h:64
otb::Multi3DMapToDEMFilter::MapPixelType
InputMapType::PixelType MapPixelType
Definition: otbMulti3DMapToDEMFilter.h:113
otb::Multi3DMapToDEMFilter::SetOutputParametersFromImage
void SetOutputParametersFromImage()
Definition: otbMulti3DMapToDEMFilter.hxx:143
otb::Multi3DMapToDEMFilter::GetDEMOutput
const TOutputDEMImage * GetDEMOutput() const
Definition: otbMulti3DMapToDEMFilter.hxx:122
otb::SpatialReference::FromWGS84
static SpatialReference FromWGS84()
otb::Multi3DMapToDEMFilter::ThreadedGenerateData
void ThreadedGenerateData(const RegionType &outputRegionForThread, itk::ThreadIdType threadId) override
Definition: otbMulti3DMapToDEMFilter.hxx:529
otb::Image::Pointer
itk::SmartPointer< Self > Pointer
Definition: otbImage.h:97
otb::MetaDataKey::ProjectionRefKey
OTBMetadata_EXPORT char const * ProjectionRefKey
otb::Multi3DMapToDEMFilter::SetMaskInput
void SetMaskInput(unsigned int index, const TMaskImage *mask)
Definition: otbMulti3DMapToDEMFilter.hxx:91
otb::Multi3DMapToDEMFilter::AfterThreadedGenerateData
void AfterThreadedGenerateData() override
Definition: otbMulti3DMapToDEMFilter.hxx:721
Generated at Wed Feb 28 2024 15:23:05 for OTB with   doxygen 1.8.17