Transform stack using eigenvectors

[Michael Doube]

Hi all

I have written a plugin that calculates the 3D moments of inertia of a
bone imaged in CT, using the Jama package and eigen decomposition.

Now that I have a set of eigenvalues and eigenvectors I want to align my
object such that its principal axes are parallel with the x,y,z axes of
a new stack.

Can anyone shed some light as to how to best achieve this?

Mike

[Thomas Boudier]

Hi Michael,

Maybe you can use the align stack plugin :
http://www.med.harvard.edu/JPNM/ij/plugins/AlignStacks.html

a bit complex to use but quite powerful !

Thomas

> Hi all
>
> I have written a plugin that calculates the 3D moments of inertia of a
> bone imaged in CT, using the Jama package and eigen decomposition.
>
> Now that I have a set of eigenvalues and eigenvectors I want to align
> my object such that its principal axes are parallel with the x,y,z
> axes of a new stack.
>
> Can anyone shed some light as to how to best achieve this?
>
> Mike
>
>

--
/*****************************************************/
 Thomas Boudier, MCU Université Pierre et Marie Curie
 UMR 7101 / IFR 83. Bat A 328, Campus Jussieu
 Tél : [phone]  Fax : [phone]
/*****************************************************/


 

[Michael Doube]

To answer my own question:

The aligned position is the 3 dot products of the original coordinate
(x,y,z) and the 3 eigenvectors (the principal axes).  Each eigenvector
is a normal to the plane that contains the other 2 eigenvectors because
they are orthogonal.  A plane is defined by the vector of its normal and
the position of the plane on that normal.  The dot product of a plane
and a point is the distance between plane and point, assuming that
(0,0,0) is a solution for the plane.  So if you subtract the centroid
from the point and work out the 3 dot products corresponding to the 3
eigenplanes, you have the distance along each of the eigenvectors which
defines the point in the principal axis coordinate frame, and hence
gives you 'aligned' coordinates.  Then the aligned coordinates can be
drawn in ordinary image coordinates and the object is aligned according
to its principal axes.

If that makes no sense or is blatantly wrong, please email me, otherwise
it's going in a plugin tonight...

Mike

> -----Original Message-----
> From: ImageJ Interest Group [mailto:[hidden email]] On Behalf Of
> Michael Doube
> Sent: Friday, August 15, 2008 5:57 PM
> To: [hidden email]
> Subject: Transform stack using eigenvectors
>
> Hi all
>
> I have written a plugin that calculates the 3D moments of inertia of a
> bone imaged in CT, using the Jama package and eigen decomposition.
>
> Now that I have a set of eigenvalues and eigenvectors I want to align my
> object such that its principal axes are parallel with the x,y,z axes of
> a new stack.
>
> Can anyone shed some light as to how to best achieve this?
>
> Mike
>
>  

[Stephan Saalfeld]

Hi Mike,

once you have the axes x'=(x1,y1,z1), y'=(x2,y2,z2), z'=(x3,y3,z3) of
the new coordinate frame you can easily write down the linear component
of the homogeneous matrix that transfers your coordinates into this
coordinate frame as:

 x' y' z'

 |  |  |
 v  v  v

 x1 x2 x3
 y1 y2 y3  } = A
 z1 z2 z3

such that Aa = a'.  If you want to map an image into the `normalized'
coordinate space, you will typically iterate over the pixels in the
target image and pick the corresponding coordinates from the
source---for this you will have to invert A.

I have in mind, that eigenvectors are only defined up to a rotation of
180deg such that you would have to find out where is top and where is
bottom for each of them.  I would be glad if someone could prove me
wrong with this, particularly because I don't find any reference about
it any more...

Best regards,
Stephan


On Thu, 2008-08-21 at 16:19 +0100, Michael Doube wrote:

> To answer my own question:
>
> The aligned position is the 3 dot products of the original coordinate
> (x,y,z) and the 3 eigenvectors (the principal axes).  Each eigenvector
> is a normal to the plane that contains the other 2 eigenvectors because
> they are orthogonal.  A plane is defined by the vector of its normal and
> the position of the plane on that normal.  The dot product of a plane
> and a point is the distance between plane and point, assuming that
> (0,0,0) is a solution for the plane.  So if you subtract the centroid
> from the point and work out the 3 dot products corresponding to the 3
> eigenplanes, you have the distance along each of the eigenvectors which
> defines the point in the principal axis coordinate frame, and hence
> gives you 'aligned' coordinates.  Then the aligned coordinates can be
> drawn in ordinary image coordinates and the object is aligned according
> to its principal axes.
>
> If that makes no sense or is blatantly wrong, please email me, otherwise
> it's going in a plugin tonight...
>
> Mike
> > -----Original Message-----
> > From: ImageJ Interest Group [mailto:[hidden email]] On Behalf Of
> > Michael Doube
> > Sent: Friday, August 15, 2008 5:57 PM
> > To: [hidden email]
> > Subject: Transform stack using eigenvectors
> >
> > Hi all
> >
> > I have written a plugin that calculates the 3D moments of inertia of a
> > bone imaged in CT, using the Jama package and eigen decomposition.
> >
> > Now that I have a set of eigenvalues and eigenvectors I want to align my
> > object such that its principal axes are parallel with the x,y,z axes of
> > a new stack.
> >
> > Can anyone shed some light as to how to best achieve this?
> >
> > Mike
> >
> >  

[Michael Doube]

Hi all

My implementation of this is here:
http://doube.org/plugins.html

It's set up for my particular situation: 16-bit, HU calibrated stacks of
CT images (dry bones in air), in which the specimen is already pretty
much aligned with its long axis parallel to the stack's Z axis (i.e.
most rotation is in the x-y plane).  Only nearest neighbour
interpolation is used, and the output stack has the same dimensions as
the input stack: data falling outside those dimensions are clipped - so
there's lots of room for improvement!

Cheers

Mike

Stephan Saalfeld wrote:

> Hi Mike,
>
> once you have the axes x'=(x1,y1,z1), y'=(x2,y2,z2), z'=(x3,y3,z3) of
> the new coordinate frame you can easily write down the linear component
> of the homogeneous matrix that transfers your coordinates into this
> coordinate frame as:
>
>  x' y' z'
>
>  |  |  |
>  v  v  v
>
>  x1 x2 x3
>  y1 y2 y3  } = A
>  z1 z2 z3
>
> such that Aa = a'.  If you want to map an image into the `normalized'
> coordinate space, you will typically iterate over the pixels in the
> target image and pick the corresponding coordinates from the
> source---for this you will have to invert A.
>
> I have in mind, that eigenvectors are only defined up to a rotation of
> 180deg such that you would have to find out where is top and where is
> bottom for each of them.  I would be glad if someone could prove me
> wrong with this, particularly because I don't find any reference about
> it any more...
>
> Best regards,
> Stephan
>
>
> On Thu, 2008-08-21 at 16:19 +0100, Michael Doube wrote:
>  
>> To answer my own question:
>>
>> The aligned position is the 3 dot products of the original coordinate
>> (x,y,z) and the 3 eigenvectors (the principal axes).  Each eigenvector
>> is a normal to the plane that contains the other 2 eigenvectors because
>> they are orthogonal.  A plane is defined by the vector of its normal and
>> the position of the plane on that normal.  The dot product of a plane
>> and a point is the distance between plane and point, assuming that
>> (0,0,0) is a solution for the plane.  So if you subtract the centroid
>> from the point and work out the 3 dot products corresponding to the 3
>> eigenplanes, you have the distance along each of the eigenvectors which
>> defines the point in the principal axis coordinate frame, and hence
>> gives you 'aligned' coordinates.  Then the aligned coordinates can be
>> drawn in ordinary image coordinates and the object is aligned according
>> to its principal axes.
>>
>> If that makes no sense or is blatantly wrong, please email me, otherwise
>> it's going in a plugin tonight...
>>
>> Mike
>>    
>>> -----Original Message-----
>>> From: ImageJ Interest Group [mailto:[hidden email]] On Behalf Of
>>> Michael Doube
>>> Sent: Friday, August 15, 2008 5:57 PM
>>> To: [hidden email]
>>> Subject: Transform stack using eigenvectors
>>>
>>> Hi all
>>>
>>> I have written a plugin that calculates the 3D moments of inertia of a
>>> bone imaged in CT, using the Jama package and eigen decomposition.
>>>
>>> Now that I have a set of eigenvalues and eigenvectors I want to align my
>>> object such that its principal axes are parallel with the x,y,z axes of
>>> a new stack.
>>>
>>> Can anyone shed some light as to how to best achieve this?
>>>
>>> Mike
>>>
>>>  
>>>      

[Thomas Boudier]

Hi Michael,

I've looked at your implementation of 3D moments calculation, I must
admit that I do not fully understand the Icxx, Icyy and Iczz calculation
especially the term you add :
 
                            // OK
                            Icxx += Math.pow(y * vH - centY3D, 2) +
Math.pow(z * vD - centZ3D, 2);
                            Icyy += Math.pow(x * vW - centX3D, 2) +
Math.pow(z * vD - centZ3D, 2);
                            Iczz += Math.pow(y * vH - centY3D, 2) +
Math.pow(x * vW - centX3D, 2);
                            Icxy += (x * vW - centX3D) * (y * vH - centY3D);
                            Icxz += (x * vW - centX3D) * (z * vD - centZ3D);
                            Icyz += (y * vH - centY3D) * (z * vD - centZ3D);

                            // ??
                            double shapeVolume = cstack * vW * vH * vD;
                            Icxx += shapeVolume * (vH * vH + vD * vD) / 12;
                            Icyy += shapeVolume * (vW * vW + vD * vD) / 12;
                            Iczz += shapeVolume * (vH * vH + vW * vW) / 12;

Is it because the computation is done relatively to the barycenter or
something like that (I must admit I'm not fully familiar with inertia
tensor) ? and why  /12 ?

thanks

Thomas

PS : By the way nice plugin  ;-)

> Hi all
>
> My implementation of this is here:
> http://doube.org/plugins.html
>
> It's set up for my particular situation: 16-bit, HU calibrated stacks
> of CT images (dry bones in air), in which the specimen is already
> pretty much aligned with its long axis parallel to the stack's Z axis
> (i.e. most rotation is in the x-y plane).  Only nearest neighbour
> interpolation is used, and the output stack has the same dimensions as
> the input stack: data falling outside those dimensions are clipped -
> so there's lots of room for improvement!
>
> Cheers
>
> Mike
>
> Stephan Saalfeld wrote:
>> Hi Mike,
>>
>> once you have the axes x'=(x1,y1,z1), y'=(x2,y2,z2), z'=(x3,y3,z3) of
>> the new coordinate frame you can easily write down the linear component
>> of the homogeneous matrix that transfers your coordinates into this
>> coordinate frame as:
>>
>>  x' y' z'
>>
>>  |  |  |
>>  v  v  v
>>
>>  x1 x2 x3
>>  y1 y2 y3  } = A
>>  z1 z2 z3
>>
>> such that Aa = a'.  If you want to map an image into the `normalized'
>> coordinate space, you will typically iterate over the pixels in the
>> target image and pick the corresponding coordinates from the
>> source---for this you will have to invert A.
>>
>> I have in mind, that eigenvectors are only defined up to a rotation of
>> 180deg such that you would have to find out where is top and where is
>> bottom for each of them.  I would be glad if someone could prove me
>> wrong with this, particularly because I don't find any reference about
>> it any more...
>>
>> Best regards,
>> Stephan
>>
>>
>> On Thu, 2008-08-21 at 16:19 +0100, Michael Doube wrote:
>>  
>>> To answer my own question:
>>>
>>> The aligned position is the 3 dot products of the original
>>> coordinate (x,y,z) and the 3 eigenvectors (the principal axes).  
>>> Each eigenvector is a normal to the plane that contains the other 2
>>> eigenvectors because they are orthogonal.  A plane is defined by the
>>> vector of its normal and the position of the plane on that normal.  
>>> The dot product of a plane and a point is the distance between plane
>>> and point, assuming that (0,0,0) is a solution for the plane.  So if
>>> you subtract the centroid from the point and work out the 3 dot
>>> products corresponding to the 3 eigenplanes, you have the distance
>>> along each of the eigenvectors which defines the point in the
>>> principal axis coordinate frame, and hence gives you 'aligned'
>>> coordinates.  Then the aligned coordinates can be drawn in ordinary
>>> image coordinates and the object is aligned according to its
>>> principal axes.
>>>
>>> If that makes no sense or is blatantly wrong, please email me,
>>> otherwise it's going in a plugin tonight...
>>>
>>> Mike
>>>    
>>>> -----Original Message-----
>>>> From: ImageJ Interest Group [mailto:[hidden email]] On Behalf Of
>>>> Michael Doube
>>>> Sent: Friday, August 15, 2008 5:57 PM
>>>> To: [hidden email]
>>>> Subject: Transform stack using eigenvectors
>>>>
>>>> Hi all
>>>>
>>>> I have written a plugin that calculates the 3D moments of inertia
>>>> of a bone imaged in CT, using the Jama package and eigen
>>>> decomposition.
>>>>
>>>> Now that I have a set of eigenvalues and eigenvectors I want to
>>>> align my object such that its principal axes are parallel with the
>>>> x,y,z axes of a new stack.
>>>>
>>>> Can anyone shed some light as to how to best achieve this?
>>>>
>>>> Mike
>>>>
>>>>        
>
>

--
/*****************************************************/
   Thomas Boudier, MCU Université Paris 6,
   UMR 7101 / IFR 83. Bat A 328, Jussieu.
   Tel : [phone]  Fax : [phone]
/****************************************************/

[Michael Doube]

Hi Thomas

It's not immediately apparent when you read the formal definition of
moments of inertia that you have to consider each voxel as a
parallelepiped with its own volume, mass, density and hence inertia -
usually the definition treats objects as integrals (infinite sums of
infinitely small volumes), or treats the masses as 'point masses'.  The
implication for the pixelated image is that you have to treat each pixel
(or voxel) as a finite element with a non-0 area or volume, so it is no
longer a (very) long sum of infinitesimal volumes, but a much shorter
sum of a finite number of real volumes.

An alternative way to view this is to consider the moment of the whole
to be the sum of the moments of the elements plus their moment relative
to the centroid of the whole, which are being calculated with the
parallel axis theorem.

So the moments of inertia are the sums of all the distances^2 *plus* the
sums of the moments of inertia of the parallelepipeds, which are the
(voxel's width^2 * the voxels height^2) / 12 for rotation around the
voxel's z-axis, as an example.

http://en.wikipedia.org/wiki/List_of_moments_of_inertia

A similar adjustment is necessary in 2D; the borrowed source of my
original moment of area macro did not have an adjustment for the area of
the pixels and assumed infinitesimally small pixels, so the result was
always a bit small when I tested it against a known area e.g. a filled
rectangle (which should evaluate to give moments of w*h^3 / 12 and h*w^3
/ 12).

Why /12?  I don't know the proof I'm sorry (can I play the "I'm just a
veterinarian" card now please?).

Anyway, rerunning that macro with the correct calculation is the current
task at hand so I'd better get to it.

Hope that helps

Mike

Thomas Boudier wrote:

> Hi Michael,
>
> I've looked at your implementation of 3D moments calculation, I must
> admit that I do not fully understand the Icxx, Icyy and Iczz calculation
> especially the term you add :
>
>                            // OK
>                            Icxx += Math.pow(y * vH - centY3D, 2) +
> Math.pow(z * vD - centZ3D, 2);
>                            Icyy += Math.pow(x * vW - centX3D, 2) +
> Math.pow(z * vD - centZ3D, 2);
>                            Iczz += Math.pow(y * vH - centY3D, 2) +
> Math.pow(x * vW - centX3D, 2);
>                            Icxy += (x * vW - centX3D) * (y * vH - centY3D);
>                            Icxz += (x * vW - centX3D) * (z * vD - centZ3D);
>                            Icyz += (y * vH - centY3D) * (z * vD - centZ3D);
>
>                            // ??
>                            double shapeVolume = cstack * vW * vH * vD;
>                            Icxx += shapeVolume * (vH * vH + vD * vD) / 12;
>                            Icyy += shapeVolume * (vW * vW + vD * vD) / 12;
>                            Iczz += shapeVolume * (vH * vH + vW * vW) / 12;
>
> Is it because the computation is done relatively to the barycenter or
> something like that (I must admit I'm not fully familiar with inertia
> tensor) ? and why  /12 ?
>
> thanks
>
> Thomas
>
> PS : By the way nice plugin  ;-)
>
>> Hi all
>>
>> My implementation of this is here:
>> http://doube.org/plugins.html
>>
>> It's set up for my particular situation: 16-bit, HU calibrated stacks
>> of CT images (dry bones in air), in which the specimen is already
>> pretty much aligned with its long axis parallel to the stack's Z axis
>> (i.e. most rotation is in the x-y plane).  Only nearest neighbour
>> interpolation is used, and the output stack has the same dimensions as
>> the input stack: data falling outside those dimensions are clipped -
>> so there's lots of room for improvement!
>>
>> Cheers
>>
>> Mike
>>