Hi, I'm researching camera calibration and uncertainty propagation. I'm trying to understand how calibrateCameraExtended estimates stdDeviationsIntrinsics and stdDeviationsExtrinsics. The docs say very little.
So I go straight to the source code. Check this lines of code where I understand the calculation is made.
It starts with sigma2 (the "deviation of the noise") calculated as norm(allErrors, NORM_L2SQR) / (total - nparams_nz);
which is just the formula for the unbiased estimator of the variance. Ok so far.
And then it calculates each s-element of the vector of standard deviations stdDevsM by
stdDevsM.at<double>(s) = std::sqrt(JtJinv.at<double>(j,j) * sigma2);
Where JtJinv is the pseudo-inverse of the jacobian calculated a few lines above from _JtJ which in turn comes from the
LM solver invoked in previous lines.
First question: what exactly is _JtJ? I assume it must be the 1xN Jacobian of the projection error with respect to the parameters (there are N parameters). I've tried to trace the calculation of all the way to its origin, I got
this far
, but I'm not sure.
Second question: The moore-penrose of a 1xN matrix is a Nx1 matrix. So calling JtJinv.at<double>(j,j) with two indices j confuses me.
Third question: I couldn't make sense of the formula itself and I think It's wrong. The code has the comment
//see any papers about variance of the least squares estimator for
//detailed description of the variance estimation methods
But my understanding from simple uncertainty propagation for the case of uncorrelated parameters is that the Jacobian vector _JtJ, the variance of the projection error sigma2 and the vector of parameters standard deviations stdDevsM should follow (in matlab-like pseudocode):
sigma2 = dot(_JtJ.^2, stdDevsM.^2)
Which is ill-conditioned, the simplest solution would be
stdDevsM = sqrt(sigma2) .*_JtJ ./ norm(_JtJ.^2)
Is my reasoning correct? and, where does the calculation implemented in OpenCV come from?