OTB  6.7.0
Orfeo Toolbox
otbSpectralResponse.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 otbSpectralResponse_hxx
22 #define otbSpectralResponse_hxx
23 
24 #include "itkNumericTraits.h"
25 
26 #include "otbSpectralResponse.h"
27 
28 #include <algorithm>
29 #include <fstream>
30 
31 namespace otb
32 {
33 
34 template<class TPrecision, class TValuePrecision>
36 {
37  m_SensitivityThreshold = 0.01;
38  m_IntervalComputed = false;
39  m_PosGuess = 0;
40  m_UsePosGuess=false;
41 }
42 
43 template<class TPrecision, class TValuePrecision>
44 void SpectralResponse<TPrecision, TValuePrecision>::Load(const std::string & filename,
45  ValuePrecisionType coefNormalization)
46 {
47  //Parse JPL file spectral response (ASCII file)
48  //Begin line 27
49  std::ifstream fin(filename);
50  if (fin.fail())
51  {
52  itkExceptionMacro(<<"Error opening file" << filename);
53  }
54 
55  int NumLigne = 26; // Go to the line 27
56  //Ignore first 26th lines which are metadatas information
57  for (int i = 0; i < NumLigne; ++i)
58  fin.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
59 
60  while (!fin.eof())
61  {
62  //For each
63  std::pair<TPrecision, TValuePrecision> currentPair;
64 
65  fin >> currentPair.first;
66  fin >> currentPair.second;
67  currentPair.second = currentPair.second / coefNormalization;
68  if (currentPair.first != itk::NumericTraits<TPrecision>::ZeroValue() && currentPair.second != itk::NumericTraits<
69  TValuePrecision>::ZeroValue())
70  //Add not null pair of values to the vector
71  m_Response.push_back(currentPair);
72  }
73  fin.close();
74  //Sort the vector using the specific functor sort_pair
75  std::sort(m_Response.begin(), m_Response.end(), sort_pair());
76 
77  m_IntervalComputed = false;
78 }
79 
80 template<class TPrecision, class TValuePrecision>
82 {
83  m_Response.clear();
84  m_IntervalComputed = false;
85  return true;
86 }
87 
88 template<class TPrecision, class TValuePrecision>
90 {
91  return m_Response.size();
92 }
93 
94 
95 template<class TPrecision, class TValuePrecision>
97 {
98  m_PosGuess = 0;
99  if (m_Response.size() <= 1)
100  {
101  itkExceptionMacro(<<"ERROR spectral response need at least 2 value to perform interpolation.");
102  }
103 
104  TPrecision lambdaMax = this->GetInterval().second;
105  if (lambda > lambdaMax) return;
106  typename VectorPairType::const_iterator it = m_Response.begin();
107 
108  while (((*it).first < lambda))
109  {
110  m_PosGuess++;
111  ++it;
112  if (it == (m_Response.end())) return;
113  }
114 
115  if (m_PosGuess > 0) m_PosGuess--;
116  return;
117 }
118 
119 
120 template<class TPrecision, class TValuePrecision>
122  TValuePrecision>::operator()(const PrecisionType & lambda)
123 {
124 
125  //Suppose that the vector is sorted
126 
127  //Guess a starting lambda
128  if (m_Response.size() <= 1)
129  {
130  itkExceptionMacro(<<"ERROR spectral response need at least 2 value to perform interpolation.");
131  }
132 
133  typename VectorPairType::const_iterator beg = m_Response.begin();
134  typename VectorPairType::const_iterator last = m_Response.end();
135  --last;
136 
137  TPrecision lambdaMin = this->GetInterval().first;
138  TPrecision lambdaMax = this->GetInterval().second;
139 
140  if (lambda < lambdaMin) return static_cast<TValuePrecision> (0.0);
141  if (lambda > lambdaMax) return static_cast<TValuePrecision> (0.0);
142 
143  typename VectorPairType::const_iterator it;
144 
145  if(m_UsePosGuess)
146  it= beg + m_PosGuess;
147  else
148  it= beg;
149 
150  TPrecision lambda1 = (*beg).first;
151  TValuePrecision SR1 = static_cast<TValuePrecision> (0.0);
152  while (((*it).first < lambda))
153  {
154 
155  lambda1 = (*it).first;
156  SR1 = (*it).second;
157  ++it;
158  if (it == (m_Response.end()))
159  {
160  return static_cast<TValuePrecision> (0.0);
161  }
162  }
163  TPrecision lambda2 = (*it).first;
164 
165  // if the guess is just right
166  if (lambda2 == lambda)
167  {
168  return (*it).second;
169  }
170  else
171  {
172 
173  TPrecision lambdaDist = lambda - lambda1;
174  TPrecision ratio = lambdaDist / (lambda2 - lambda1);
175 
176  TValuePrecision SR2 = (*it).second;
177 
178  return static_cast<TValuePrecision> (ratio * SR1 + (1 - ratio) * SR2);
179 
180  }
181 
182 }
183 
184 template<class TPrecision, class TValuePrecision>
186  ImagePointerType image) const
187 {
188  typename ImageType::IndexType start;
189  start[0] = 0;
190  start[1] = 0;
191 
192  typename ImageType::SizeType size;
193  size[0] = 1;
194  size[1] = 1;
195 
196  typename ImageType::PointType origin;
197  origin[0] = 0;
198  origin[1] = 0;
199 
200  typename ImageType::SpacingType spacing;
201  spacing[0] = 1;
202  spacing[1] = 1;
203 
204  typename ImageType::RegionType region;
205  region.SetSize(size);
206  region.SetIndex(start);
207 
208  image->SetRegions(region);
209  image->SetNumberOfComponentsPerPixel(this->Size());
210  image->Allocate();
211 
212  typename ImageType::IndexType idx;
213  typename ImageType::PixelType pixel;
214  pixel.SetSize(this->Size());
215 
216  for (unsigned int j = 0; j < this->Size(); ++j)
217  {
218  pixel[j] = m_Response[j].second;
219  }
220  idx[0] = 0;
221  idx[1] = 0;
222  image->SetPixel(idx, pixel);
223  return image;
224 }
225 
226 template<class TPrecision, class TValuePrecision>
228 {
229 
230  typename ImageType::IndexType idx;
231  idx[0] = 0;
232  idx[1] = 0;
233 
234  for (unsigned int j = 0; j < this->Size(); ++j)
235  {
236  m_Response[j].second = image->GetPixel(idx)[j];
237  }
238  m_IntervalComputed = false;
239 }
240 
241 template<class TPrecision, class TValuePrecision>
243  TValuePrecision>::GetFilterFunctionValues(double step)
244 {
245 
246  //Assume that the SR is sorted
247  typename FilterFunctionValuesType::ValuesVectorType valuesVector;
248  Self& responseCalculator = *this;
249  for (double i = m_Response.front().first; i <= m_Response.back().first; i += step)
250  {
251  valuesVector.push_back(responseCalculator(i));
252  }
253  FilterFunctionValuesPointerType functionValues = FilterFunctionValuesType::New();
254 
255  functionValues->SetFilterFunctionValues(valuesVector);
256  functionValues->SetMinSpectralValue(m_Response.front().first);
257  functionValues->SetMaxSpectralValue(m_Response.back().first);
258  functionValues->SetUserStep(step);
259 
260  return functionValues;
261 }
262 
263 template<class TPrecision, class TValuePrecision>
265 {
266  typename VectorPairType::const_iterator it = m_Response.begin();
267 
268  while ((*it).second <= m_SensitivityThreshold)
269  {
270  ++it;
271  if (it == (m_Response.end() - 1))
272  {
273  m_Interval.first = static_cast<TPrecision> (0.0);
274  m_Interval.second = static_cast<TPrecision> (0.0);
275  m_IntervalComputed = true;
276  return;
277  }
278 
279  }
280  m_Interval.first = (*it).first;
281  it = (m_Response.end() - 1);
282  while ((*it).second <= m_SensitivityThreshold)
283  {
284  if (it == (m_Response.begin()))
285  {
286  m_Interval.second = (*it).first;
287  m_IntervalComputed = true;
288  return;
289  }
290  --it;
291 
292  }
293 
294  m_Interval.second = (*it).first;
295  m_IntervalComputed = true;
296 }
297 
298 template<class TPrecision, class TValuePrecision>
300 {
301  Superclass::PrintSelf(os, indent);
302  os << std::endl;
303  os << indent << "[Wavelength (micrometers), Reflectance (percent)]" << std::endl;
304  for (typename VectorPairType::const_iterator it = m_Response.begin(); it != m_Response.end(); ++it)
305  {
306  os << indent << "Num " << it - m_Response.begin() << ": [" << (*it).first << "," << (*it).second << "]"
307  << std::endl;
308  }
309 }
310 
311 } // end namespace otb
312 
313 #endif
Superclass::SpacingType SpacingType
Superclass::IndexType IndexType
ImagePointerType GetImage(ImagePointerType image) const
virtual unsigned int Size() const
void SetPosGuessMin(const PrecisionType &lambda)
void SetFromImage(ImagePointerType image)
std::vector< WavelengthSpectralBandType > ValuesVectorType
Superclass::SizeType SizeType
Superclass::PointType PointType
void Load(const std::string &filename, ValuePrecisionType coefNormalization=1.0)
void PrintSelf(std::ostream &os, itk::Indent indent) const override
TValuePrecision ValuePrecisionType
This class represents the spectral response of an object (or a satellite band).
void SetSize(unsigned int sz, TReallocatePolicy reallocatePolicy, TKeepValuesPolicy keepValues)
Superclass::RegionType RegionType