jacobian in LM optimization in findhomography

edit

I'm trying to understand the code of findhomography function, while for the pose refine after ransac/prosac, ie. the createLMSolver(makePtr<HomographyRefineCallback>(src, dst), 10)->run(H8); line in findHomography() function, I'm confused about the computation of the jacobian (in the compute() function of HomographyRefineCallback when implementing the LM solver:

Jptr[0] = Mx*ww; Jptr[1] = My*ww; Jptr[2] = ww;
Jptr[3] = Jptr[4] = Jptr[5] = 0.;
Jptr[6] = -Mx*ww*xi; Jptr[7] = -My*ww*xi;
Jptr[8] = Jptr[9] = Jptr[10] = 0.;
Jptr[11] = Mx*ww; Jptr[12] = My*ww; Jptr[13] = ww;
Jptr[14] = -Mx*ww*yi; Jptr[15] = -My*ww*yi;

The source and dest points are 2d points, and the homography matrix to refine is a 3*3 matrix (opencv sets the last element to 1, so h has 8 elements here), so why the jacobian here has 16 elements? Could anyone explain it, or provide the original paper on this solution here?

Thanks!