Orfeo Toolbox  4.0
otbProspectModel.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 
19 #include "itkNumericTraits.h"
20 
21 #include "otbProspectModel.h"
22 #include <boost/math/special_functions/expint.hpp>
23 #include <boost/shared_ptr.hpp>
24 #include "otbMath.h"
25 
26 //TODO check EPSILON matlab
27 #define EPSILON 0.0000000000000000000000001
28 
29 namespace otb
30 {
31 
35 {
36  this->ProcessObject::SetNumberOfRequiredInputs(1);
37  this->ProcessObject::SetNumberOfRequiredOutputs(2);
38 
39  SpectralResponseType::Pointer outputRefl = static_cast<SpectralResponseType *>(this->MakeOutput(0).GetPointer());
40  this->itk::ProcessObject::SetNthOutput(0, outputRefl.GetPointer());
41 
42  SpectralResponseType::Pointer outputTrans = static_cast<SpectralResponseType *>(this->MakeOutput(1).GetPointer());
43  this->itk::ProcessObject::SetNthOutput(1, outputTrans.GetPointer());
44 }
45 
49 {}
50 
52 void
55 {
56  this->itk::ProcessObject::SetNthInput(0, const_cast<LeafParametersType *>(object));
57 }
58 
63 {
64  if(this->GetNumberOfInputs() != 1)
65  {
66  //exit
67  return 0;
68  }
69  return static_cast<LeafParametersType *>(this->itk::ProcessObject::GetInput(0));
70 }
71 
75 ::MakeOutput(unsigned int)
76 {
77  return static_cast<itk::DataObject *>(SpectralResponseType::New().GetPointer());
78 }
79 
84 {
85  if(this->GetNumberOfOutputs() < 2)
86  {
87  //exit
88  return 0;
89  }
90  return static_cast<SpectralResponseType *>(this->itk::ProcessObject::GetOutput(0));
91 }
92 
97 {
98  if(this->GetNumberOfOutputs() < 2)
99  {
100  //exit
101  return 0;
102  }
103  return static_cast<SpectralResponseType *>(this->itk::ProcessObject::GetOutput(1));
104 }
105 
106 
108 void
110 ::SetInput(const ParametersType & params)
111 {
112 // m_Parameters = params;
113  if(params.Size()!=6) itkExceptionMacro( << "Must have 6 parameters in that order : Cab, Car, CBrown, Cw, Cm, N" );
114  LeafParametersType::Pointer leafParams = LeafParametersType::New();
115  leafParams->SetCab(params[0]);
116  leafParams->SetCar(params[1]);
117  leafParams->SetCBrown(params[2]);
118  leafParams->SetCw(params[3]);
119  leafParams->SetCm(params[4]);
120  leafParams->SetN(params[5]);
121 
122  this->itk::ProcessObject::SetNthInput(0, leafParams);
123 }
124 
125 
127 void
130 {
131 
132  LeafParametersType::Pointer leafParameters = this->GetInput();
133  SpectralResponseType::Pointer outRefl = this->GetReflectance();
134  SpectralResponseType::Pointer outTrans = this->GetTransmittance();
135 
136  unsigned int alpha=40;
137  double lambda, n, k, trans, t12, temp, t21, r12, r21, x, y, ra, ta, r90, t90;
138  double delta, beta, va, vb, vbNN, vbNNinv, vainv, s1, s2, s3, RN, TN;
139  double N, Cab, Car, CBrown, Cw, Cm;
140 
141  N = leafParameters->GetN();
142  Cab = leafParameters->GetCab();
143  Car = leafParameters->GetCar();
144  CBrown = leafParameters->GetCBrown();
145  Cw = leafParameters->GetCw();
146  Cm = leafParameters->GetCm();
147 
148  int nbdata = sizeof(DataSpecP5B) / sizeof(DataSpec);
149  for (int i = 0; i < nbdata; ++i)
150  {
151  lambda = DataSpecP5B[i].lambda;
152  n = DataSpecP5B[i].refLeafMatInd;
153 
155  k = k + Cm*DataSpecP5B[i].dryAbsCoef;
156  k = k / N;
158 
159  trans=(1.-k)*exp(-k)+k*k*boost::math::expint(1, k);
160 
161  t12 = this->Tav(alpha, n);
162  temp = this->Tav(90, n);
163 
164 
165  t21 = temp/(n*n);
166  r12 = 1.-t12;
167  r21 = 1.-t21;
168  x = t12/temp;
169  y = x*(temp-1)+1-t12;
170 
171  ra = r12+(t12*t21*r21*(trans*trans))/(1.-r21*r21*trans*trans);
172  ta = (t12*t21*trans)/(1.-r21*r21*trans*trans);
173  r90 = (ra-y)/x;
174  t90 = ta/x;
175 
176  delta = (t90*t90-r90*r90-1.)*(t90*t90-r90*r90-1.) - 4.*r90*r90;
177  if(delta < 0) delta = EPSILON;
178  else delta=vcl_sqrt(delta);
179 
180  beta = (1.+r90*r90-t90*t90-delta)/(2.*r90);
181  va=(1.+r90*r90-t90*t90+delta)/(2.*r90);
182  if ((beta-r90)<=0)
183  vb=vcl_sqrt(beta*(va-r90)/(va*EPSILON));
184  else
185  vb=vcl_sqrt(beta*(va-r90)/(va*(beta-r90)));
186 
187  vbNN = vcl_pow(vb, N-1.);
188  vbNNinv = 1./vbNN;
189  vainv = 1./va;
190  s1=ta*t90*(vbNN-vbNNinv);
191  s2=ta*(va-vainv);
192  s3=va*vbNN-vainv*vbNNinv-r90*(vbNN-vbNNinv);
193 
194  RN=ra+s1/s3;
195  TN=s2/s3;
196 
197 
200  refl.first=lambda/1000.0;
201  refl.second=RN;
202  trans.first=lambda/1000.0;
203  trans.second=TN;
204  outRefl->GetResponse().push_back(refl);
205  outTrans->GetResponse().push_back(trans);
206  }
207 }
208 
209 
210 double
212 ::Tav(const int theta, double ref)
213 {
214 
215  double theta_rad = theta*CONST_PI/180;
216  double r2, rp, rm, a, k, ds, k2, rm2, res, b1, b2, b;
217  double ts, tp1, tp2, tp3, tp4, tp5, tp;
218 
219  r2=ref*ref;
220  rp=r2+1;
221  rm=r2-1;
222  a=(ref+1)*(ref+1)/2;
223  k=-(r2-1)*(r2-1)/4;
224  ds=sin(theta_rad);
225 
226  k2=k*k;
227  rm2=rm*rm;
228 
229  if(theta_rad==0) res=4*ref/((ref+1)*(ref+1));
230  else
231  {
232  if(theta_rad==CONST_PI/2) b1=itk::NumericTraits<double>::ZeroValue();
233  else b1=vcl_sqrt((ds*ds-rp/2)*(ds*ds-rp/2)+k);
234 
235  b2=ds*ds-rp/2;
236  b=b1-b2;
237  ts=(k2/(6*vcl_pow(b, 3))+k/b-b/2)-(k2/(6*vcl_pow(a, 3))+k/a-a/2);
238  tp1=-2*r2*(b-a)/(rp*rp);
239  tp2=-2*r2*rp*log(b/a)/rm2;
240  tp3=r2*(1./b-1./a)/2;
241  tp4=16*r2*r2*(r2*r2+1)*log((2*rp*b-rm2)/(2*rp*a-rm2))/(vcl_pow(rp, 3)*rm2);
242  tp5=16*vcl_pow(r2, 3)*(1./(2*rp*b-rm2)-1./(2*rp*a-rm2))/vcl_pow(rp, 3);
243  tp=tp1+tp2+tp3+tp4+tp5;
244  res=(ts+tp)/(2*ds*ds);
245  }
246  return res;
247 
248 
249 }
250 
251 void
253 ::PrintSelf(std::ostream& os, itk::Indent indent) const
254 {
255  Superclass::PrintSelf(os, indent);
256 
257 }
258 } // end namespace otb
259 

Generated at Sat Mar 8 2014 16:14:19 for Orfeo Toolbox with doxygen 1.8.3.1