theta angle - Githubissues

zilinyan commented 5 years ago

Hi @dlegland , thanks for the codes. I want to measure 3D ellipsoid particle orientation in 3D space. To test your code, I made around 2000 ellipsoids digitally which are randomly distributed in space. I sorted their theta angles (-90-90) with a bin size of 5 degrees. It seems most of the particles have a -45 and 45-degree theta angle. Theoretical calculation suggests that particles that around theta=0 are the most, as reported in this figure. I also tried real randomly packed particles which are obtained from X-ray CT, the results are similar.

I am wondering if these results have something to do with these lines in your codes? thanks

// extract |cos(theta)| Matrix mat = svd.getU(); double tmp = hypot(mat.get(0, 0), mat.get(1, 0)); double phi, theta, psi;

        // avoid dividing by 0
        if (tmp > 16 * Double.MIN_VALUE) 
        {
            // normal case: theta <> 0
            psi     = atan2( mat.get(2, 1), mat.get(2, 2));
            theta   = atan2(-mat.get(2, 0), tmp);
            phi     = atan2( mat.get(1, 0), mat.get(0, 0));
        }
        else 
        {
            // theta is around 0 
            psi     = atan2(-mat.get(1, 2), mat.get(1,1));
            theta   = atan2(-mat.get(2, 0), tmp);
            phi     = 0;
        }

dlegland commented 5 years ago

Hi,

thanks for this interesting observation.

The lines of code you refer too simply consist in converting 3D rotation matrix into three Euler angles. This code is rather common, and I'm surprised there is such a difference with theory.

There may be some discretization effects, but I suppose you considered ellipsoids "large enough"?

Have you noticed similar biaswith other library/software?

I will try to investigate this effect.

Best, David

dlegland commented 5 years ago

Hi,

I have made some additionnal checks, and did not exhibit any suspect behaviour.

However, I think I understood why you obtained such curve. It depends on the way you generate the orientation of ellipsoids. Let us suppose that ellipsoids are elongated, with equal second and third radius (i.e., "prolate" ellipsoids). Then, the orientation is solely determined by the (phi, theta) pair.

If the (phi,theta) pairs are distributed uniformly over the surface of the sphere, one should obtain the blue curve in your graph.
if you generate (phi,theta) pairs by choosing a 3D point inside a centered cube, and converting cartesian coordinates to spherical, then you over-sample points in the corners and around the edges of the cube. As the cube edges have angles around 45 degrees, this should cause the effect observed on the orange curve.

Can you confirm? If yes we can close.

zilinyan commented 5 years ago

Hi,

I have made some additionnal checks, and did not exhibit any suspect behaviour.

However, I think I understood why you obtained such curve. It depends on the way you generate the orientation of ellipsoids. Let us suppose that ellipsoids are elongated, with equal second and third radius (i.e., "prolate" ellipsoids). Then, the orientation is solely determined by the (phi, theta) pair.

If the (phi,theta) pairs are distributed uniformly over the surface of the sphere, one should obtain the blue curve in your graph.

if you generate (phi,theta) pairs by choosing a 3D point inside a centered cube, and converting cartesian coordinates to spherical, then you over-sample points in the corners and around the edges of the cube. As the cube edges have angles around 45 degrees, this should cause the effect observed on the orange curve.

Can you confirm? If yes we can close.

Hi @dlegland , Thanks for your prompt help. i do use "prolate" ellipsoids (a>b=c). But I also tried on natural particles and also caught a strange behavior at 45 degree. I cannot provide the details at the moment on how I generated the ellipsoid as I used a third-party software. However, I compared the result by another library, 3D Suite ImageJ plugin (http://imagejdocu.tudor.lu/doku.php?id=plugin:stacks:3d_ij_suite:start) with previously reported results. The 3D suite calculates the theta angle in the range (0, 90). In order to do the comparison, i replaced the negative theta in MorphlibJ with its absolute value (the orientation should be symmetric about XY plane? ).

Bellowing are the comparisons. I think 2000 particles are kind of sufficient albeit there are some deviations (gray line) compared to the theoretical values (orange line).

zilinyan commented 5 years ago

Hi,

thanks for this interesting observation.

The lines of code you refer too simply consist in converting 3D rotation matrix into three Euler angles. This code is rather common, and I'm surprised there is such a difference with theory.

There may be some discretization effects, but I suppose you considered ellipsoids "large enough"?

Have you noticed similar biaswith other library/software?

I will try to investigate this effect.

Best, David

Hi @dlegland ,

Could it be because of the different definitions of inertia of moments?

In your codes: https://github.com/ijpb/MorphoLibJ/blob/820628a09d3570275d5804069371e2a7fe03f237/src/main/java/inra/ijpb/measure/region3d/InertiaEllipsoid.java#L257

               // convert coordinates relative to centroid 
                double x2 = x * sx - moment.cx;
                double y2 = y * sy - moment.cy;
                double z2 = z * sz - moment.cz;

                // update coefficients of inertia matrix
                moment.Ixx += x2 * x2;
                moment.Iyy += y2 * y2;
                moment.Izz += z2 * z2;
                moment.Ixy += x2 * y2;
                moment.Ixz += x2 * z2;
                moment.Iyz += y2 * z2;

In 3D Suite https://github.com/mcib3d/mcib3d-core/blob/7edc6aa3fb35dcc3770e76d753351a0382c5da1b/src/main/java/mcib3d/geom/Object3DVoxels.java#L507

for (Voxel3D voxel : voxels) { vox = voxel; i = vox.getX(); j = vox.getY(); k = vox.getZ(); s200 += resXY resXY (j - by) (j - by) + resZ resZ (k - bz) (k - bz); s020 += resXY resXY (i - bx) (i - bx) + resZ resZ (k - bz) (k - bz); s002 += resXY resXY (i - bx) (i - bx) + resXY resXY (j - by) (j - by); s110 += resXY resXY (i - bx) (j - by); s101 += resXY resZ (i - bx) (k - bz); s011 += resXY resZ (j - by) * (k - bz);

                final double xvWcX = x * vW - cX;
                final double yvHcY = y * vH - cY;
                final double zvDcZ = z * vD - cZ;
                Icxx += (yvHcY * yvHcY + zvDcZ * zvDcZ + voxVhVd) * voxMass;
                Icyy += (xvWcX * xvWcX + zvDcZ * zvDcZ + voxVwVd) * voxMass;
                Iczz += (yvHcY * yvHcY + xvWcX * xvWcX + voxVhVw) * voxMass;
                Icxy += xvWcX * yvHcY * voxMass;
                Icxz += xvWcX * zvDcZ * voxMass;
                Icyz += yvHcY * zvDcZ * voxMass;

Best regards, Yan

dlegland commented 5 years ago

Well, different definitions may lead to different results, yes!

I used definition of statistical moments, widely described in image analysis books, and also implemented in Matlab. It seems 3D suite and BoneJ use some king of mechanical related defintion. I do not know how the two relate to each other...

For information, I attached a small report that gathers the references, and some investigations on the normalization coefficient for ellipsoid radius. inertiaEllipsoid.pdf

The explanations for conversion of angles can be obtained here: http://www.gregslabaugh.net/publications/euler.pdf

zilinyan commented 5 years ago

Hi, Thanks for your prompt reply and for sharing the math behind the codes. I still don't catch why we can use second order moments to resolve the problem. Here I found another reference which is also using the statistical moments you mentioned. They use different terms (https://uk.mathworks.com/matlabcentral/mlc-downloads/downloads/submissions/48913/versions/4/screenshot.jpg) which seem to the same as the Bonej and ImageSuite implementation except for the mass term. .

I am also interested to know how much the results are different by the two different treatments. It is very easy to adapt your current codes. Can you teach me how to compile your repo with ImageJ package? Thanks in advance! Yan

dlegland commented 5 years ago

Thanks for the link!

to compile the repo, we use Eclipse plateform, and the maven configuration file. When properly installed, using "File->import->maven->Existing maven project", and selecting the pom.xml file of MorphoLibJ should install all dependencies (well... mainly ImageJ itslef, and the JAMA matrix package).

David

dlegland commented 5 years ago

Hi, The latest version of MorphoLibJ (release v1.4.1) provides ellipsoid as "Equivalent Ellipsoid" rather than "Inertia Ellipsoid". This should help avoid problems. It is also possible to extract the components of the eigen vectors. Therefore it is possible to investigate results without necessarily re-compiling MorphoLibJ. User manual was updated as well.

Best, David

ijpb / MorphoLibJ

theta angle #30