<br clear="all">Hi all,<div><br></div><div>The question is at the end of this letter, I would like to sketch the problem before.</div><div><br></div><div><div>Well known, that the affine transformation that moves x voxel to y position in 3D can be represented as a 4x4 homogenous matrix.</div>
<div><br></div><div> ( a11 a12 a13 a14</div><div>A = a21 a22 a23 a24</div><div> a31 a32 a33 a34</div><div> 0 0 0 1 )</div></div><div><br></div><div><div>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.</div>
<div>u1 = (x1 y1 z1 1), ..., un = (xn, yn, zn 1)</div><div><br></div><div>The same is true about the pairs of the voxels.</div><div>v1 = (x1' y1' z1' 1), ..., vn = (xn', yn', zn' 1)</div></div><div>
<br></div><div>I know all the u_i and v_i coordinates ( i=1,..,N ), what I need is the transformation between them.</div><div><br></div><div>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):</div>
<div><br></div><div>a11*x_i + a12*y_i + a13*z_i + a14 = x_i'</div><div>a21*x_i + a22*y_i + a23*z_i + a24 = y_i'</div><div>a31*x_i + a32*y_i + a33*z_i + a34 = z_i'</div><div><br></div><div>This is for one pair of voxels and if I want to write all equations in one formula I should do that like this:</div>
<div><br></div><div><div>B = ( x1 y1 z1 1 0 0 0 0 0 0 0 0 x = (a11 y = ( x1'</div><div> 0 0 0 0 x1 y1 z1 1 0 0 0 0 a12 y1'</div>
<div> 0 0 0 0 0 0 0 0 x1 y1 z1 1 a13 z1'</div><div> a14</div><div>
........ a21 ...</div><div> a22</div>
<div> xn yn zn 1 0 0 0 0 0 0 0 0 a23 xn' </div><div> 0 0 0 0 xn yn zn 1 0 0 0 0 a24 yn'</div>
<div> 0 0 0 0 0 0 0 0 xn yn zn 1 ) a31 zn' ) </div></div><div> a32</div>
<div> a33</div><div> a34 )</div><div><br></div>
<div>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.</div>
<div><br></div><div>Sorry for the too much lines. Would you be so kind to give me some idea?</div><div><br></div><div>Regards,<br>-- <br>Zoli<br>
</div>