[Insight-users] linear solver for affine parameters

[Zoltan Seress]

Hi all,

The question is at the end of this letter, I would like to sketch the
problem before.

Well known, that the affine transformation that moves x voxel to y position
in 3D can be represented as a 4x4 homogenous matrix.

      ( a11 a12 a13 a14
A =   a21 a22 a23 a24
        a31 a32 a33 a34
         0     0    0     1  )

I set a region R in the image, the area I want to register. This R has N
voxel. A voxel can be represented by its coordinates, therefore I have N
3-elements vectors.
u1 = (x1 y1 z1 1), ..., un = (xn, yn, zn 1)

The same is true about the pairs of the voxels.
v1 = (x1' y1' z1' 1), ..., vn = (xn', yn', zn' 1)

I know all the u_i and v_i coordinates ( i=1,..,N ), what I need is the
transformation between them.

This can be write in matrix form for every voxel: A*u_i=v_i. If I modify
this formula that the 12 independent entries of the affine matrix in
homogeneous coordinates be the unknown parameters. For this I split up every
equation to 3 parts (because of the 3D), then I get (if Bx=y):

a11*x_i + a12*y_i + a13*z_i + a14 = x_i'
a21*x_i + a22*y_i + a23*z_i + a24 = y_i'
a31*x_i + a32*y_i + a33*z_i + a34 = z_i'

This is for one pair of voxels and if I want to write all equations in one
formula I should do that like this:

B = ( x1   y1   z1   1   0     0     0    0    0    0     0    0      x =
(a11       y = ( x1'
        0     0     0    0   x1   y1   z1   1    0    0     0    0
  a12               y1'
        0     0     0    0    0     0    0    0   x1   y1   z1   1
  a13               z1'

             a14
                                      ........
             a21               ...

             a22
        xn   yn   zn   1    0    0    0    0    0     0     0    0
  a23                xn'
         0     0    0    0   xn  yn   zn   1    0     0     0    0
  a24                yn'
         0     0    0    0    0    0    0    0   xn    yn   zn   1  )
   a31                zn' )

             a32

             a33

             a34 )

Here B is a 3N*12 matrix, x is a 12*1 vector and y is a 3N*1 vector. It is
obvious that B is almost never symmetric and it contain 33% non-zero
elements. The least trimmed squares method was proposed to solve this kind
of linear systems, but I couldn't find any c++ implementation (and I don't
intend to implement one), so would you suggest some reliable solvers for
this over-constrained system? Additional info that matrix B and vector y are
weighted with scalars come from former image analysis, but there is a chance
that some of the ( x_i, y_i ) voxel pairings are not the best (problem of
outliers), so for example the traditional least squares surely is not robust
enough.

Sorry for the too much lines. Would you be so kind to give me some idea?

Regards,
-- 
Zoli
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20100810/268f158a/attachment.htm>