Orfeo Toolbox  3.16
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 
32 SailModel
33 ::SailModel()
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 
56 SailModel
57 ::~SailModel()
58 {}
59 
61 void
62 SailModel
63 ::SetReflectance(const SpectralResponseType * object)
64 {
65  this->itk::ProcessObject::SetNthInput(0, const_cast<SpectralResponseType *>(object));
66 }
67 
68 SailModel::SpectralResponseType *
69 SailModel
70 ::GetReflectance()
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
82 SailModel
83 ::SetTransmittance(const SpectralResponseType * object)
84 {
85  this->itk::ProcessObject::SetNthInput(1, const_cast<SpectralResponseType *>(object));
86 }
87 
88 SailModel::SpectralResponseType *
89 SailModel
90 ::GetTransmittance()
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 
109 SailModel::SpectralResponseType *
110 SailModel
111 ::GetViewingReflectance()
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 
122 SailModel::SpectralResponseType *
123 SailModel
124 ::GetHemisphericalReflectance()
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 
154 const SailModel::ParametersType
155 SailModel
156 ::GetInput()
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
178 ::GenerateData()
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 
433  SpectralResponseType::PairType tmp1, tmp2;
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 
Generated at Sun Feb 3 2013 00:45:16 for Orfeo Toolbox with doxygen 1.8.1.1