Orfeo ToolBox  4.2
Orfeo ToolBox is not a black box
otbSailModel.cxx
Go to the documentation of this file.
1 /*=========================================================================
2 
3  Program: ORFEO Toolbox
4  Language: C++
5  Date: $Date$
6  Version: $Revision$
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 "itkNumericTraits.h"
19 
20 #include "otbSailModel.h"
21 #include <boost/math/special_functions/expint.hpp>
22 #include <boost/shared_ptr.hpp>
23 #include "otbMath.h"
24 
25 //TODO check EPSILON matlab
26 #define EPSILON 0.0000000000000000000000001
27 
28 namespace otb
29 {
30 
34 {
35  this->ProcessObject::SetNumberOfRequiredInputs(2);
36  this->ProcessObject::SetNumberOfRequiredOutputs(2);
37 
38  SpectralResponseType::Pointer vRefl = static_cast<SpectralResponseType *>(this->MakeOutput(0).GetPointer());
39  this->itk::ProcessObject::SetNthOutput(0, vRefl.GetPointer());
40 
41  SpectralResponseType::Pointer hRefl = static_cast<SpectralResponseType *>(this->MakeOutput(1).GetPointer());
42  this->itk::ProcessObject::SetNthOutput(1, hRefl.GetPointer());
43 
44  //default values
45  m_LAI=2;
46  m_Angl=50;
47  m_PSoil=1;
48  m_Skyl=70;
49  m_HSpot=0.2;
50  m_TTS=30;
51  m_TTO=0;
52  m_PSI=0;
53 }
54 
58 {}
59 
61 void
64 {
65  this->itk::ProcessObject::SetNthInput(0, const_cast<SpectralResponseType *>(object));
66 }
67 
71 {
72  if(this->GetNumberOfInputs() != 2)
73  {
74  //exit
75  return 0;
76  }
77  return static_cast<SpectralResponseType *>(this->itk::ProcessObject::GetInput(0));
78 }
79 
81 void
84 {
85  this->itk::ProcessObject::SetNthInput(1, const_cast<SpectralResponseType *>(object));
86 }
87 
91 {
92  if(this->GetNumberOfInputs() != 2)
93  {
94  //exit
95  return 0;
96  }
97  return static_cast<SpectralResponseType *>(this->itk::ProcessObject::GetInput(1));
98 }
99 
101 SailModel::DataObjectPointer
102 SailModel
103 ::MakeOutput(unsigned int)
104 {
105  return static_cast<itk::DataObject*>(SpectralResponseType::New().GetPointer());
106 }
107 
110 SailModel
112 {
113  if(this->GetNumberOfOutputs() < 2)
114  {
115  //exit
116  return 0;
117  }
118  return static_cast<SpectralResponseType *>(this->itk::ProcessObject::GetOutput(0));
119 }
120 
123 SailModel
125 {
126  if(this->GetNumberOfOutputs() < 2)
127  {
128  //exit
129  return 0;
130  }
131  return static_cast<SpectralResponseType *>(this->itk::ProcessObject::GetOutput(1));
132 }
133 
134 
136 void
137 SailModel
138 ::SetInput(const ParametersType & params)
139 {
140 
141  if(params.Size()!=8) itkExceptionMacro( << "Must have 8 parameters in that order : LAI, Angl, PSoil, Skyl, HSpot, TTS, TTO, PSI" );
142  this->SetParameters(params);
143  m_LAI=params[0];
144  m_Angl=params[1];
145  m_PSoil=params[2];
146  m_Skyl=params[3];
147  m_HSpot=params[4];
148  m_TTS=params[5];
149  m_TTO=params[6];
150  m_PSI=params[7];
151 }
152 
155 SailModel
157 {
158  ParametersType parameters=this->GetParameters();
159  if(parameters.Size()!=8)
160  {
161  parameters[0]=m_LAI;
162  parameters[1]=m_Angl;
163  parameters[2]=m_PSoil;
164  parameters[3]=m_Skyl;
165  parameters[4]=m_HSpot;
166  parameters[5]=m_TTS;
167  parameters[6]=m_TTO;
168  parameters[7]=m_PSI;
169  this->SetParameters(parameters);
170  }
171  return this->GetParameters();
172 }
173 
174 
176 void
177 SailModel
179 {
180 
181  SpectralResponseType::Pointer inRefl = this->GetReflectance();
182  SpectralResponseType::Pointer inTrans = this->GetTransmittance();
183  SpectralResponseType::Pointer outVRefl = this->GetViewingReflectance();
184  SpectralResponseType::Pointer outHRefl = this->GetHemisphericalReflectance();
185 
186  // LEAF ANGLE DISTRIBUTION
187  double rd = CONST_PI/180;
188  VectorType lidf;
189  this->Calc_LIDF(m_Angl, lidf);
190 
191  double cts, cto, ctscto, tants, tanto, cospsi, dso;
192  cts = vcl_cos(rd*m_TTS);
193  cto = vcl_cos(rd*m_TTO);
194  ctscto = cts*cto;
195  tants = vcl_tan(rd*m_TTS);
196  tanto = vcl_tan(rd*m_TTO);
197  cospsi = vcl_cos(rd*m_PSI);
198  dso = vcl_sqrt(tants*tants+tanto*tanto-2.*tants*tanto*cospsi);
199 
200  // angular distance, compensation of shadow length
201  // Calculate geometric factors associated with extinction and scattering
202  // Initialise sums
203  double ks = 0;
204  double ko = 0;
205  double bf = 0;
206  double sob = 0;
207  double sof = 0;
208  double ttl, ctl, ksli, koli, sobli, sofli, bfli;
209  double chi_s, chi_o, frho, ftau;
210  VectorType result(4);
211 
212  // Weighted sums over LIDF
213  for(unsigned int i=0; i<lidf.size(); ++i)
214  {
215  ttl = 2.5+5*i; // leaf inclination discrete values
216  ctl = vcl_cos(rd*ttl);
217  // SAIL volume scattering phase function gives interception and portions to be
218  // multiplied by rho and tau
219 
220  this->Volscatt(m_TTS, m_TTO, m_PSI, ttl, result);
221  chi_s = result[0];
222  chi_o = result[1];
223  frho = result[2];
224  ftau = result[3];
225 
226  //********************************************************************************
227  //* SUITS SYSTEM COEFFICIENTS
228  //*
229  //* ks : Extinction coefficient for direct solar flux
230  //* ko : Extinction coefficient for direct observed flux
231  //* att : Attenuation coefficient for diffuse flux
232  //* sigb : Backscattering coefficient of the diffuse downward flux
233  //* sigf : Forwardscattering coefficient of the diffuse upward flux
234  //* sf : Scattering coefficient of the direct solar flux for downward diffuse flux
235  //* sb : Scattering coefficient of the direct solar flux for upward diffuse flux
236  //* vf : Scattering coefficient of upward diffuse flux in the observed direction
237  //* vb : Scattering coefficient of downward diffuse flux in the observed direction
238  //* w : Bidirectional scattering coefficient
239  //********************************************************************************
240 
241  // Extinction coefficients
242  ksli = chi_s/cts;
243  koli = chi_o/cto;
244 
245  // Area scattering coefficient fractions
246  sobli = frho*CONST_PI/ctscto;
247  sofli = ftau*CONST_PI/ctscto;
248  bfli = ctl*ctl;
249  ks = ks+ksli*lidf[i];
250  ko = ko+koli*lidf[i];
251  bf = bf+bfli*lidf[i];
252  sob = sob+sobli*lidf[i];
253  sof = sof+sofli*lidf[i];
254  }
255 
256  // Geometric factors to be used later with rho and tau
257  double sdb, sdf, dob, dof, ddb, ddf;
258  sdb = 0.5*(ks+bf);
259  sdf = 0.5*(ks-bf);
260  dob = 0.5*(ko+bf);
261  dof = 0.5*(ko-bf);
262  ddb = 0.5*(1.+bf);
263  ddf = 0.5*(1.-bf);
264 
265  double lambda, Es, Ed, Rsoil1, Rsoil2, rsoil0, rho, tau, PARdiro, PARdifo;
266  double sigb, sigf, att, m2, m, sb, sf, vb, vf, w;
267  double tss, too, tsstoo, rdd, tdd, rsd, tsd, rdo, tdo, rsos, rsod;
268  double rddt, rsdt, rdot, rsodt, rsost, rsot, dn;
269  double e1, e2, rinf, rinf2, re, denom, J1ks, J2ks, J1ko, J2ko;
270  double Ps, Qs, Pv, Qv, z, g1, g2, Tv1, Tv2, T1, T2, T3;
271  double alf, sumint, fhot, x1, y1, f1, fint, x2, y2, f2;
272  double resh, resv;
273 
274  int nbdata = sizeof(DataSpecP5B) / sizeof(DataSpec);
275  for (int i = 0; i < nbdata; ++i)
276  {
277  lambda = DataSpecP5B[i].lambda;
278  Es = DataSpecP5B[i].directLight; //8
279  Ed = DataSpecP5B[i].diffuseLight; //9
280  Rsoil1 = DataSpecP5B[i].drySoil; //10
281  Rsoil2 = DataSpecP5B[i].wetSoil; //11
282  rho = inRefl->GetResponse()[i].second; //rho = LRT[1][i];
283  tau = inTrans->GetResponse()[i].second; //tau = LRT[2][i];
284 
285  // direct/diffuse light
286  //Es = direct
287  //Ed = diffuse
288  PARdiro = (1-m_Skyl/100.)*Es;
289  PARdifo = (m_Skyl/100.)*Ed;
290 
291  // Soil Reflectance Properties
292  //rsoil1 = dry soil
293  //rsoil2 = wet soil
294  rsoil0 = m_PSoil*Rsoil1+(1-m_PSoil)*Rsoil2;
295 
296  // Here rho and tau come in
297  sigb = ddb*rho+ddf*tau;
298  sigf = ddf*rho+ddb*tau;
299  att = 1.-sigf;
300  m2 = (att+sigb)*(att-sigb);
301  if(m2<=0) m2 = 0;
302  m = vcl_sqrt(m2);
303 
304 
305  sb = sdb*rho+sdf*tau;
306  sf = sdf*rho+sdb*tau;
307  vb = dob*rho+dof*tau;
308  vf = dof*rho+dob*tau;
309  w = sob*rho+sof*tau;
310 
311  // Here the LAI comes in
312  // Outputs for the case LAI = 0
313  if (m_LAI<0)
314  {
315  tss = 1;
316  too = 1;
317  tsstoo = 1;
318  rdd = 0;
319  tdd = 1;
320  rsd = 0;
321  tsd = 0;
322  rdo = 0;
323  tdo = 0;
324  //rso = 0;
325  rsos = 0;
326  rsod = 0;
327 
328  rddt = rsoil0;
329  rsdt = rsoil0;
330  rdot = rsoil0;
331  rsodt = 0;
332  rsost = rsoil0;
333  rsot = rsoil0;
334  }
335 
336  // Other cases (LAI > 0)
337  e1 = exp(-m*m_LAI);
338  e2 = e1*e1;
339  rinf = (att-m)/sigb;
340  rinf2 = rinf*rinf;
341  re = rinf*e1;
342  denom = 1.-rinf2*e2;
343 
344  J1ks=Jfunc1(ks, m, m_LAI);
345  J2ks=Jfunc2(ks, m, m_LAI);
346  J1ko=Jfunc1(ko, m, m_LAI);
347  J2ko=Jfunc2(ko, m, m_LAI);
348 
349  Ps = (sf+sb*rinf)*J1ks;
350  Qs = (sf*rinf+sb)*J2ks;
351  Pv = (vf+vb*rinf)*J1ko;
352  Qv = (vf*rinf+vb)*J2ko;
353 
354  rdd = rinf*(1.-e2)/denom;
355  tdd = (1.-rinf2)*e1/denom;
356  tsd = (Ps-re*Qs)/denom;
357  rsd = (Qs-re*Ps)/denom;
358  tdo = (Pv-re*Qv)/denom;
359  rdo = (Qv-re*Pv)/denom;
360 
361  tss = exp(-ks*m_LAI);
362  too = exp(-ko*m_LAI);
363  z = Jfunc3(ks, ko, m_LAI);
364  g1 = (z-J1ks*too)/(ko+m);
365  g2 = (z-J1ko*tss)/(ks+m);
366 
367  Tv1 = (vf*rinf+vb)*g1;
368  Tv2 = (vf+vb*rinf)*g2;
369  T1 = Tv1*(sf+sb*rinf);
370  T2 = Tv2*(sf*rinf+sb);
371  T3 = (rdo*Qs+tdo*Ps)*rinf;
372 
373  // Multiple scattering contribution to bidirectional canopy reflectance
374  rsod = (T1+T2-T3)/(1.-rinf2);
375 
376  // Treatment of the hotspot-effect
377  alf=1e6;
378  // Apply correction 2/(K+k) suggested by F.-M. Bron
379  if (m_HSpot>0) alf=(dso/m_HSpot)*2./(ks+ko);
380  if (alf>200) alf=200;
381  if (alf==0)
382  {
383  // The pure hotspot - no shadow
384  tsstoo = tss;
385  sumint = (1-tss)/(ks*m_LAI);
386  }
387  else
388  {
389  // Outside the hotspot
390  fhot=m_LAI*vcl_sqrt(ko*ks);
391  // Integrate by exponential Simpson method in 20 steps
392  // the steps are arranged according to equal partitioning
393  // of the slope of the joint probability function
394  x1=0;
395  y1=0;
396  f1=1;
397  fint=(1.-exp(-alf))*0.05;
398  sumint=0;
399 
400  for(unsigned int j=1; j<=20; ++j)
401  {
402  if (j<20) x2 = -vcl_log(1.-j*fint)/alf;
403  else x2 = 1;
404  y2 = -(ko+ks)*m_LAI*x2+fhot*(1.-exp(-alf*x2))/alf;
405  f2 = exp(y2);
406  sumint = sumint+(f2-f1)*(x2-x1)/(y2-y1);
407  x1=x2;
408  y1=y2;
409  f1=f2;
410  }
411  tsstoo=f1;
412  }
413 
414  // Bidirectional reflectance
415  // Single scattering contribution
416  rsos = w*m_LAI*sumint;
417  // Total canopy contribution
418  // rso=rsos+rsod;
419  //Interaction with the soil
420  dn=1.-rsoil0*rdd;
421 
422  rddt = rdd+tdd*rsoil0*tdd/dn;
423  rsdt = rsd+(tsd+tss)*rsoil0*tdd/dn;
424  rdot = rdo+tdd*rsoil0*(tdo+too)/dn;
425 
426  rsodt = rsod+((tss+tsd)*tdo+(tsd+tss*rsoil0*rdd)*too)*rsoil0/dn;
427  rsost = rsos+tsstoo*rsoil0;
428  rsot = rsost+rsodt;
429 
430  resh = (rddt*PARdifo+rsdt*PARdiro)/(PARdiro+PARdifo);
431  resv = (rdot*PARdifo+rsot*PARdiro)/(PARdiro+PARdifo);
432 
434  tmp1.first=lambda/1000.0;
435  tmp1.second=resh;
436  tmp2.first=lambda/1000.0;
437  tmp2.second=resv;
438 
439  outVRefl->GetResponse().push_back(tmp2);
440  outHRefl->GetResponse().push_back(tmp1);
441 
442  }
443 }
444 
445 
446 void
447 SailModel
448 ::Calc_LIDF(const double a, VectorType &lidf)
449 {
450  int ala=a;
451  VectorType freq;
452  Campbell(ala, freq);
453  lidf=freq;
454 
455 }
456 
457 
458 void
459 SailModel
460 ::Campbell(const double ala, VectorType &freq)
461 {
462  unsigned int n=18;
463  double excent = exp(-1.6184e-5*vcl_pow(ala, 3)+2.1145e-3*ala*ala-1.2390e-1*ala+3.2491);
464  double sum=0;
465  unsigned int tx2, tx1;
466  double tl1, tl2, x1, x2, v, alpha, alpha2, x12, x22, alpx1, alpx2, dum, almx1, almx2;
467  VectorType temp;
468 
469  for(unsigned int i=0; i<n; ++i)
470  {
471  tx2 = 5*i;
472  tx1 = 5*(i+1);
473  tl1 = tx1*CONST_PI/180;
474  tl2 = tx2*CONST_PI/180;
475 
476 
477  x1 = excent/sqrt(1.+excent*excent*vcl_tan(tl1)*vcl_tan(tl1));
478  x2 = excent/sqrt(1.+excent*excent*vcl_tan(tl2)*vcl_tan(tl2));
479  if(excent==1)
480  {
481  v = abs(cos(tl1)-cos(tl2));
482  temp.push_back( v );
483  sum = sum + v;
484  }
485  else
486  {
487  alpha = excent/vcl_sqrt(vcl_abs(1.-excent*excent));
488  alpha2 = alpha*alpha;
489  x12 = x1*x1;
490  x22 = x2*x2;
491  if(excent>1)
492  {
493  alpx1 = vcl_sqrt(alpha2+x12);
494  alpx2 = vcl_sqrt(alpha2+x22);
495  dum = x1*alpx1+alpha2*log(x1+alpx1);
496  v = vcl_abs(dum-(x2*alpx2+alpha2*log(x2+alpx2)));
497  temp.push_back( v );
498  sum = sum + v;
499  }
500  else
501  {
502  almx1 = sqrt(alpha2-x12);
503  almx2 = sqrt(alpha2-x22);
504  dum = x1*almx1+alpha2*asin(x1/alpha);
505  v = vcl_abs(dum-(x2*almx2+alpha2*asin(x2/alpha)));
506  temp.push_back( v );
507  sum = sum + v;
508  }
509  }
510  }
511 
512  for(unsigned int i=0; i<n; ++i)
513  {
514  freq.push_back(temp[i]/sum);
515  }
516 
517 }
518 
519 
520 void
521 SailModel
522 ::Volscatt(const double tts, const double tto, const double psi, const double ttl, VectorType &result)
523 {
524 
525  double rd = CONST_PI/180;
526  double costs = vcl_cos(rd*tts);
527  double costo = vcl_cos(rd*tto);
528  double sints = vcl_sin(rd*tts);
529  double sinto = vcl_sin(rd*tto);
530  double cospsi = vcl_cos(rd*psi);
531  double psir = rd*psi;
532  double costl = vcl_cos(rd*ttl);
533  double sintl = vcl_sin(rd*ttl);
534  double cs = costl*costs;
535  double co = costl*costo;
536  double ss = sintl*sints;
537  double so = sintl*sinto;
538 
539  // ..............................................................................
540  // betas -bts- and betao -bto- computation
541  // Transition angles (beta) for solar (betas) and view (betao) directions
542  // if thetav+thetal>pi/2, bottom side of the leaves is observed for leaf azimut
543  // interval betao+phi<leaf azimut<2pi-betao+phi.
544  // if thetav+thetal<pi/2, top side of the leaves is always observed, betao=pi
545  // same consideration for solar direction to compute betas
546  // ..............................................................................
547  double cosbts, cosbto, bts, ds, chi_s, bto, doo, chi_o;
548  double btran1, btran2, bt1, bt2, bt3, t1, t2 , denom, frho, ftau;
549 
550  cosbts = 5;
551  if (vcl_abs(ss)>1e-6) cosbts = -cs/ss;
552 
553  cosbto=5;
554  if (vcl_abs(so)>1e-6) cosbto = -co/so;
555 
556 
557  if (vcl_abs(cosbts)<1)
558  {
559  bts = vcl_acos(cosbts);
560  ds = ss;
561  }
562  else
563  {
564  bts = CONST_PI;
565  ds = cs;
566  }
567 
568  chi_s = 2./CONST_PI*((bts-CONST_PI*0.5)*cs+vcl_sin(bts)*ss);
569 
570  if (vcl_abs(cosbto)<1)
571  {
572  bto = vcl_acos(cosbto);
573  doo = so;
574  }
575  else if(tto<90)
576  {
577  bto = CONST_PI;
578  doo = co;
579  }
580  else
581  {
582  bto = 0;
583  doo = -co;
584  }
585  chi_o = 2./CONST_PI*((bto-CONST_PI*0.5)*co+vcl_sin(bto)*so);
586 
587  // ..............................................................................
588  // Computation of auxiliary azimut angles bt1, bt2, bt3 used
589  // for the computation of the bidirectional scattering coefficient w
590  // .............................................................................
591 
592  btran1 = vcl_abs(bts-bto);
593  btran2 = CONST_PI - vcl_abs(bts+bto-CONST_PI);
594 
595  if (psir<=btran1)
596  {
597  bt1=psir;
598  bt2=btran1;
599  bt3=btran2;
600  }
601  else
602  {
603  bt1=btran1;
604  if (psir<=btran2)
605  {
606  bt2=psir;
607  bt3=btran2;
608  }
609  else
610  {
611  bt2=btran2;
612  bt3=psir;
613  }
614  }
615 
616  t1 = 2.*cs*co+ss*so*cospsi;
617  t2 = 0;
618  if (bt2>0) t2=sin(bt2)*(2.*ds*doo+ss*so*cos(bt1)*cos(bt3));
619  denom = 2.*CONST_PI*CONST_PI;
620  frho = ((CONST_PI-bt2)*t1+t2)/denom;
621  ftau = (-bt2*t1+t2)/denom;
622 
623  if (frho<0) frho = 0;
624  if (ftau<0) ftau = 0;
625 
626  result[0] = chi_s;
627  result[1] = chi_o;
628  result[2] = frho;
629  result[3] = ftau;
630 
631 }
632 
633 
634 double
635 SailModel
636 ::Jfunc1(const double k, const double l, const double t)
637 {
638  //J1 function with avoidance of singularity problem
639  double v;
640  double del=(k-l)*t;
641  if(vcl_abs(del)>1e-3)
642  {
643  v = (exp(-l*t)-exp(-k*t))/(k-l);
644  return v;
645  }
646  else
647  {
648  v = 0.5*t*(exp(-k*t)+exp(-l*t))*(1.-del*del/12.);
649  return v;
650  }
651 }
652 
653 
654 double
655 SailModel
656 ::Jfunc2(const double k, const double l, const double t)
657 {
658  double v;
659  v = (1.-exp(-(k+l)*t))/(k+l);
660  return v;
661 }
662 
663 
664 double
665 SailModel
666 ::Jfunc3(const double k, const double l, const double t)
667 {
668  double v;
669  v = (1.-exp(-(k+l)*t))/(k+l);
670  return v;
671 }
672 
673 
674 void
675 SailModel
676 ::PrintSelf(std::ostream& os, itk::Indent indent) const
677 {
678  Superclass::PrintSelf(os, indent);
679 
680 }
681 } // end namespace otb
682 
const double CONST_PI
Definition: otbMath.h:45
std::vector< double > VectorType
Definition: otbSailModel.h:50
void Campbell(const double ala, VectorType &freq)
const DataSpec DataSpecP5B[]
std::pair< TPrecision, TValuePrecision > PairType
void SetTransmittance(const SpectralResponseType *object)
ObjectType * GetPointer() const
double Jfunc3(const double k, const double l, const double t)
double directLight
virtual SpectralResponseType * GetHemisphericalReflectance()
virtual ~SailModel()
virtual void GenerateData()
void SetInput(const ParametersType &)
virtual SpectralResponseType * GetViewingReflectance()
void Volscatt(const double tts, const double tto, const double psi, const double ttl, VectorType &result)
const ParametersType GetInput()
void PrintSelf(std::ostream &os, itk::Indent indent) const
void Calc_LIDF(const double a, VectorType &lidf)
void SetReflectance(const SpectralResponseType *object)
double Jfunc1(const double k, const double l, const double t)
SpectralResponseType * GetReflectance()
This class represents the spectral response of an object (or a satellite band).
double diffuseLight
SizeValueType Size(void) const
This struct contains the specific absorption coefficients of each constituent and the refractive inde...
double Jfunc2(const double k, const double l, const double t)
virtual DataObjectPointer MakeOutput(unsigned int)
SpectralResponseType * GetTransmittance()