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