Posted by
Michael Doube on
Sep 03, 2008; 3:10pm
URL: http://imagej.273.s1.nabble.com/Transform-stack-using-eigenvectors-tp3695184p3695189.html
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_inertiaA 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
>>