Vtyjkyuk

I have two point sets of about 800 3D-points that are matched. That means, I know the corresponding points in the two point clouds.
Now I want to minimize the distance between the corresponding points using a projective transformation because I know that the transformation can't be completly described by an affine transformation. So after using an affine transformation I still have an error that is a little too high.

Now my question is how I can compute a viable projective transformation that maps the points of one set to the points of the other set? I know how to calculate least square estimates of a rigid transformation and of an affine transformation but I can't find any sources how to handle two large point sets and compute something like a least squares solution of a projective transform. Is that possible?

At the moment my idea is the following:
1. Compute an affine transformation, using the least squares solution of:

$
beginpmatrix
x' \
y' \
z' \
1
endpmatrix
=
beginpmatrix
a_1 & a_2 & a_3 & a_4 \
b_1 & b_2 & b_3 & b_4 \
c_1 & c_2 & c_3 & c_4 \
0 & 0 & 0 & 1
endpmatrix
cdot
beginpmatrix
x \
y \
z \
1
endpmatrix
$

This matrix is used as an start vector for a minimization procedure using levenberg-marquardt on the following matrix:

$
T =
beginpmatrix
a_1 & a_2 & a_3 & a_4 \
b_1 & b_2 & b_3 & b_4 \
c_1 & c_2 & c_3 & c_4 \
d_1 & d_2 & d_3 & 1
endpmatrix
$

The idea is to minimize the cost function of the reprojection error:

$
sum | Tcdot x_i - x_i' |
$

The problem in my eyes is that I don't use a parametrization of the matrix that preserves non singularity in the optimization process which could be a problem because projective transformation matrices have to be non singular and in addition I thought that I don't need a non-linear optimization routine (because the projective transformation is still linear)!
Do you think this approach is viable or is there any better and easier way to determine a projective transformation between two 3D-point sets?

Actually I know that I only have to solve the following equation system:
$
beginpmatrix
wx' \
wy' \
wz' \
w
endpmatrix
=
beginpmatrix
a_1 & a_2 & a_3 & a_4 \
b_1 & b_2 & b_3 & b_4 \
c_1 & c_2 & c_3 & c_4 \
d_1 & d_2 & d_3 & 1
endpmatrix
cdot
beginpmatrix
x \
y \
z \
1
endpmatrix
$

So the difference to the affine transformations are the homogeneous coordinates. But how is this possible in a sense of regarding 800 points in a least squares way?

Best regards!