Orfeo Toolbox  3.16
otbSIXSTraits.cxx
Go to the documentation of this file.
1 /*=========================================================================
2 
3  Program: ORFEO Toolbox
4  Language: C++
5  Date: 2008-05-30
6  Version: 2.2.0
7 
8 
9  Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
10  See OTBCopyright.txt for details.
11 
12 
13  This software is distributed WITHOUT ANY WARRANTY; without even
14  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
15  PURPOSE. See the above copyright notices for more information.
16 
17 =========================================================================*/
18 #include "otbSIXSTraits.h"
19 
20 #include "otb_6S.h"
21 #include "main_6s.h"
22 #include "otbMacro.h"
23 
24 #include <iomanip>
25 
26 namespace otb
27 {
28 
29 void
30 SIXSTraits::ComputeAtmosphericParameters(
31  const double SolarZenithalAngle,
32  const double SolarAzimutalAngle,
33  const double ViewingZenithalAngle,
34  const double ViewingAzimutalAngle,
35  const unsigned int Month,
36  const unsigned int Day,
37  const double AtmosphericPressure,
38  const double WaterVaporAmount,
39  const double OzoneAmount,
40  const AerosolModelType& AerosolModel,
41  const double AerosolOptical,
42  WavelengthSpectralType* WavelengthSpectralBand,
44  double& AtmosphericReflectance,
45  double& AtmosphericSphericalAlbedo,
46  double& TotalGaseousTransmission,
47  double& DownwardTransmittance,
48  double& UpwardTransmittance,
49  double& UpwardDiffuseTransmittance,
50  double& UpwardDirectTransmittance,
51  double& UpwardDiffuseTransmittanceForRayleigh,
52  double& UpwardDiffuseTransmittanceForAerosol
53  )
54 {
55 // geometrical conditions
56  otb_6s_doublereal asol(static_cast<otb_6s_doublereal>(SolarZenithalAngle));
57  otb_6s_doublereal phi0(static_cast<otb_6s_doublereal>(SolarAzimutalAngle));
58  otb_6s_doublereal avis(static_cast<otb_6s_doublereal>(ViewingZenithalAngle));
59  otb_6s_doublereal phiv(static_cast<otb_6s_doublereal>(ViewingAzimutalAngle));
60  otb_6s_integer month(static_cast<otb_6s_integer>(Month));
61  otb_6s_integer jday(static_cast<otb_6s_integer>(Day));
62  otb_6s_doublereal pressure(static_cast<otb_6s_doublereal>(AtmosphericPressure));
63  otb_6s_doublereal uw(static_cast<otb_6s_doublereal>(WaterVaporAmount));
64  otb_6s_doublereal uo3(static_cast<otb_6s_doublereal>(OzoneAmount));
65 // atmospheric model
66  otb_6s_integer iaer(static_cast<otb_6s_integer>(AerosolModel));
67  otb_6s_doublereal taer55(static_cast<otb_6s_doublereal>(AerosolOptical));
68 
69  // Init output parameters
70  AtmosphericReflectance = 0.;
71  AtmosphericSphericalAlbedo = 0.;
72  TotalGaseousTransmission = 0.;
73  DownwardTransmittance = 0.;
74  UpwardTransmittance = 0.;
75  UpwardDiffuseTransmittance = 0.;
76  UpwardDirectTransmittance = 0.;
77  UpwardDiffuseTransmittanceForRayleigh = 0.;
78  UpwardDiffuseTransmittanceForAerosol = 0.;
79 
80  otb_6s_doublereal wlinf(0.), wlsup(0.);
81  otb_6s_doublereal otb_ratm__(0.), sast(0.), tgasm(0.), sdtott(0.), sutott(0.);
82  otb_6s_doublereal tdif_up(0.), tdir_up(0.), tdif_up_ray(0.), tdif_up_aer(0.);
83 
84  // 6S official Wavelength Spectral Band step value
85  const float SIXSStepOfWavelengthSpectralBandValues = .0025;
86  // Generate 6s Wavelength Spectral Band with the offcicial step value
87  ComputeWavelengthSpectralBandValuesFor6S(SIXSStepOfWavelengthSpectralBandValues,
88  WavelengthSpectralBand // Update
89  );
90  try
91  {
92 
93  // 6S official tab size Wavelength Spectral
94  const otb_6s_integer S_6S_SIZE = 1501;
95  // Generate WavelengthSpectralBand in 6S compatible buffer s[1501]
96  wlinf = static_cast<otb_6s_doublereal>(WavelengthSpectralBand->GetMinSpectralValue());
97  wlsup = static_cast<otb_6s_doublereal>(WavelengthSpectralBand->GetMaxSpectralValue());
98  otb_6s_integer iinf =
99  static_cast<otb_6s_integer>((wlinf - (float) .25) / SIXSStepOfWavelengthSpectralBandValues + (float) 1.5);
100  otb_6s_integer cpt = iinf - 1;
101  otb_6s_doublereal * s(NULL);
102  s = new otb_6s_doublereal[S_6S_SIZE];
103  memset(s, 0, S_6S_SIZE * sizeof(otb_6s_doublereal));
104  const ValuesVectorType& FilterFunctionValues6S = WavelengthSpectralBand->GetFilterFunctionValues6S();
105  // Set the values of FilterFunctionValues6S in s between [iinf-1; isup]
106  for (unsigned int i = 0; i < FilterFunctionValues6S.size() && cpt < S_6S_SIZE; ++i)
107  {
108  s[cpt] = FilterFunctionValues6S[i];
109  ++cpt;
110  }
111 
112  // Call 6s main function
113  otbMsgDevMacro(<< "Start call 6S main function ...");
114  otb_6s_ssssss_otb_main_function(&asol, &phi0, &avis, &phiv, &month, &jday,
115  &pressure, &uw, &uo3,
116  &iaer,
117  &taer55,
118  &wlinf, &wlsup,
119  s,
120  &otb_ratm__,
121  &sast,
122  &tgasm,
123  &sdtott,
124  &sutott,
125  &tdif_up,
126  &tdir_up,
127  &tdif_up_ray,
128  &tdif_up_aer);
129  otbMsgDevMacro(<< "Done call 6S main function!");
130  delete[] s;
131  s = NULL;
132  }
133  catch (std::bad_alloc& err)
134  {
135  itkGenericExceptionMacro(<< "Exception bad_alloc in SIXSTraits class: " << (char*) err.what());
136  }
137  catch (...)
138  {
139  itkGenericExceptionMacro(<< "Unknown exception in SIXSTraits class (catch(...)");
140  }
141 
142  // Set outputs parameters
143  AtmosphericReflectance = static_cast<double>(otb_ratm__);
144  AtmosphericSphericalAlbedo = static_cast<double>(sast);
145  TotalGaseousTransmission = static_cast<double>(tgasm);
146  DownwardTransmittance = static_cast<double>(sdtott);
147  UpwardTransmittance = static_cast<double>(sutott);
148  UpwardDiffuseTransmittance = static_cast<double>(tdif_up);
149  UpwardDirectTransmittance = static_cast<double>(tdir_up);
150  UpwardDiffuseTransmittanceForRayleigh = static_cast<double>(tdif_up_ray);
151  UpwardDiffuseTransmittanceForAerosol = static_cast<double>(tdif_up_aer);
152 }
153 
154 void
155 SIXSTraits::ComputeWavelengthSpectralBandValuesFor6S(
156  const double SIXSStepOfWavelengthSpectralBandValues,
157  WavelengthSpectralType* WavelengthSpectralBand
158  )
159 {
160  const double epsilon(.000001);
161  const double L_min = static_cast<double>(WavelengthSpectralBand->GetMinSpectralValue());
162  const double L_max = static_cast<double>(WavelengthSpectralBand->GetMaxSpectralValue());
163  const double L_userStep = static_cast<double>(WavelengthSpectralBand->GetUserStep());
164  const ValuesVectorType& FilterFunctionValues = WavelengthSpectralBand->GetFilterFunctionValues();
165  unsigned int i = 1;
166  unsigned int j = 1;
167  const double invStep = static_cast<double>(1. / L_userStep);
168  double value(0.);
169 
170  if (FilterFunctionValues.size() <= 1)
171  {
172  itkGenericExceptionMacro(<< "The FilterFunctionValues vector must have more than 1 values !");
173  }
174  if (! (L_min + static_cast<double>(FilterFunctionValues.size() - 1) * L_userStep < (L_max + epsilon) ) )
175  {
176  itkGenericExceptionMacro(
177  << "The following condition: " << L_min << "+(" << FilterFunctionValues.size() << "-1)*" << L_userStep <<
178  " < (" << L_max << "+" << epsilon << ") is not respected !" << "val1 " << L_min + static_cast<double>(FilterFunctionValues.size() - 1) * L_userStep << " val2 " << L_max - epsilon);
179  }
180 
181  // Generate WavelengthSpectralBand if the step is not the offical 6S step value
182  if (vcl_abs(L_userStep - SIXSStepOfWavelengthSpectralBandValues) > epsilon)
183  {
184  ValuesVectorType values(1, FilterFunctionValues[0]); //vector size 1 with the value vect[0]
185 
186  // Stop the interpolation at the max spectral value.
187  value = SIXSStepOfWavelengthSpectralBandValues;
188  while (L_min + value <= L_max)
189  {
190  // Search the User interval that surround the StepOfWavelengthSpectralBandValues current value.
191 
192  // removed the <= here, might be wrong
193  while (j * L_userStep < value)
194  {
195  ++j;
196  }
197 
198  // Check if we are not out of bound
199  if (j >= FilterFunctionValues.size())
200  {
201  itkGenericExceptionMacro(
202  << "Index " << j << " out of bound for FilterFunctionValues vector (size: " << FilterFunctionValues.size() <<
203  ")." << " and value is " << value << " and SIXSStepOfWavelengthSpectralBandValues is " << SIXSStepOfWavelengthSpectralBandValues<< " and i is " << i);
204  }
205 
206  double valueTemp;
207  valueTemp = static_cast<double>(FilterFunctionValues[j - 1])
208  + ((static_cast<double>(FilterFunctionValues[j]) -
209  static_cast<double>(FilterFunctionValues[j - 1])) * invStep)
210  * (value - L_userStep * (j - 1));
211  values.push_back(static_cast<WavelengthSpectralBandType>(valueTemp));
212 
213  ++i;
214  value += SIXSStepOfWavelengthSpectralBandValues;
215  }
216 
217  if (L_min + (i - 1) * SIXSStepOfWavelengthSpectralBandValues != L_max)
218  {
219  values.push_back(0);
220  }
221  // Store this values
222  WavelengthSpectralBand->SetFilterFunctionValues6S(values);
223  // Store the new Max MaxSpectralValue
224  WavelengthSpectralBand->SetMaxSpectralValue(static_cast<WavelengthSpectralBandType>(L_min + i *
225  SIXSStepOfWavelengthSpectralBandValues));
226 
227  }
228  else
229  {
230  // Init with copy of FilterFunctionValues input vector values
231  WavelengthSpectralBand->SetFilterFunctionValues6S(FilterFunctionValues);
232  }
233 }
234 
235 void
236 SIXSTraits::ComputeEnvironmentalContribution(const double diffuseTransmittanceForRayleighScattering,
237  const double diffuseTransmittanceForAerosolScattering,
238  const double radiusInKilometers,
239  const double altitude,
240  const double cosineOfViewingAngle,
241  double& rayleighEstimation,
242  double& aerosolEstimation,
243  double& globalEstimation)
244 {
245  otb_6s_doublereal difr(static_cast<otb_6s_doublereal>(diffuseTransmittanceForRayleighScattering));
246  otb_6s_doublereal difa(static_cast<otb_6s_doublereal>(diffuseTransmittanceForAerosolScattering));
247  otb_6s_doublereal rad(static_cast<otb_6s_doublereal>(radiusInKilometers));
248  otb_6s_doublereal palt(static_cast<otb_6s_doublereal>(altitude));
249  otb_6s_doublereal xmuv(static_cast<otb_6s_doublereal>(cosineOfViewingAngle));
250  otb_6s_doublereal fra(0.), fae(0.), fr(0.);
251  otb_6s_enviro_(&difr, &difa, &rad, &palt, &xmuv, &fra, &fae, &fr);
252  rayleighEstimation = static_cast<double>(fra);
253  aerosolEstimation = static_cast<double>(fae);
254  globalEstimation = static_cast<double>(fr);
255 }
256 
257 } // namespace otb
Generated at Sun Feb 3 2013 00:47:38 for Orfeo Toolbox with doxygen 1.8.1.1