OTB  6.7.0
Orfeo Toolbox
otbDEMCaracteristicsExtractor.hxx
Go to the documentation of this file.
1 /*
2  * Copyright (C) 1999-2011 Insight Software Consortium
3  * Copyright (C) 2005-2019 Centre National d'Etudes Spatiales (CNES)
4  *
5  * This file is part of Orfeo Toolbox
6  *
7  * https://www.orfeo-toolbox.org/
8  *
9  * Licensed under the Apache License, Version 2.0 (the "License");
10  * you may not use this file except in compliance with the License.
11  * You may obtain a copy of the License at
12  *
13  * http://www.apache.org/licenses/LICENSE-2.0
14  *
15  * Unless required by applicable law or agreed to in writing, software
16  * distributed under the License is distributed on an "AS IS" BASIS,
17  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
18  * See the License for the specific language governing permissions and
19  * limitations under the License.
20  */
21 
22 #ifndef otbDEMCaracteristicsExtractor_hxx
23 #define otbDEMCaracteristicsExtractor_hxx
24 
26 #include "itkImageRegionIterator.h"
27 #include "otbMath.h"
28 
29 namespace otb
30 {
31 
32 template <class TInputImage, class TOutputImage>
35 {
36  this->SetNumberOfRequiredInputs(1);
37  this->SetNumberOfRequiredOutputs(3);
38 
39  this->SetNthOutput(0, OutputImageType::New());
40  this->SetNthOutput(1, OutputImageType::New());
41  this->SetNthOutput(2, OutputImageType::New());
42  this->SetNthOutput(3, OutputImageType::New());
43 
44  m_SolarAngle = 0;
45  m_SolarAzimut = 0;
46  m_ViewAngle = 0;
47  m_ViewAzimut = 0;
48 }
49 
50 template <class TInputImage, class TOutputImage>
53 {
54 }
55 
59 template <class TInputImage, class TOutputImage>
60 void
63 {
64  // Input and output pointer definition
65  typename InputImageType::Pointer inputPtr = const_cast<InputImageType *>(this->GetInput());
66  typename OutputImageType::Pointer SlopOutputPtr = this->GetSlopOutput();
67  typename OutputImageType::Pointer AspectOutputPtr = this->GetAspectOutput();
68  typename OutputImageType::Pointer IncidenceOutputPtr = this->GetIncidenceOutput();
69  typename OutputImageType::Pointer ExitanceOutputPtr = this->GetExitanceOutput();
71 
72  // Gradient Magnitude Image Filter used to compute the slope.
73  typename GradientMagnitudeFilterType::Pointer GradientMagnitudeFilter = GradientMagnitudeFilterType::New();
74  // Gradient Recursive Gaussian Image Filter used to compute the aspect.
75  typename GradientRecursiveGaussianImageFilterType::Pointer GradientRecursiveGaussianFilter =
76  GradientRecursiveGaussianImageFilterType::New();
77  // Atan used to compute the slop
78  typename AtanFilterType::Pointer AtanFilter = AtanFilterType::New();
79  // Atan2 Image Filter used to compute the aspect.
80  typename Atan2FilterType::Pointer AspectFilter = Atan2FilterType::New();
81  // Inverse cosinus Image filter used to compute the incidence image
82  typename AcosImageFilterType::Pointer IncidenceFilter = AcosImageFilterType::New();
83  // Inverse cosinus Image filter used to compute the exitance image
84  typename AcosImageFilterType::Pointer ExitanceFilter = AcosImageFilterType::New();
85 
86  // Degrees To Radian _-> Radian To Degree coefficient
87  double rad2degCoef;
88  rad2degCoef = 180 / CONST_PI;
89 
90  // Slop calculation
91  GradientMagnitudeFilter->SetInput(inputPtr);
92  AtanFilter->SetInput(GradientMagnitudeFilter->GetOutput());
93  // Transform values from radian to degrees.
94  typename MultiplyByScalarImageFilterType::Pointer rad2DegFilter = MultiplyByScalarImageFilterType::New();
95  //rad2DegFilter->SetInput( GradientMagnitudeFilter->GetOutput() );
96  rad2DegFilter->SetInput(AtanFilter->GetOutput());
97  rad2DegFilter->SetCoef(rad2degCoef);
98  rad2DegFilter->GraftOutput(SlopOutputPtr);
99  rad2DegFilter->Update();
100  this->GraftNthOutput(0, rad2DegFilter->GetOutput());
101 
102  // Aspect calcultation
103  GradientRecursiveGaussianFilter->SetInput(inputPtr);
104  GradientRecursiveGaussianFilter->Update();
105 
106  // // Extract the X and the Y gradient
107  typename AdaptorType::Pointer XAdaptator = AdaptorType::New();
108  typename AdaptorType::Pointer YAdaptator = AdaptorType::New();
109  XAdaptator->SetImage(GradientRecursiveGaussianFilter->GetOutput());
110  YAdaptator->SetImage(GradientRecursiveGaussianFilter->GetOutput());
111  XAdaptator->SelectNthElement(0);
112  YAdaptator->SelectNthElement(1);
113  // // Compute Arctan
114  AspectFilter->SetInput1(XAdaptator);
115  AspectFilter->SetInput2(YAdaptator);
116  // // Transform values from radian to degres.
117  typename MultiplyByScalarImageFilterType::Pointer rad2DegFilter1 = MultiplyByScalarImageFilterType::New();
118  rad2DegFilter1->SetInput(AspectFilter->GetOutput());
119  rad2DegFilter1->SetCoef(rad2degCoef);
120  rad2DegFilter1->GraftOutput(AspectOutputPtr);
121  rad2DegFilter1->Update();
122  this->GraftNthOutput(1, rad2DegFilter1->GetOutput());
123 
124  // Angle calculation :
125  // sin(slop)
126  typename SinImageFilterType::Pointer sinS = SinImageFilterType::New();
127  sinS->SetInput(GradientMagnitudeFilter->GetOutput());
128  // cos (slop)
129  typename CosImageFilterType::Pointer cosS = CosImageFilterType::New();
130  cosS->SetInput(GradientMagnitudeFilter->GetOutput());
131  // -aspect
132  typename MultiplyByScalarImageFilterType::Pointer oppositeFilter = MultiplyByScalarImageFilterType::New();
133  oppositeFilter->SetInput(AspectFilter->GetOutput());
134  oppositeFilter->SetCoef(-1);
135 
136  // Incidence calculation
137  typename ShiftScaleImageFilterType::Pointer addAzimut = ShiftScaleImageFilterType::New();
138  addAzimut->SetScale(1.);
139  addAzimut->SetShift(m_SolarAzimut / rad2degCoef);
140  addAzimut->SetInput(oppositeFilter->GetOutput());
141 
142  typename CosImageFilterType::Pointer cosAAzimut = CosImageFilterType::New();
143  cosAAzimut->SetInput(addAzimut->GetOutput());
144 
145  typename MultiplyByScalarImageFilterType::Pointer sinSsinSolarAngleFilter = MultiplyByScalarImageFilterType::New();
146  sinSsinSolarAngleFilter->SetCoef(std::sin(m_SolarAngle / rad2degCoef));
147  sinSsinSolarAngleFilter->SetInput(sinS->GetOutput());
148 
149  typename MultiplyImageFilterType::Pointer cosAAzimuthsinSsinAngle = MultiplyImageFilterType::New();
150  cosAAzimuthsinSsinAngle->SetInput1(sinSsinSolarAngleFilter->GetOutput());
151  cosAAzimuthsinSsinAngle->SetInput2(cosAAzimut->GetOutput());
152 
153  typename MultiplyByScalarImageFilterType::Pointer cosScosSolarAngleFilter = MultiplyByScalarImageFilterType::New();
154  cosScosSolarAngleFilter->SetCoef(std::cos(m_SolarAngle / rad2degCoef));
155  cosScosSolarAngleFilter->SetInput(cosS->GetOutput());
156 
157  typename AddImageFilterType::Pointer cosIncidence = AddImageFilterType::New();
158  cosIncidence->SetInput1(cosAAzimuthsinSsinAngle->GetOutput());
159  cosIncidence->SetInput2(cosScosSolarAngleFilter->GetOutput());
160 
161  IncidenceFilter->SetInput(cosIncidence->GetOutput());
162 
163  // // Change radians in degres
164  typename MultiplyByScalarImageFilterType::Pointer rad2DegFilter2 = MultiplyByScalarImageFilterType::New();
165  rad2DegFilter2->SetInput(IncidenceFilter->GetOutput());
166  rad2DegFilter2->SetCoef(rad2degCoef);
167  // // Link to the output
168  rad2DegFilter2->GraftOutput(IncidenceOutputPtr);
169  rad2DegFilter2->Update();
170  this->GraftNthOutput(2, rad2DegFilter2->GetOutput());
171 
172  // Exitance calculation
173  typename ShiftScaleImageFilterType::Pointer addAzimut2 = ShiftScaleImageFilterType::New();
174  addAzimut2->SetScale(1.);
175  addAzimut2->SetShift(m_ViewAzimut / rad2degCoef);
176  addAzimut2->SetInput(oppositeFilter->GetOutput());
177 
178  typename CosImageFilterType::Pointer cosAAzimut2 = CosImageFilterType::New();
179  cosAAzimut2->SetInput(addAzimut2->GetOutput());
180 
181  typename MultiplyByScalarImageFilterType::Pointer sinSsinSolarAngleFilter2 = MultiplyByScalarImageFilterType::New();
182  sinSsinSolarAngleFilter2->SetCoef(std::sin(m_ViewAngle / rad2degCoef));
183  sinSsinSolarAngleFilter2->SetInput(sinS->GetOutput());
184 
185  typename MultiplyImageFilterType::Pointer cosAAzimuthsinSsinAngle2 = MultiplyImageFilterType::New();
186  cosAAzimuthsinSsinAngle2->SetInput1(sinSsinSolarAngleFilter2->GetOutput());
187  cosAAzimuthsinSsinAngle2->SetInput2(cosAAzimut2->GetOutput());
188 
189  typename MultiplyByScalarImageFilterType::Pointer cosScosSolarAngleFilter2 = MultiplyByScalarImageFilterType::New();
190  cosScosSolarAngleFilter2->SetCoef(std::cos(m_ViewAngle / rad2degCoef));
191  cosScosSolarAngleFilter2->SetInput(cosS->GetOutput());
192 
193  typename AddImageFilterType::Pointer cosIncidence2 = AddImageFilterType::New();
194  cosIncidence2->SetInput1(cosAAzimuthsinSsinAngle2->GetOutput());
195  cosIncidence2->SetInput2(cosScosSolarAngleFilter2->GetOutput());
196 
197  ExitanceFilter->SetInput(cosIncidence2->GetOutput());
198 
199  // // Change radians in degres
200  typename MultiplyByScalarImageFilterType::Pointer rad2DegFilter3 = MultiplyByScalarImageFilterType::New();
201  rad2DegFilter3->SetInput(ExitanceFilter->GetOutput());
202  rad2DegFilter3->SetCoef(rad2degCoef);
203  // // Link to the output
204  rad2DegFilter3->GraftOutput(ExitanceOutputPtr);
205  rad2DegFilter3->Update();
206  this->GraftNthOutput(3, rad2DegFilter3->GetOutput());
207 
208 }
209 
211 template <class TInputImage, class TOutputImage>
212 void
214 ::PrintSelf(std::ostream& os, itk::Indent indent) const
215 {
216  Superclass::PrintSelf(os, indent);
217  os << indent << "Solar Angle: " << m_SolarAngle << std::endl;
218  os << indent << "Solar Azimut: " << m_SolarAzimut << std::endl;
219  os << indent << "View Angle: " << m_ViewAngle << std::endl;
220  os << indent << "View Azimut: " << m_ViewAzimut << std::endl;
221 }
223 
224 } // end namespace otb
225 
226 #endif
constexpr double CONST_PI
Definition: otbMath.h:48
void PrintSelf(std::ostream &os, itk::Indent indent) const override
TInputImage InputImageType