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

FIX: correct OLS regression formula (switch X=Rainfall and Y=NDVI)

parent 343dcabc
......@@ -420,7 +420,7 @@ public:
// pvalue
OutputPixelValueType pvalue = 0.0;
if (vcl_abs(correlation) < 1)
if (vcl_abs(correlation) < 1 && nf > 2) // Degree of freedom > 0
{
// degrees of freedom
OutputPixelValueType df = nf - 2;
......@@ -586,7 +586,10 @@ public:
typedef typename TNDVIPixel::ValueType NDVIPixelValueType;
typedef typename TRainfallPixel::ValueType RainfallPixelValueType;
RainfallEstimatedNDVIResiduesFunctor(){}
RainfallEstimatedNDVIResiduesFunctor(){
m_NDVINoDataValue = 0;
m_RainfallNoDataValue = 0;
}
~RainfallEstimatedNDVIResiduesFunctor(){}
// Purposely not implemented
......@@ -609,6 +612,29 @@ public:
{
unsigned int n = inputNDVIPixel.Size();
// // TODO: delete this ////
// std::cout << "Operator ()" << std::endl;
// std::cout << "Input pixels:" << std::endl;
// std::cout << "ndvi=[";
// for (unsigned int i = 0 ; i < n ; i++)
// {
// if (inputNDVIPixel[i] != m_NDVINoDataValue && inputRainfallPixel[i] != m_RainfallNoDataValue)
// {
// std::cout << inputNDVIPixel[i] << ", ";
// }
// }
// std::cout << "];" << std::endl;
// std::cout << "rain=[";
// for (unsigned int i = 0 ; i < n ; i++)
// {
// if (inputNDVIPixel[i] != m_NDVINoDataValue && inputRainfallPixel[i] != m_RainfallNoDataValue)
// {
// std::cout << inputRainfallPixel[i] << ", ";
// }
// }
// std::cout << "];" << std::endl;
// /////////////////////////
// Prepare output pixel
TNDVIPixel residues;
residues.SetSize(n);
......@@ -616,7 +642,7 @@ public:
// 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++)
......@@ -627,7 +653,7 @@ public:
NDVIPixelValueType rainVal = static_cast<NDVIPixelValueType>(inputRainfallPixel[i]);
ndviSum += ndviVal;
rainSum += rainVal;
ndviSqSum += ndviVal*ndviVal;
rainSqSum += rainVal*rainVal;
coSum += ndviVal*rainVal;
count++;
}
......@@ -637,18 +663,18 @@ public:
{
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;
// OLS Y = slope * X + zeroy
NDVIPixelValueType slope = cov / var_ndvi;
NDVIPixelValueType zeroy = nfInv *( ndviSum - slope * rainSum );
NDVIPixelValueType slope = cov / var_rain;
NDVIPixelValueType zeroy = nfInv * ( slope * rainSum - ndviSum );
// 2. Compute residues
for (unsigned int i = 0 ; i < n ; i++)
......@@ -669,6 +695,17 @@ public:
residues.Fill(m_NDVINoDataValue);
}
// TODO: delete this
// std::cout << "res=[";
// for (unsigned int i = 0 ; i < n ; i++)
// {
// if (inputNDVIPixel[i] != m_NDVINoDataValue && inputRainfallPixel[i] != m_RainfallNoDataValue)
// {
// std::cout << residues[i] << ", ";
// }
// }
// std::cout << "];" << std::endl;
return residues;
}
......
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