[ITK-users] How to implement VersorRigid3DTransform ?

[suguru KONDO]

Hello everyone.

I have a couple of questions about VersorRigid3DTransform.
I will sequentially consider registration process using VersorRigid3DTransform. 
If, my concept is wrong, please point it.

1.Transform moving image to fixed image

I think the formulation of VersorRigid3DTransform is below
and then all input points in moving image are transformed by this formulation.

[x']   [q1^2 - q2^2 - q3^2 + q4^2, 2(q1q2 - q3q4), 2(q1q3 + q2q4)]  [x - xc]   [tx]
[y'] = [2(q1q2 + q3q4), -q1^2 + q2^2 - q3^2 + q4^2, 2(q2q3 - q1q4)] [y - yc] + [ty]
[z']   [2(q1q3 - q2q4), 2(q2q3 + q1q4), -q1^2 - q2^2 + q3^2 + q4^2] [z - zc]   [tz]

P' = R * (P - C) + T

Versor : Q = q1 * i + q2 * j + q3 * k + q4 (q1^2 + q2^2 + q3^2 + q4^2 = 1)
Input point : P = (x, y, z)
Output point : P' = (x', y', z')
Center of rotation : C = (xc, yc, zc)
Translation vector : T = (tx, ty, tz) 
Rotation matrix : R 
Optimization parameter : {q1, q2, q3, tx, ty, tz}

2.Compute metric (cost function) and its derivative
For simplicity, I assumed that interpolation completed and metric is mean squares metric.  

In this case, I think metric is formulated by this formula.

MS(f(X), Tm(X)) = 1/N * Σ(i=1...N)(f(X)_i  - Tm(X)_i )^2

All points in fixed image region : X
Metric : MS(f(X), Tm(X))
Fixed image : f(X)
Transformed moving image : Tm(X)
the i-th pixel of Image f(X) : f(X)_i 
the i-th pixel of Image Tm(X) : Tm(X)_i
the number of pixels considered : N = the number of pixels in fixed image (when calculation region is GetLargestPossibleRegion())

I think Tm(X)_i is represented by using optimization parameters {q1, q2, q3, tx, ty, tz}.
So, I think its jacobian matrix can be represented like this.

    [∂Tm(X)_1/∂q1, ∂Tm(X)_1/∂q2, ∂Tm(X)_1/∂q3, ∂Tm(X)_1/∂tx, ∂Tm(X)_1/∂ty, ∂Tm(X)_1/∂tz]
    [∂Tm(X)_2/∂q1, ∂Tm(X)_2/∂q2, ∂Tm(X)_2/∂q3, ∂Tm(X)_2/∂tx, ∂Tm(X)_2/∂ty, ∂Tm(X)_2/∂tz]
    [        .               .               .               .               .               .     ]
 J =[        .               .               .               .               .               .     ]
    [        .               .               .               .               .               .     ]
    [∂Tm(X)_N/∂q1, ∂Tm(X)_N/∂q2, ∂Tm(X)_N/∂q3, ∂Tm(X)_N/∂tx, ∂Tm(X)_N/∂ty, ∂Tm(X)_N/∂tz]

Here, I have a couple of questions.

q1, q2, q3 are versor components in versor space but I deal them same as tx, ty, tz.
Is my operation correct?   

Then, if this jacobian matrix is correct, we need to compute metric derivative by using this
because we need to decide search direction in optimization parameter space. 
But, I don't know how to compute metric derivative using jacobian matrix.
So, is there any publication outlining the exact calculations leading 
to the metric derivative implemented in the VersorRigid3DTransform?

I hope that somebody finds the time to answer.

Thanks

Suguru Kondo