OTB  7.2.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-2020 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>
34 {
35  this->SetNumberOfRequiredInputs(1);
36  this->SetNumberOfRequiredOutputs(3);
37 
38  this->SetNthOutput(0, OutputImageType::New());
39  this->SetNthOutput(1, OutputImageType::New());
40  this->SetNthOutput(2, OutputImageType::New());
41  this->SetNthOutput(3, OutputImageType::New());
42 
43  m_SolarAngle = 0;
44  m_SolarAzimut = 0;
45  m_ViewAngle = 0;
46  m_ViewAzimut = 0;
47 }
48 
49 template <class TInputImage, class TOutputImage>
51 {
52 }
53 
57 template <class TInputImage, class TOutputImage>
59 {
60  // Input and output pointer definition
61  typename InputImageType::Pointer inputPtr = const_cast<InputImageType*>(this->GetInput());
62  typename OutputImageType::Pointer SlopOutputPtr = this->GetSlopOutput();
63  typename OutputImageType::Pointer AspectOutputPtr = this->GetAspectOutput();
64  typename OutputImageType::Pointer IncidenceOutputPtr = this->GetIncidenceOutput();
65  typename OutputImageType::Pointer ExitanceOutputPtr = this->GetExitanceOutput();
67 
68  // Gradient Magnitude Image Filter used to compute the slope.
69  typename GradientMagnitudeFilterType::Pointer GradientMagnitudeFilter = GradientMagnitudeFilterType::New();
70  // Gradient Recursive Gaussian Image Filter used to compute the aspect.
71  typename GradientRecursiveGaussianImageFilterType::Pointer GradientRecursiveGaussianFilter = GradientRecursiveGaussianImageFilterType::New();
72  // Atan used to compute the slop
73  typename AtanFilterType::Pointer AtanFilter = AtanFilterType::New();
74  // Atan2 Image Filter used to compute the aspect.
75  typename Atan2FilterType::Pointer AspectFilter = Atan2FilterType::New();
76  // Inverse cosinus Image filter used to compute the incidence image
77  typename AcosImageFilterType::Pointer IncidenceFilter = AcosImageFilterType::New();
78  // Inverse cosinus Image filter used to compute the exitance image
79  typename AcosImageFilterType::Pointer ExitanceFilter = AcosImageFilterType::New();
80 
81  // Degrees To Radian _-> Radian To Degree coefficient
82  double rad2degCoef;
83  rad2degCoef = 180 / CONST_PI;
84 
85  // Slop calculation
86  GradientMagnitudeFilter->SetInput(inputPtr);
87  AtanFilter->SetInput(GradientMagnitudeFilter->GetOutput());
88  // Transform values from radian to degrees.
89  typename MultiplyByScalarImageFilterType::Pointer rad2DegFilter = MultiplyByScalarImageFilterType::New();
90  // rad2DegFilter->SetInput( GradientMagnitudeFilter->GetOutput() );
91  rad2DegFilter->SetInput(AtanFilter->GetOutput());
92  rad2DegFilter->SetCoef(rad2degCoef);
93  rad2DegFilter->GraftOutput(SlopOutputPtr);
94  rad2DegFilter->Update();
95  this->GraftNthOutput(0, rad2DegFilter->GetOutput());
96 
97  // Aspect calcultation
98  GradientRecursiveGaussianFilter->SetInput(inputPtr);
99  GradientRecursiveGaussianFilter->Update();
100 
101  // // Extract the X and the Y gradient
102  typename AdaptorType::Pointer XAdaptator = AdaptorType::New();
103  typename AdaptorType::Pointer YAdaptator = AdaptorType::New();
104  XAdaptator->SetImage(GradientRecursiveGaussianFilter->GetOutput());
105  YAdaptator->SetImage(GradientRecursiveGaussianFilter->GetOutput());
106  XAdaptator->SelectNthElement(0);
107  YAdaptator->SelectNthElement(1);
108  // // Compute Arctan
109  AspectFilter->SetInput1(XAdaptator);
110  AspectFilter->SetInput2(YAdaptator);
111  // // Transform values from radian to degres.
112  typename MultiplyByScalarImageFilterType::Pointer rad2DegFilter1 = MultiplyByScalarImageFilterType::New();
113  rad2DegFilter1->SetInput(AspectFilter->GetOutput());
114  rad2DegFilter1->SetCoef(rad2degCoef);
115  rad2DegFilter1->GraftOutput(AspectOutputPtr);
116  rad2DegFilter1->Update();
117  this->GraftNthOutput(1, rad2DegFilter1->GetOutput());
118 
119  // Angle calculation :
120  // sin(slop)
121  typename SinImageFilterType::Pointer sinS = SinImageFilterType::New();
122  sinS->SetInput(GradientMagnitudeFilter->GetOutput());
123  // cos (slop)
124  typename CosImageFilterType::Pointer cosS = CosImageFilterType::New();
125  cosS->SetInput(GradientMagnitudeFilter->GetOutput());
126  // -aspect
127  typename MultiplyByScalarImageFilterType::Pointer oppositeFilter = MultiplyByScalarImageFilterType::New();
128  oppositeFilter->SetInput(AspectFilter->GetOutput());
129  oppositeFilter->SetCoef(-1);
130 
131  // Incidence calculation
132  typename ShiftScaleImageFilterType::Pointer addAzimut = ShiftScaleImageFilterType::New();
133  addAzimut->SetScale(1.);
134  addAzimut->SetShift(m_SolarAzimut / rad2degCoef);
135  addAzimut->SetInput(oppositeFilter->GetOutput());
136 
137  typename CosImageFilterType::Pointer cosAAzimut = CosImageFilterType::New();
138  cosAAzimut->SetInput(addAzimut->GetOutput());
139 
140  typename MultiplyByScalarImageFilterType::Pointer sinSsinSolarAngleFilter = MultiplyByScalarImageFilterType::New();
141  sinSsinSolarAngleFilter->SetCoef(std::sin(m_SolarAngle / rad2degCoef));
142  sinSsinSolarAngleFilter->SetInput(sinS->GetOutput());
143 
144  typename MultiplyImageFilterType::Pointer cosAAzimuthsinSsinAngle = MultiplyImageFilterType::New();
145  cosAAzimuthsinSsinAngle->SetInput1(sinSsinSolarAngleFilter->GetOutput());
146  cosAAzimuthsinSsinAngle->SetInput2(cosAAzimut->GetOutput());
147 
148  typename MultiplyByScalarImageFilterType::Pointer cosScosSolarAngleFilter = MultiplyByScalarImageFilterType::New();
149  cosScosSolarAngleFilter->SetCoef(std::cos(m_SolarAngle / rad2degCoef));
150  cosScosSolarAngleFilter->SetInput(cosS->GetOutput());
151 
152  typename AddImageFilterType::Pointer cosIncidence = AddImageFilterType::New();
153  cosIncidence->SetInput1(cosAAzimuthsinSsinAngle->GetOutput());
154  cosIncidence->SetInput2(cosScosSolarAngleFilter->GetOutput());
155 
156  IncidenceFilter->SetInput(cosIncidence->GetOutput());
157 
158  // // Change radians in degres
159  typename MultiplyByScalarImageFilterType::Pointer rad2DegFilter2 = MultiplyByScalarImageFilterType::New();
160  rad2DegFilter2->SetInput(IncidenceFilter->GetOutput());
161  rad2DegFilter2->SetCoef(rad2degCoef);
162  // // Link to the output
163  rad2DegFilter2->GraftOutput(IncidenceOutputPtr);
164  rad2DegFilter2->Update();
165  this->GraftNthOutput(2, rad2DegFilter2->GetOutput());
166 
167  // Exitance calculation
168  typename ShiftScaleImageFilterType::Pointer addAzimut2 = ShiftScaleImageFilterType::New();
169  addAzimut2->SetScale(1.);
170  addAzimut2->SetShift(m_ViewAzimut / rad2degCoef);
171  addAzimut2->SetInput(oppositeFilter->GetOutput());
172 
173  typename CosImageFilterType::Pointer cosAAzimut2 = CosImageFilterType::New();
174  cosAAzimut2->SetInput(addAzimut2->GetOutput());
175 
176  typename MultiplyByScalarImageFilterType::Pointer sinSsinSolarAngleFilter2 = MultiplyByScalarImageFilterType::New();
177  sinSsinSolarAngleFilter2->SetCoef(std::sin(m_ViewAngle / rad2degCoef));
178  sinSsinSolarAngleFilter2->SetInput(sinS->GetOutput());
179 
180  typename MultiplyImageFilterType::Pointer cosAAzimuthsinSsinAngle2 = MultiplyImageFilterType::New();
181  cosAAzimuthsinSsinAngle2->SetInput1(sinSsinSolarAngleFilter2->GetOutput());
182  cosAAzimuthsinSsinAngle2->SetInput2(cosAAzimut2->GetOutput());
183 
184  typename MultiplyByScalarImageFilterType::Pointer cosScosSolarAngleFilter2 = MultiplyByScalarImageFilterType::New();
185  cosScosSolarAngleFilter2->SetCoef(std::cos(m_ViewAngle / rad2degCoef));
186  cosScosSolarAngleFilter2->SetInput(cosS->GetOutput());
187 
188  typename AddImageFilterType::Pointer cosIncidence2 = AddImageFilterType::New();
189  cosIncidence2->SetInput1(cosAAzimuthsinSsinAngle2->GetOutput());
190  cosIncidence2->SetInput2(cosScosSolarAngleFilter2->GetOutput());
191 
192  ExitanceFilter->SetInput(cosIncidence2->GetOutput());
193 
194  // // Change radians in degres
195  typename MultiplyByScalarImageFilterType::Pointer rad2DegFilter3 = MultiplyByScalarImageFilterType::New();
196  rad2DegFilter3->SetInput(ExitanceFilter->GetOutput());
197  rad2DegFilter3->SetCoef(rad2degCoef);
198  // // Link to the output
199  rad2DegFilter3->GraftOutput(ExitanceOutputPtr);
200  rad2DegFilter3->Update();
201  this->GraftNthOutput(3, rad2DegFilter3->GetOutput());
202 }
203 
205 template <class TInputImage, class TOutputImage>
206 void DEMCaracteristicsExtractor<TInputImage, TOutputImage>::PrintSelf(std::ostream& os, itk::Indent indent) const
207 {
208  Superclass::PrintSelf(os, indent);
209  os << indent << "Solar Angle: " << m_SolarAngle << std::endl;
210  os << indent << "Solar Azimut: " << m_SolarAzimut << std::endl;
211  os << indent << "View Angle: " << m_ViewAngle << std::endl;
212  os << indent << "View Azimut: " << m_ViewAzimut << std::endl;
213 }
215 
216 } // end namespace otb
217 
218 #endif
constexpr double CONST_PI
Definition: otbMath.h:49
void PrintSelf(std::ostream &os, itk::Indent indent) const override
The "otb" namespace contains all Orfeo Toolbox (OTB) classes.