[Insight-users] Similarity3DTransform : Jacobian computation

[Luis Ibanez]

Hi Heike,


            Thanks for pointing this out.


In my previous email I put derivatives out of place in the
description of operations.

I derived prematurely the second part of that expression.


Let me describe the operations here again,
in a cleaner way:


If we take the [0][1] term of the Rotation matrix:

          f(x,y,z) =  2xy - 2wz

where    w = sqrt( 1 - x.x - y.y - z.z )

then

   f(x,y,z) =  2xy - 2 sqrt(1-x.x-y.y-z.z) z

   f(x,y,z) =  2xy - 2 (1-x.x-y.y-z.z)^(1/2) z

   f(x,y,z) =  2xy - 2z (1-x.x-y.y-z.z)^(1/2)


since

  d( G(p)^(1/2) )/dp = (1/2)( (G(p)^(-1/2)).d(G(p))/dp

then

df(x,y,z)/dx = 2y - 2z . (1/2).(1-x.x-y.y-z.z)^(-1/2).( -2x )


since w = sqrt(1-x.x-y.y-z.z)= (1-x.x-y.y-z.z)^(1/2)

then  (1-x.x-y.y-z.z)^(-1/2) = 1/w


therefore

   df(x,y,z)/dx = 2y - 2z . (1/2).(1/w).( -2x )

   df(x,y,z)/dx = 2y - 2z . (-2x)/(2w)

   df(x,y,z)/dx = 2y - 2z . (-x/w)

   df(x,y,z)/dx = 2y + 2 zx / w

   df(x,y,z)/dx = 2yw/w + 2zx/w

   df(x,y,z)/dx = ( 2yw + 2zx ) / w

   df(x,y,z)/dx = 2 ( yw + zx ) / w


Which matches the first term of line  174 in
Insight/Code/Common/itkVersorRigid3DTransform.txx

this->m_Jacobian[0][0] = 2.0 * ( (vyw+vxz)*py + (vzw-vxy)*pz) 
          / vw;

That is, the multiplier of "py" is

       2 ( vyw + vxz ) / vw

where

vyw = vy.vw
vxz = vx.vz

and {vx,vy,vz,vw} are the Versor components expressed
as a Quaternion.


------


BTW, In my previous email, I also missed the sign in the step
from:

 >>                      =  2xy - 2z (1/2) ( -2 x ) / w
 >>

to:

 >>
 >>                         2y - 2xz / w


It should have been (as shown in this new email):

                        =   2y + 2xz / w

                        =  ( 2yw + 2xz ) / w



    My apologies for the mishap,




       Regards,


          Luis


-----------------
Heike Moll wrote:
> Hi Luis,
> 
> thanks for your quick respond.
> 
> Please can you discribe the following step of your calculation in detail.
> 
>   2xy - 2 sqrt(1-x.x-y.y-z.z) z =  2xy - 2z (1/2) ( -2 x ) / w
> 
> Thx Heike
> 
> ----- Ursprüngliche Nachricht -----
> Von: Luis Ibanez <luis.ibanez at kitware.com>
> Datum: Mittwoch, Juni 4, 2008 10:58 pm
> Betreff: Re: [Insight-users] Similarity3DTransform
> An: Heike Moll <Heike.Moll at rwth-aachen.de>
> Cc: insight users <insight-users at itk.org>
> 
>>Hi Heike,
>>
>>
>>1) We will discourage you from including the coordinates of the center
>>   in the list of transform parameters. We have found in the past 
>>that    this make the optimization unstable...   :-/
>>
>>
>>2) About the Jacobian computation:
>>   The rotation matrix is the one in page 61 of
>>   http://www.itk.org/CourseWare/Training/QuaternionsI.pdf
>>
>>
>>   The important feature to keep in mind when deriving with respect
>>   to Versor components, as opposed to deriving with respect to
>>   Quaternion components, is that in the Versor, the "w" component
>>   is not independent.  That is the four values {x,y,z,w} are
>>   actually {x,y,z,sqrt(1-xx-yy-zz)} (where xx = x*x).
>>
>>
>>   Therefore, the derivative with respect to "x" of a term such as
>>
>>                         2xy - 2wz
>>
>>   (the x,y term of the Rotation matrix).
>>
>>   that in the case of a Quaternion will yield:
>>
>>                         2y
>>
>>   will instead, in the case of a Versor requires to express W
>>   in terms of x,y,z as w = sqrt( 1 - x.x - y.y - z.z )
>>
>>                         2xy - 2wz
>>                      =  2xy - 2 sqrt(1-x.x-y.y-z.z) z
>>                      =  2xy - 2z (1/2) ( -2 x ) / w
>>
>>     deriving this with respect to x will yield:
>>
>>                         2y - 2xz / w
>>
>>      which is the same as
>>
>>                        ( 2wy - 2xz ) / w
>>
>>      as you will find in the first part of line 214 in
>>
>>
>>                itkVersorRigid3DTransform.txx
>>
>>
>>
>>
>>   Please let us know if you find any flaw in the math,
>>
>>
>>      Thanks
>>
>>
>>
>>         Luis
>>
>>
>>
>>-------------------
>>Heike Moll wrote:
>>
>>>Hi,
>>>
>>>I am trying to implement a class named 
>>
>>CenteredSimilarity3DTransform. 
>>
>>>This class is derived from Similarity3DTransform and is suposed 
>>
>>to optimize the center of transformation (like 
>>CenteredSimilarity2DTransform).> My problem is the method 
>>GetJacobian( const InputPointType & p )!!
>>
>>>I have to rewrite the Jacobian matrix for 3D, so I tried to 
>>
>>understand the calculation in the VersorRigid3DTransform class.
>>
>>>You used the equation P`=R*P+T.
>>>
>>>But if I recalculate the jacobian matrix (whose elements are the 
>>
>>partial derivatives of the output point with respect to the 
>>transform parameters) with above equation and the (quaternion) 
>>rotation matrix R:
>>
>>>R=1-2(y²+z²), 2xy-2wz, 2xz+2wy;
>>>     2xy+2wz, 1-2(x²+z²), 2yz-2wx;
>>>     2xz-2wy, 2yz+2wx, 1-2(x²+y²),
>>>
>>>I get another result as the one used in VersorRigid3DTransform. 
>>
>>The versor differs in all positions accept in the diagonal. 
>>
>>>Why? Did I use the wrong rotation matrix R? Or have I a logical 
>>
>>mistake in my thoughts?
>>
>>>Thx Heike
>>>_______________________________________________
>>>Insight-users mailing list
>>>Insight-users at itk.org
>>>http://www.itk.org/mailman/listinfo/insight-users
>>>
>>
>