void fillInConstraintMatrix( )
	{
		int N = listHInv.size();
		A.reshape(N*numberOfConstraints(),6);

		if( A.numRows < notZeros.size ) {
			throw new IllegalArgumentException("More unknowns than equations");
		}

		double b = knownAspectRatio ? aspectRatio*aspectRatio : 1;
		for (int i = 0,idx=0; i < N; i++) {
			Homography2D_F64 H = listHInv.get(i);

			if( zeroSkew ) { // w12
				A.data[idx++] = H.a11*H.a12;
				A.data[idx++] = H.a12*H.a21 + H.a11*H.a22;
				A.data[idx++] = H.a12*H.a31 + H.a11*H.a32;
				A.data[idx++] = H.a21*H.a22;
				A.data[idx++] = H.a22*H.a31 + H.a21*H.a32;
				A.data[idx++] = H.a31*H.a32;
			}
			if( principlePointOrigin ) { // w13 = 0, w23 = 0
				A.data[idx++] = H.a11*H.a13;
				A.data[idx++] = H.a13*H.a21 + H.a11*H.a23;
				A.data[idx++] = H.a13*H.a31 + H.a11*H.a33;
				A.data[idx++] = H.a21*H.a23;
				A.data[idx++] = H.a23*H.a31 + H.a21*H.a33;
				A.data[idx++] = H.a31*H.a33;

				A.data[idx++] = H.a12*H.a13;
				A.data[idx++] = H.a13*H.a22 + H.a12*H.a23;
				A.data[idx++] = H.a13*H.a32 + H.a12*H.a33;
				A.data[idx++] = H.a22*H.a23;
				A.data[idx++] = H.a23*H.a32 + H.a22*H.a33;
				A.data[idx++] = H.a32*H.a33;
			}

			if( knownAspectRatio ) { // w11 - b*w22 = 0
				A.data[idx++] = H.a11*H.a11 - H.a12*H.a12*b ;
				A.data[idx++] = 2*H.a11*H.a21 - 2*H.a12*H.a22*b;
				A.data[idx++] = 2*H.a11*H.a31 - 2*H.a12*H.a32*b ;
				A.data[idx++] = H.a21*H.a21 - H.a22*H.a22*b;
				A.data[idx++] = 2*H.a21*H.a31 - 2*H.a22*H.a32*b;
				A.data[idx++] = H.a31*H.a31 - H.a32*H.a32*b;
			}
		}
	}
void fillInConstraintMatrix ( ) { int N = listHInv . size ( ) ; A . reshape ( N * numberOfConstraints ( ) , 6 ) ; if ( A . numRows < notZeros . size ) { throw new IllegalArgumentException ( "More unknowns than equations" ) ; } double b = knownAspectRatio ? aspectRatio * aspectRatio : 1 ; for ( int i = 0 , idx = 0 ; i < N ; i ++ ) { Homography2D_F64 H = listHInv . get ( i ) ; if ( zeroSkew ) { A . data [ idx ++ ] = H . a11 * H . a12 ; A . data [ idx ++ ] = H . a12 * H . a21 + H . a11 * H . a22 ; A . data [ idx ++ ] = H . a12 * H . a31 + H . a11 * H . a32 ; A . data [ idx ++ ] = H . a21 * H . a22 ; A . data [ idx ++ ] = H . a22 * H . a31 + H . a21 * H . a32 ; A . data [ idx ++ ] = H . a31 * H . a32 ; } if ( principlePointOrigin ) { A . data [ idx ++ ] = H . a11 * H . a13 ; A . data [ idx ++ ] = H . a13 * H . a21 + H . a11 * H . a23 ; A . data [ idx ++ ] = H . a13 * H . a31 + H . a11 * H . a33 ; A . data [ idx ++ ] = H . a21 * H . a23 ; A . data [ idx ++ ] = H . a23 * H . a31 + H . a21 * H . a33 ; A . data [ idx ++ ] = H . a31 * H . a33 ; A . data [ idx ++ ] = H . a12 * H . a13 ; A . data [ idx ++ ] = H . a13 * H . a22 + H . a12 * H . a23 ; A . data [ idx ++ ] = H . a13 * H . a32 + H . a12 * H . a33 ; A . data [ idx ++ ] = H . a22 * H . a23 ; A . data [ idx ++ ] = H . a23 * H . a32 + H . a22 * H . a33 ; A . data [ idx ++ ] = H . a32 * H . a33 ; } if ( knownAspectRatio ) { A . data [ idx ++ ] = H . a11 * H . a11 - H . a12 * H . a12 * b ; A . data [ idx ++ ] = 2 * H . a11 * H . a21 - 2 * H . a12 * H . a22 * b ; A . data [ idx ++ ] = 2 * H . a11 * H . a31 - 2 * H . a12 * H . a32 * b ; A . data [ idx ++ ] = H . a21 * H . a21 - H . a22 * H . a22 * b ; A . data [ idx ++ ] = 2 * H . a21 * H . a31 - 2 * H . a22 * H . a32 * b ; A . data [ idx ++ ] = H . a31 * H . a31 - H . a32 * H . a32 * b ; } } }
Fill in A using all the constraints

Sage Math used for derivation
<pre>
H = matrix(SR, 3, 3, var('H11,H12,H13,H21,H22,H23,H31,H32,H33'))
var('W11,W12,W13,W22,W23,W33')
W = matrix(SR, 3, 3, [W11,W12,W13,W12,W22,W23,W13,W23,W33])
a=H.transpose()*W*H
a[0,1].expand()
</pre>