OTB  6.7.0
Orfeo Toolbox
otbReduceSpectralResponse.hxx
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2005-2019 Centre National d'Etudes Spatiales (CNES)
3  *
4  * This file is part of Orfeo Toolbox
5  *
6  * https://www.orfeo-toolbox.org/
7  *
8  * Licensed under the Apache License, Version 2.0 (the "License");
9  * you may not use this file except in compliance with the License.
10  * You may obtain a copy of the License at
11  *
12  * http://www.apache.org/licenses/LICENSE-2.0
13  *
14  * Unless required by applicable law or agreed to in writing, software
15  * distributed under the License is distributed on an "AS IS" BASIS,
16  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
17  * See the License for the specific language governing permissions and
18  * limitations under the License.
19  */
20 
21 #ifndef otbReduceSpectralResponse_hxx
22 #define otbReduceSpectralResponse_hxx
23 
24 #include <algorithm>
25 #include "itkNumericTraits.h"
27 
28 namespace otb
29 {
30 
31 template <class TSpectralResponse , class TRSR>
32 ReduceSpectralResponse <TSpectralResponse , TRSR>
33 ::ReduceSpectralResponse() : m_ReflectanceMode(false)
34 {
35  m_ReduceResponse = InputSpectralResponseType::New();
36 }
37 
38 template <class TSpectralResponse , class TRSR>
39 bool
42 {
43  return ( m_InputSatRSR->Clear() & m_InputSpectralResponse->Clear() );
44 }
45 
46 template <typename T>
47 inline
48 T trapezoid_area(T x1, T x2, T y1, T y2)
49 {
50  /* We compute the area of the trapezoid
51  by computing the lower square and the
52  upper triangle.
53  /- |y2
54  /- |
55  /- |
56  /- |
57  y1 +----------+
58  | |
59  | |
60  | |
61  +----------+
62  x1 x2
63 
64  (x2-x1)*min(y1,y2) + (x2-x1)*abs(y2-y1)/2
65  */
66  return (x2-x1)*std::min(y1,y2) + (x2-x1)*fabs(y2-y1)*0.5;
67 }
68 
69 template <class TSpectralResponse , class TRSR>
70 inline typename ReduceSpectralResponse<TSpectralResponse , TRSR>::ValuePrecisionType
72 ::operator()(const unsigned int numBand)
73 {
74  if (numBand >= m_InputSatRSR->GetNbBands())
75  {
76  itkExceptionMacro(<< "There is no band num " << numBand
77  << " in the RSR vector!(Size of the current RSR vector is "
78  << m_InputSatRSR->GetNbBands() << ")" );
79  }
80  else
81  {
83  typename InputRSRType::SpectralResponseType* solarIrradiance = this->m_InputSatRSR->GetSolarIrradiance();
84 
85  if(m_ReflectanceMode)
86  {
87  if (solarIrradiance == nullptr)
88  {
89  itkExceptionMacro(<<"Error occurs getting solar irradiance. Solar irradiance is mandatory using the reflectance mode.");
90  }
91  }
92  typename VectorPairType::const_iterator it;
93  VectorPairType pairs = (m_InputSatRSR->GetRSR())[numBand]->GetResponse();
94  it = pairs.begin();
95  ValuePrecisionType totalArea(0);
96  //start with the second value for the numerical integration
97  ++it;
98  while (it != pairs.end())
99  {
100  ValuePrecisionType rsr1 = (*(it-1)).second;
101  ValuePrecisionType rsr2 = (*it).second;
102  if(rsr1 > 0 || rsr2 >0)
103  {
104  PrecisionType lambda1 = (*(it-1)).first;
105  PrecisionType lambda2 = (*it).first;
106 // PrecisionType deltaLambda = lambda2-lambda1;
107  ValuePrecisionType spectrum1 = (*m_InputSpectralResponse)(lambda1);
108  ValuePrecisionType spectrum2 = (*m_InputSpectralResponse)(lambda2);
109  /*
110  In order to simplify the computation for the reflectance mode,
111  we introduce the solar irradiance in the general formula with
112  a value of 1.0 for the radiance case.
113 
114  In this way the formula is the same if we weight the RSR by
115  the solar irradiance before the integration.
116  */
117  ValuePrecisionType solarIrradiance1(1.0);
118  ValuePrecisionType solarIrradiance2(1.0);
119  if(m_ReflectanceMode)
120  {
121  solarIrradiance1 = (*solarIrradiance)(lambda1);
122  solarIrradiance2 = (*solarIrradiance)(lambda2);
123  }
124  rsr1 *= solarIrradiance1;
125  rsr2 *= solarIrradiance2;
126  res += trapezoid_area(lambda1, lambda2,
127  rsr1*spectrum1,
128  rsr2*spectrum2);
129  totalArea += trapezoid_area(lambda1, lambda2,
130  rsr1,
131  rsr2);
132  }
133  ++it;
134  }
135  return res / totalArea;
136  }
137 }
138 
139 template <class TSpectralResponse , class TRSR>
140 void
143 {
144  m_ReduceResponse->Clear();
145  //Compute the reduce response for each band of the sensor
146  for (unsigned int i=0; i<m_InputSatRSR->GetNbBands(); ++i)
147  {
148  m_InputSpectralResponse->SetPosGuessMin((this->m_InputSatRSR->GetRSR())[i]->GetInterval().first);
149  m_InputSpectralResponse->SetUsePosGuess(true);
150  PairType pair;
151  //pair.first = center wavelength of the band
152  pair.first = ((this->m_InputSatRSR->GetRSR())[i]->GetInterval().first
153  + (this->m_InputSatRSR->GetRSR())[i]->GetInterval().second);
154  pair.first = pair.first / 2.0;
155  pair.second = (*this)(i);
156  m_ReduceResponse->GetResponse().push_back(pair);
157  m_InputSpectralResponse->SetUsePosGuess(false);
158  }
159 
160 }
161 
162 
163 template <class TSpectralResponse , class TRSR>
164 void
166 ::LoadInputsFromFiles (const std::string & spectralResponseFile, const std::string & RSRFile, const unsigned int nbRSRBands, ValuePrecisionType coefNormSpectre, ValuePrecisionType coefNormRSR)
167 {
168  //Instantiation
169  m_InputSpectralResponse=InputSpectralResponseType::New();
171  m_InputSpectralResponse->Load(spectralResponseFile, coefNormSpectre);
172 
173  m_InputSatRSR=InputRSRType::New();
175  m_InputSatRSR->SetNbBands(nbRSRBands);
176 
178  m_InputSatRSR->Load(RSRFile, coefNormRSR);
179 }
180 
181 
182 template <class TSpectralResponse , class TRSR>
183 void
185 ::PrintSelf(std::ostream& os, itk::Indent indent) const
186 {
187  Superclass::PrintSelf(os, indent);
188  os<<std::endl;
189  os << "spectre " << m_InputSpectralResponse <<std::endl;
190  os << "Sat RSR " << m_InputSatRSR <<std::endl;
191  os<<std::endl;
192 
193  if(m_ReflectanceMode)
194  {
195  os << "Solar irradiance " << std::endl;
196  this->m_InputSatRSR->GetSolarIrradiance()->PrintSelf(os, indent);
197  os<<std::endl;
198  os <<indent << "[Center Wavelength (micrometers), Reflectance (percent)]" << std::endl;
199  }
200  else{
201  os <<indent << "[Center Wavelength (micrometers), Radiance (percent)]" << std::endl;
202  }
203 
204  for(typename VectorPairType::const_iterator it = m_ReduceResponse->GetResponse().begin(); it != m_ReduceResponse->GetResponse().end(); ++it)
205  {
206  os <<indent << "Band Nb : "<< it - m_ReduceResponse->GetResponse().begin() << ": [" << (*it).first << ","<< (*it).second << "]" << std::endl;
207  }
208 }
209 
210 
211 } // end namespace otb
212 
213 #endif
void PrintSelf(std::ostream &os, itk::Indent indent) const override
InputSpectralResponseType::VectorPairType VectorPairType
static T ZeroValue()
T trapezoid_area(T x1, T x2, T y1, T y2)
ValuePrecisionType operator()(const unsigned int numBand)
InputRSRType::PrecisionType PrecisionType
InputSpectralResponsePointerType m_ReduceResponse
This class computes the reduced spectral response of each band of a sensor.
InputRSRType::ValuePrecisionType ValuePrecisionType
void LoadInputsFromFiles(const std::string &spectralResponseFile, const std::string &RSRFile, const unsigned int nbRSRBands, ValuePrecisionType coefNormSpectre=1.0, ValuePrecisionType coefNormRSR=1.0)
InputSpectralResponseType::PairType PairType