/usr/include/OTB-5.8/otbTimeSeriesLeastSquareFittingFunctor.h is in libotb-dev 5.8.0+dfsg-3.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 | /*=========================================================================
Program: ORFEO Toolbox
Language: C++
Date: $Date$
Version: $Revision$
Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
See OTBCopyright.txt for details.
This software is distributed WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE. See the above copyright notices for more information.
=========================================================================*/
#ifndef otbTimeSeriesLSFF_h
#define otbTimeSeriesLSFF_h
#include "vnl/algo/vnl_matrix_inverse.h"
#include "vnl/vnl_transpose.h"
#include "vnl/vnl_matrix.h"
namespace otb
{
namespace Functor
{
/** \class TimeSeriesLeastSquareFittingFunctor
* \brief Implements a least squares fitting of a time profile
*
* This functor
* implements a least squares fitting of a time profile. The time
* series as well as the date series are supposed to accept the []
* syntax to get the values and the Size() method to get their
* length. The fitting is performed using a time function which is
* just a weighted sum of basis functions : \f$ f(t) = c_0 *
* \phi_0(t) + c_1 * \phi_1(t) + c_1 * \phi_1(t) + \cdots \f$ Using
* a matrix notation, this can be written as follows :
* \f$ F = \Phi c \f$
* \f$ F = ( f(t_i) ) \f$
* \f$ c = ( c_j ) \f$
* \f$ \Phi_{ij} = ( \phi_j(t_i) ) \f$
* If information about the error
* of each measure is available, it can be provided using a weight
* series \f$\sigma_i\f$ (the higher the error, the higher the \f$ \sigma_i \f$. The problem then becomes:
* \f$ b = A c \f$
* \f$ b = (\frac{ f(t_i) }{\sigma_i}) \f$
* \f$ A_{ij} = \frac{\Phi_{ij}}{\sigma_i} \f$
*
*
*
* \ingroup OTBTimeSeries
*/
template <class TSeriesType, class TTimeFunction, class TDateType = TSeriesType, class TWeightType = TSeriesType>
class TimeSeriesLeastSquareFittingFunctor
{
public:
typedef typename TTimeFunction::CoefficientsType CoefficientsType;
/// Constructor
TimeSeriesLeastSquareFittingFunctor()
{
for(unsigned int i=0; i<m_WeightSeries.Size(); ++i)
m_WeightSeries[i] = 1.0;
}
/// Destructor
virtual ~TimeSeriesLeastSquareFittingFunctor() {}
inline TSeriesType operator ()(const TSeriesType& series)
{
TTimeFunction estFunction = this->EstimateTimeFunction(series);
TSeriesType outSeries;
for(unsigned int i = 0; i < m_DoySeries.Size(); ++i)
outSeries[i] = estFunction.GetValue( m_DoySeries[i] );
return outSeries;
}
inline void SetDates(const TDateType& doy)
{
for(unsigned int i = 0; i < doy.Size(); ++i)
m_DoySeries[i] = doy[i];
}
inline void SetWeights(const TWeightType& weights)
{
for(unsigned int i=0; i<weights.Size(); ++i)
m_WeightSeries[i] = weights[i];
}
inline CoefficientsType GetCoefficients(const TSeriesType& series) const
{
return (this->EstimateTimeFunction(series)).GetCoefficients();
}
inline TTimeFunction EstimateTimeFunction(const TSeriesType& series) const
{
TTimeFunction estFunction;
unsigned int nbDates = m_DoySeries.Size();
unsigned int nbCoefs = estFunction.GetCoefficients().Size();
// b = A * c
vnl_matrix<double> A(nbDates, nbCoefs);
vnl_matrix<double> b(nbDates, 1);
// fill the matrices
typename TTimeFunction::CoefficientsType tmpCoefs;
for(unsigned int j = 0; j < nbCoefs; ++j)
tmpCoefs[j] = 0.0;
for(unsigned int i = 0; i < nbDates; ++i)
{
b.put(i, 0, series[i] / m_WeightSeries[i]);
for(unsigned int j = 0; j < nbCoefs; ++j)
{
tmpCoefs[j] = 1.0;
estFunction.SetCoefficients(tmpCoefs);
A.put(i, j, estFunction.GetValue(m_DoySeries[i]) / m_WeightSeries[i]);
tmpCoefs[j] = 0.0;
}
}
// solve the problem c = (At * A)^-1*At*b
vnl_matrix<double> atainv = vnl_matrix_inverse<double>(vnl_transpose(A) * A);
vnl_matrix<double> atainvat = atainv * vnl_transpose(A);
vnl_matrix<double> c = atainvat * b;
for(unsigned int j = 0; j < nbCoefs; ++j)
tmpCoefs[j] = c.get(j, 0);
estFunction.SetCoefficients(tmpCoefs);
return estFunction;
}
private:
///
TDateType m_DoySeries;
TWeightType m_WeightSeries;
};
}
} //namespace otb
#endif
|