Commit 8c69473e authored by remi cresson's avatar remi cresson
Browse files

ADD: new functors + 1 test app

parent 91b01c27
cmake_minimum_required (VERSION 2.8)
OTB_CREATE_APPLICATION(NAME NDVITimeSeriesIndices
SOURCES otbNDVITimeSeriesIndices.cxx
LINK_LIBRARIES OTBCommon)
#include "itkFixedArray.h"
#include "itkObjectFactory.h"
// Elevation handler
#include "otbWrapperElevationParametersHandler.h"
#include "otbWrapperApplicationFactory.h"
// Application engine
#include "otbStandardFilterWatcher.h"
#include "itkFixedArray.h"
// Functors
#include "otbNDVITimeSeriesFunctor.h"
namespace otb
{
namespace Wrapper
{
class NDVITimeSeriesIndices : public Application
{
public:
/** Standard class typedefs. */
typedef NDVITimeSeriesIndices Self;
typedef Application Superclass;
typedef itk::SmartPointer<Self> Pointer;
typedef itk::SmartPointer<const Self> ConstPointer;
/** Standard macro */
itkNewMacro(Self);
itkTypeMacro(NDVITimeSeriesIndices, Application);
/** functors */
typedef otb::functor::SlopeAndPValueFunctor<FloatVectorImageType::PixelType, FloatVectorImageType::PixelType> SAPVFunctorType;
private:
void testFunctor()
{
SAPVFunctorType f;
FloatVectorImageType::PixelType p,y;
unsigned int n = 10;
p.SetSize(n);
y.SetSize(n);
for (unsigned int i = 1 ; i <= n ; i++)
{
p[i] = i*2;
}
std::cout << f(p) << std::endl;
}
void DoInit()
{
testFunctor();
SetName("NDVITimeSeriesIndices");
SetDescription("Perform NDVITimeSeriesIndices");
// Documentation
SetDocName("NDVITimeSeriesIndices");
SetDocLongDescription("This application performs NDVITimeSeriesIndices of images");
SetDocLimitations("None");
SetDocAuthors("Remi Cresson");
SetDocSeeAlso(" ");
// Input image
AddParameter(ParameterType_InputImageList, "il", "Input images");
// Output image
AddParameter(ParameterType_OutputImage, "out", "Output image");
SetDefaultOutputPixelType("out", ImagePixelType_uint8);
AddRAMParameter();
}
void DoUpdateParameters()
{
// Nothing to do here : all parameters are independent
}
void DoExecute()
{
} // DOExecute()
void AfterExecuteAndWriteOutputs()
{
// Nothing to do
}
};
}
}
OTB_APPLICATION_EXPORT( otb::Wrapper::NDVITimeSeriesIndices )
/*=========================================================================
Copyright (c) IRSTEA. All rights reserved.
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 __otbNDVITimeSeriesFunctor_h
#define __otbNDVITimeSeriesFunctor_h
#include "vnl/vnl_vector.h"
#include "boost/date_time/gregorian/gregorian.hpp"
#include "boost/math/distributions/students_t.hpp"
#include <vector>
#include "assert.h"
namespace otb
{
namespace functor
{
class SingleDate
{
public:
SingleDate(int yyyy, int mm, int dd) {
year = yyyy;
day = dd;
month = mm;
// Compute julianday
boost::gregorian::date d(yyyy, mm, dd);
julianday = d.julian_day();
}
~SingleDate (){}
bool IsInRange(const int& startMonth,const int& startDay,const int& endMonth, const int& endDay) const
{
if (startMonth < month && month < endMonth)
{
return true;
}
else if (startMonth == month && startDay <= day)
{
return true;
}
else if (endMonth == month && day <= endDay)
{
return true;
}
else
{
return false;
}
}
int day;
int month;
int year;
int julianday;
};
/**
* \class TSReduceFunctorBase
* \brief Base class for reducing N-components pixel into M-component pixel (N=nb. of images and M=nb. of years)
*
* Dates must be provided in chronological order, else the "SetDates()" method will throw an exception.
*/
template< class TInputPixel, class TOutputPixel>
class TSReduceFunctorBase
{
public:
typedef std::vector<SingleDate> DatesType;
typedef typename TInputPixel::ValueType PixelValueType;
TSReduceFunctorBase() {}
virtual ~TSReduceFunctorBase() {}
bool operator!=( const TSReduceFunctorBase & ) const { return false; }
bool operator==( const TSReduceFunctorBase & other ) const { return !(*this != other); }
virtual inline TOutputPixel operator()( const TInputPixel & A ) const = 0;
inline unsigned int GetOutputSize(){ return nbOfYears; }
/* Set dates */
void SetDates(DatesType dates)
{
unsigned int n = dates.size();
// Check that dates are ordered
for (unsigned int i = 0 ; i < n - 1 ; i++)
{
assert(dates[i].julianday < dates[i+1].julianday);
}
// Assign dates
m_Dates = dates;
// Compute the number of years
nbOfYears = m_Dates[n-1].year - m_Dates[0].year + 1;
}
/* Get dates */
DatesType GetDates(){ return m_Dates; }
/* Set / Get input no-data value */
void SetInputNoDataValue(PixelValueType value) { m_InputNoDataValue = value; }
PixelValueType GetInputNoDataValue(){ return m_InputNoDataValue; }
private:
unsigned int nbOfYears;
DatesType m_Dates;
PixelValueType m_InputNoDataValue;
}; // TSReduceFunctorBase
/**
* \class TSReduceFromRangeFunctor
* \brief Class for reducing N-components pixel into M-component pixel (N=nb. of images and M=nb. of years)
* given a range (dd/mm--->dd'/mm')
*/
template< class TInputPixel, class TOutputPixel>
class TSReduceFromRangeFunctor : TSReduceFunctorBase<TInputPixel,TOutputPixel>
{
public:
typedef typename TSReduceFunctorBase<TInputPixel,TOutputPixel>::DatesType DatesType;
typedef typename TInputPixel::ValueType PixelValueType;
TSReduceFromRangeFunctor() {}
~TSReduceFromRangeFunctor() {}
bool operator!=( const TSReduceFromRangeFunctor & ) const {
return false;
}
bool operator==( const TSReduceFromRangeFunctor & other ) const {
return !(*this != other);
}
inline TOutputPixel operator ()(const TInputPixel& input) const
{
TOutputPixel output;
output.SetSize(this->GetOutputSize());
output.Fill(0);
DatesType m_Dates = this->GetDates();
int firstYear = m_Dates[0].year;
for (unsigned int i = 0; i < input.Size(); i++)
{
int index = m_Dates[i].year - firstYear;
int m = m_Dates[i].month;
PixelValueType inValue = input[i];
// Range check
if (inValue != this->GetInputNoDataValue())
{
if (m_Month1 <= m && m <= m_Month2)
{
int d = m_Dates[i].day;
if (m_Day1 <= d && d <= m_Day2)
{
output[index] += inValue;
}
}
}
}
return output;
}
void SetDay1(int day) { m_Day1 = day; }
void SetDay2(int day) { m_Day2 = day; }
void SetMonth1(int month) { m_Month1 = month; }
void SetMonth2(int month) { m_Month2 = month; }
private:
int m_Day1;
int m_Day2;
int m_Month1;
int m_Month2;
}; // TSReduceFromRangeFunctor
template< class TInputPixel, class TOutputPixel>
class SlopeAndPValueFunctor
{
public:
typedef typename TInputPixel::ValueType InputPixelValueType;
typedef typename TOutputPixel::ValueType OutputPixelValueType;
SlopeAndPValueFunctor(){}
~SlopeAndPValueFunctor(){}
bool operator!=( const SlopeAndPValueFunctor & ) const { return false; }
bool operator==( const SlopeAndPValueFunctor & other ) const { return !(*this != other); }
inline unsigned int GetOutputSize(){ return 2; }
inline TOutputPixel operator ()(const TInputPixel& input) const
{
TOutputPixel outpixel;
outpixel.SetSize(2);
outpixel.Fill(0);
// compute sum & sum of squares
unsigned int n = input.Size();
OutputPixelValueType sum = 0;
OutputPixelValueType isum = 0;
OutputPixelValueType sqSum = 0;
OutputPixelValueType isqSum = 0;
OutputPixelValueType coSum = 0;
unsigned int count = 0;
for (unsigned int i = 0 ; i < n ; i++)
{
if (input[i] != m_InputNoDataValue)
{
OutputPixelValueType val = static_cast<OutputPixelValueType>(input[i]);
OutputPixelValueType ival = static_cast<OutputPixelValueType>(i);
sum += val;
isum += ival;
sqSum += val*val;
isqSum += ival*ival;
coSum += val*ival;
count++;
}
}
if (count > 1)
{
OutputPixelValueType nf = static_cast<OutputPixelValueType>(count);
OutputPixelValueType nfInv = -1.0 / nf;
OutputPixelValueType sum2 = sum * sum;
OutputPixelValueType isum2 = isum * isum;
OutputPixelValueType sumisum = sum * isum;
// Covariance
OutputPixelValueType cov = coSum + nfInv * sumisum;
// Variance(t)
OutputPixelValueType var_t = isqSum + nfInv * isum2;
// Variance(x)outpixel[0]
OutputPixelValueType var_x = sqSum + nfInv * sum2;
// OLS
OutputPixelValueType slope = cov / var_t;
// Pearson linear correlation coefficient
OutputPixelValueType correlation = cov / vcl_sqrt(var_t * var_x);
// pvalue
OutputPixelValueType pvalue = 0.0;
if (correlation < 0.9999999)
{
// degrees of freedom
OutputPixelValueType df = nf - 2;
OutputPixelValueType t = correlation;
t /= sqrt((1.0 - correlation * correlation) / df);
// Compute student cumulative distribution function
boost::math::students_t_distribution<OutputPixelValueType> dist(df);
OutputPixelValueType cdf = boost::math::cdf(dist, vcl_abs(t));
// Compute p-value
pvalue = 2.0 * (1.0 - cdf);
}
// output pixel
outpixel[0] = slope;
outpixel[1] = pvalue;
}
return outpixel;
}
/* Set / Get input no-data value */
void SetInputNoDataValue(InputPixelValueType value) { m_InputNoDataValue = value; }
InputPixelValueType GetInputNoDataValue(){ return m_InputNoDataValue; }
InputPixelValueType m_InputNoDataValue;
}; // SlopeAndPValueFunctor
} // namespace functor
} // namespace otb
#endif
......@@ -5,8 +5,6 @@ otb_module(NDVITimeSeries
OTBITK
OTBCommon
OTBApplicationEngine
OTBTimeSeries
MultitemporalSmoothing
TEST_DEPENDS
OTBTestKernel
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment