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