Commit c42ac2ce authored by remi cresson's avatar remi cresson
Browse files

ADD: new functor PearsonCorrelationCoefficientFunctor

parent 1ce966a8
......@@ -840,6 +840,108 @@ protected:
}; // FactorsLabelingFunctor
/*
* \class PearsonCorrelationCoefficientFunctor
* \brief Computes the Pearson Correlation Coefficient
*
*/
template< class TNDVIPixel, class TRainfallPixel, class TCorrelationPixel>
class PearsonCorrelationCoefficientFunctor
{
public:
typedef typename TNDVIPixel::ValueType NDVIPixelValueType;
typedef typename TRainfallPixel::ValueType RainfallPixelValueType;
typedef typename TCorrelationPixel::ValueType CorrelationPixelValueType;
PearsonCorrelationCoefficientFunctor(){
m_CorrelationNoDataValue = itk::NumericTraits<CorrelationPixelValueType>::NonpositiveMin();
}
~PearsonCorrelationCoefficientFunctor(){}
// Purposely not implemented
bool operator!=( const PearsonCorrelationCoefficientFunctor & ) const { return false; }
bool operator==( const PearsonCorrelationCoefficientFunctor & other ) const { return !(*this != other); }
/* Set / Get NDVI input no-data value */
void SetNDVINoDataValue(NDVIPixelValueType value) { m_NDVINoDataValue = value; }
NDVIPixelValueType GetNDVINoDataValue() const { return m_NDVINoDataValue; }
/* Set / Get Rainfall input no-data value */
void SetRainfallNoDataValue(RainfallPixelValueType value) { m_RainfallNoDataValue = value; }
RainfallPixelValueType GetRainfallNoDataValue() const { return m_RainfallNoDataValue; }
/* Set / Get correlation input no-data value */
void SetCorrelationNoDataValue(CorrelationPixelValueType value) { m_CorrelationNoDataValue = value; }
CorrelationPixelValueType GetCorrelationNoDataValue() const { return m_CorrelationNoDataValue; }
/*
* Compute output pixel.
* inputNDVIPixel and inputRainfallPixel must have the same number of components.
*/
inline TCorrelationPixel operator ()(const TNDVIPixel& inputNDVIPixel, const TRainfallPixel& inputRainfallPixel) const
{
unsigned int n = inputNDVIPixel.Size();
// Prepare output pixel
TCorrelationPixel correlation;
correlation.SetSize(n);
// 1. Compute regression between Y=NDVI and X=Rainfall
NDVIPixelValueType ndviSum = 0;
NDVIPixelValueType rainSum = 0;
NDVIPixelValueType ndviSqSum = 0;
NDVIPixelValueType rainSqSum = 0;
NDVIPixelValueType coSum = 0;
unsigned int count = 0;
for (unsigned int i = 0 ; i < n ; i++)
{
if (inputNDVIPixel[i] != m_NDVINoDataValue && inputRainfallPixel[i] != m_RainfallNoDataValue)
{
NDVIPixelValueType ndviVal = inputNDVIPixel[i];
NDVIPixelValueType rainVal = static_cast<NDVIPixelValueType>(inputRainfallPixel[i]);
ndviSum += ndviVal;
rainSum += rainVal;
ndviSqSum += ndviVal*ndviVal;
rainSqSum += rainVal*rainVal;
coSum += ndviVal*rainVal;
count++;
}
}
if (count > 1)
{
NDVIPixelValueType nf = static_cast<NDVIPixelValueType>(count);
NDVIPixelValueType nfInv = -1.0 / nf;
NDVIPixelValueType ndviSum2 = ndviSum * ndviSum;
NDVIPixelValueType rainSum2 = rainSum * rainSum;
NDVIPixelValueType prodsums = ndviSum * rainSum;
// Covariance
NDVIPixelValueType cov = coSum + nfInv * prodsums;
// Variance(ndvi)
NDVIPixelValueType var_ndvi = ndviSqSum + nfInv * ndviSum2;
// Variance(rain)
NDVIPixelValueType var_rain = rainSqSum + nfInv * rainSum2;
// Pearson correlation coefficient
correlation[0] = cov / vcl_sqrt(var_ndvi * var_rain);
}
else
{
correlation[0] = m_CorrelationNoDataValue;
}
return correlation;
}
protected:
NDVIPixelValueType m_NDVINoDataValue;
RainfallPixelValueType m_RainfallNoDataValue;
CorrelationPixelValueType m_CorrelationNoDataValue;
}; // PearsonCorrelationCoefficientFunctor
} // namespace functor
} // namespace otb
......
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