Possible error in least-squares plane fit coefficients in PlaneFit.cc
While reviewing PlaneFit.cc, I noticed what appears to be two issues in the least-squares implementation used to fit:
$$
z = ax + by + c
$$
The code computes the centered covariance terms:
xx = Σ(dx²)
xy = Σ(dx dy)
xz = Σ(dx dz)
yy = Σ(dy²)
yz = Σ(dy dz)
which leads to the normal equations:
$$
\begin{aligned}
xx,a + xy,b &= xz \\
xy,a + yy,b &= yz
\end{aligned}
$$
The inverse of the covariance matrix is:
$$
\frac{1}{xx \cdot yy - xy^2}
\begin{bmatrix}
yy & -xy \\
-xy & xx
\end{bmatrix}
$$
Multiplying by the right-hand side gives:
$$
a = \frac{yy \cdot xz - xy \cdot yz}{xx \cdot yy - xy^2}
$$
$$
b = \frac{xx \cdot yz - xy \cdot xz}{xx \cdot yy - xy^2}
$$
However, the current implementation is:
_aa = (xx * xz - xy * yz) / (xx * yy - xy * xy);
_bb = (xx * yz - xy * xz) / (xx * yy - xy * xy);
The _bb calculation matches the derived expression, but _aa appears to use xx instead of yy in the numerator:
_aa = (yy * xz - xy * yz) / (xx * yy - xy * xy);
Additionally, the intercept is currently set as:
whereas the least-squares solution for the full plane fit should be:
$$
c = \bar{z} - a\bar{x} - b\bar{y}
$$
or:
_cc = meanz - _aa * meanx - _bb * meany;
From looking at ECCO, it appears that only the fitted slopes are used, so the intercept issue likely has little practical effect there. The slope issue may also be partially masked when fitting symmetric neighborhoods, but it could affect results for asymmetric point distributions (for example, when valid points are missing from part of the fitting window). I've noticed this affects precipitation edges where non-valid echo arises. I'm by no means a math expert so I'd appreciate some double-checks. For most uses within ECCO, this does not make a meaningful difference.
Could someone confirm whether the current implementation is intentional, or whether _aa and _cc should be updated to match the standard least-squares solution? I included an image of a Pythonic port (which has the updated _aa) that I've worked on and you can see where texture slightly differs. Again, it's almost entirely edge effects and it only impacts texture by ~10% (not very meaningful, but can be interesting).
Possible error in least-squares plane fit coefficients in
PlaneFit.ccWhile reviewing
PlaneFit.cc, I noticed what appears to be two issues in the least-squares implementation used to fit:The code computes the centered covariance terms:
which leads to the normal equations:
The inverse of the covariance matrix is:
Multiplying by the right-hand side gives:
However, the current implementation is:
The
_bbcalculation matches the derived expression, but_aaappears to usexxinstead ofyyin the numerator:Additionally, the intercept is currently set as:
whereas the least-squares solution for the full plane fit should be:
or:
From looking at ECCO, it appears that only the fitted slopes are used, so the intercept issue likely has little practical effect there. The slope issue may also be partially masked when fitting symmetric neighborhoods, but it could affect results for asymmetric point distributions (for example, when valid points are missing from part of the fitting window). I've noticed this affects precipitation edges where non-valid echo arises. I'm by no means a math expert so I'd appreciate some double-checks. For most uses within ECCO, this does not make a meaningful difference.
Could someone confirm whether the current implementation is intentional, or whether
_aaand_ccshould be updated to match the standard least-squares solution? I included an image of a Pythonic port (which has the updated_aa) that I've worked on and you can see where texture slightly differs. Again, it's almost entirely edge effects and it only impacts texture by ~10% (not very meaningful, but can be interesting).