OTB  9.0.0
Orfeo Toolbox
otbTimeSeriesLeastSquareFittingFunctor.h
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2005-2022 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 otbTimeSeriesLSFF_h
22 #define otbTimeSeriesLSFF_h
23 
24 #include "vnl/algo/vnl_matrix_inverse.h"
25 #include "vnl/vnl_transpose.h"
26 #include "vnl/vnl_matrix.h"
27 
28 namespace otb
29 {
30 namespace Functor
31 {
58 template <class TSeriesType, class TTimeFunction, class TDateType = TSeriesType, class TWeightType = TSeriesType>
60 {
61 public:
62  typedef typename TTimeFunction::CoefficientsType CoefficientsType;
63 
66  {
67  for (unsigned int i = 0; i < m_WeightSeries.Size(); ++i)
68  m_WeightSeries[i] = 1.0;
69  }
72  {
73  }
74 
75  inline TSeriesType operator()(const TSeriesType& series)
76  {
77  TTimeFunction estFunction = this->EstimateTimeFunction(series);
78  TSeriesType outSeries;
79  for (unsigned int i = 0; i < m_DoySeries.Size(); ++i)
80  outSeries[i] = estFunction.GetValue(m_DoySeries[i]);
81  return outSeries;
82  }
83 
84  inline void SetDates(const TDateType& doy)
85  {
86  for (unsigned int i = 0; i < doy.Size(); ++i)
87  m_DoySeries[i] = doy[i];
88  }
89 
90  inline void SetWeights(const TWeightType& weights)
91  {
92  for (unsigned int i = 0; i < weights.Size(); ++i)
93  m_WeightSeries[i] = weights[i];
94  }
95 
96  inline CoefficientsType GetCoefficients(const TSeriesType& series) const
97  {
98  return (this->EstimateTimeFunction(series)).GetCoefficients();
99  }
100 
101  inline TTimeFunction EstimateTimeFunction(const TSeriesType& series) const
102  {
103  TTimeFunction estFunction;
104  unsigned int nbDates = m_DoySeries.Size();
105  unsigned int nbCoefs = estFunction.GetCoefficients().Size();
106 
107  // b = A * c
108  vnl_matrix<double> A(nbDates, nbCoefs);
109  vnl_matrix<double> b(nbDates, 1);
110 
111  // fill the matrices
112 
113  typename TTimeFunction::CoefficientsType tmpCoefs;
114  for (unsigned int j = 0; j < nbCoefs; ++j)
115  tmpCoefs[j] = 0.0;
116 
117  for (unsigned int i = 0; i < nbDates; ++i)
118  {
119  b.put(i, 0, series[i] / m_WeightSeries[i]);
120  for (unsigned int j = 0; j < nbCoefs; ++j)
121  {
122  tmpCoefs[j] = 1.0;
123  estFunction.SetCoefficients(tmpCoefs);
124  A.put(i, j, estFunction.GetValue(m_DoySeries[i]) / m_WeightSeries[i]);
125  tmpCoefs[j] = 0.0;
126  }
127  }
128  // solve the problem c = (At * A)^-1*At*b
129 
130  vnl_matrix<double> atainv = vnl_matrix_inverse<double>(vnl_transpose(A) * A);
131 
132  vnl_matrix<double> atainvat = atainv * vnl_transpose(A);
133  vnl_matrix<double> c = atainvat * b;
134 
135  for (unsigned int j = 0; j < nbCoefs; ++j)
136  tmpCoefs[j] = c.get(j, 0);
137  estFunction.SetCoefficients(tmpCoefs);
138 
139  return estFunction;
140  }
141 
142 private:
144  TDateType m_DoySeries;
145  TWeightType m_WeightSeries;
146 };
147 }
148 } // namespace otb
149 #endif
otb::Functor::TimeSeriesLeastSquareFittingFunctor::m_WeightSeries
TWeightType m_WeightSeries
Definition: otbTimeSeriesLeastSquareFittingFunctor.h:145
otb
The "otb" namespace contains all Orfeo Toolbox (OTB) classes.
Definition: otbJoinContainer.h:32
otb::Functor::TimeSeriesLeastSquareFittingFunctor::CoefficientsType
TTimeFunction::CoefficientsType CoefficientsType
Definition: otbTimeSeriesLeastSquareFittingFunctor.h:62
otb::Functor::TimeSeriesLeastSquareFittingFunctor::EstimateTimeFunction
TTimeFunction EstimateTimeFunction(const TSeriesType &series) const
Definition: otbTimeSeriesLeastSquareFittingFunctor.h:101
otb::Functor::TimeSeriesLeastSquareFittingFunctor::operator()
TSeriesType operator()(const TSeriesType &series)
Definition: otbTimeSeriesLeastSquareFittingFunctor.h:75
GapFilling::doy
unsigned int doy(const std::tm &d)
Return the day of year.
otb::Functor::TimeSeriesLeastSquareFittingFunctor
Implements a least squares fitting of a time profile.
Definition: otbTimeSeriesLeastSquareFittingFunctor.h:59
otb::Functor::TimeSeriesLeastSquareFittingFunctor::GetCoefficients
CoefficientsType GetCoefficients(const TSeriesType &series) const
Definition: otbTimeSeriesLeastSquareFittingFunctor.h:96
otb::Functor::TimeSeriesLeastSquareFittingFunctor::SetDates
void SetDates(const TDateType &doy)
Definition: otbTimeSeriesLeastSquareFittingFunctor.h:84
otb::Functor::TimeSeriesLeastSquareFittingFunctor::TimeSeriesLeastSquareFittingFunctor
TimeSeriesLeastSquareFittingFunctor()
Constructor.
Definition: otbTimeSeriesLeastSquareFittingFunctor.h:65
otb::Functor::TimeSeriesLeastSquareFittingFunctor::SetWeights
void SetWeights(const TWeightType &weights)
Definition: otbTimeSeriesLeastSquareFittingFunctor.h:90
otb::Functor::TimeSeriesLeastSquareFittingFunctor::~TimeSeriesLeastSquareFittingFunctor
virtual ~TimeSeriesLeastSquareFittingFunctor()
Destructor.
Definition: otbTimeSeriesLeastSquareFittingFunctor.h:71
otb::Functor::TimeSeriesLeastSquareFittingFunctor::m_DoySeries
TDateType m_DoySeries
Definition: otbTimeSeriesLeastSquareFittingFunctor.h:144