Anisotropy Eigenvectors and eigenvalues

Dear Michael,

I have a question regarding the output from the Anisotropy function from the BoneJ software. Forgive me if it seems basic but I want to double check the format and significance of the eigenvector and eigenvalue output.

For example, for one of my bone cubes I get the following output:

Fabric tensor vectors
||0.264| 0.288|-0.920||
||0.122|-0.957|-0.265||
||0.957| 0.042| 0.288||

Fabric tensor values
||0.010|0.000|0.000||
||0.000|0.024|0.000||
||0.000|0.000|0.047||

Image	DA    tDA   V1,1  V1,2   V1,3    V2,1     V2,2   V2,3    V3,1   V3,2  V3,3   D1     D2	     D3
cub102   0.786 4.676	0.264 0.288  -0.920  0.122	-0.957  -0.265  0.957	0.042 0.288 0.010 0.024 0.047

Regarding the eigenvectors, am I right in assuming that the vectors are arranged by column. Thus, in this example, the first fabric tensor vector direction is defined by the x,y,z coordinates V1.1 0.264, V2.1 0.122 V3.1 0.957 (as opposed to being defined by the values moving across each row e.g. V1.1 V1.2 V1.3).

Following on from this, am I right in assuming that each column i.e. eigenvector (or row if I am wrong!) corresponds to its associated eigenvalue i.e. eigenvector/column 1 is associated to the first eigenvalue, column 2 to the second eigenvalue and so on?

Lastly, regarding the eigenvalues, is the largest eigenvalue (0.047 in this instance) associated with the main principal direction i.e. the main trabecular direction? Or is the largest eigenvalue telling me that there are more boundaries (and therefore smaller intercept lengths) in this direction and thus it actually relates to the shortest principal direction?

Thanks for your time and for this brilliant programme!

JD

Hi @pra04dpc

Good observation. Checking the code indicates that your suspicion is correct. Each eigenvector is represented as a column in the log output. The corresponding eigenvalue is in the same column in the fabric tensor value matrix (on the diagonal).

In the Results table output, the eigenvalues are listed as v(row, column). If you find this counter-intuitive, please let me know and I’ll fix it to v(column, row), and update the documentation - which I see does not have detail on this option at the moment. Either way, I have to write some docs!

Update 2016-01-15: Edited docs to inform on the current status on this point.

As for interpretation: have a look at the Harrigan & Mann paper. If I remember it correctly, the smaller the eigenvalue, the longer the mean intercept length in that direction (because eigenvalue D relates to ellipsoid radius r as r=1/sqrt(D), and the radius relates to mean spacing in that principal direction)), which means fewer boundaries per unit of length (or inversely, the greater the spacing of repeating units, or put another way, the lower the spatial frequency in that direction). You can try this out with test images, by making a simple image with alternating black and white pixels that you space differently in x, y and z; (eg planes in xy, white every 3rd plane, and planes in yz, white every 10th plane, and planes in xz, white every 20th plane, for example).

2 Likes

Hallo,

I have a question related to this topic, and I did not want to create a new item. I would be thankful if you help me with that.

My question is about the ellipsoid representation of the fabric tensor which I obtain from the bonej-anisotropy plugin.

I build the fabric tensor by summation of dyadic products of the
eigenvalues and eigenvectors** F** = sum (fiFi) = sum (fi(fi ⊗ fi)),
i=1,2,3 where F is fabric tensor defined by eigenvalues fi and
eigenvectors fi.

I fitted an ellipsoid to the fabric tensor, and I expect that I should
get same results as I get from bonej. However, it seems the ellipsoid I obtained directly from the fabric tensor is
rotated and is not aligned with the x-y-z coordinate while I the fit
from bonej is properly aligned (the misalignment degree is about 80). I checked that with particle analyzer too, and I attached its results to this post.

In the log file, the rotation matrix is also reported, I am wondering if I consider should it in the fabric tensor. As you have answered to
the similar question before in [here],
I think we do not to need any rotation matrix because the eigenvectors
show the orientation of the ellipsoid, right?

I have another question, too. Is it possible to print MIL values directly from anisotropy plugin?

Thanks for your time,
Best regards,
Mohammad

1 Like

Dear Mohammed,

The fabric tensor is the same as the 3×3 matrix that defines the ellipsoid used to fit the mean intercept coordinates. Have a look at the the code, where you will see that the coefficients of the ellipsoid equation are used to create a new matrix that is Eigen decomposed to make the fabric tensor. The ellipsoid is centred at 0 so is the top left 3×3 matrix of the more general 4×4 quadric matrix, as described nicely here.

The red-green image that you include is not a fitted ellipsoid - it’s the product of mean intercept length in each direction and the unit vectors. These are the coordinates to which the ellipsoid is fitted.

1 Like

Hello Michael,

Thank you for such a useful tool! I have a couple more specific questions on the same topic as the original post, because I need to know exactly what the results correspond to:

1,) The MIL tensor that is calculated by the code…is that the same as the MIL tensor M as shown in the Harrigan and Mann paper? In other words, the components of the MIL tensor (which I understand are not displayed to the user) are the coefficients of the ellipsoid coefficients as published in the paper? (Unfortunately, I do not have enough background in Java to sort through the code that you linked to in one of your responses above.)

2.) What tensor exactly do the “Fabric tensor vectors” and “fabric tensor values”, which are displayed in the results, correspond to? IF you are calculating the MIL tensor M as in Harrigan and Mann, are the “fabric tensor vectors” just the eigenvectors of M? And the “fabric tensor values” just the eigenvalues of M? Where I am confused is that you only cite Harrigan and Mann, but if I remember correctly, nowhere in that paper do they discuss “fabric tensors”. If my understanding is correct, fabric tensors have been calculated by later authors (such as Cowin and Zysset) as the inverse square root of the MIL tensor. For example, if I calculate the MIL tensor M as is Harrigan and Mann, then Zysset defines a fabric tensor (let’s call it F for convenience, because his notation is different) such that F = M^(-1/2). And he further normalizes such that trace(F) = 3.
I just want to be clear…this is NOT what is happening in the BoneJ code?

3.) I want to make sure I understand the convergence process clearly. A point is picked at random somewhere in the sample…but the (x,y,z) coordinates of this point must be sufficiently far away from the boundaries to allow for a vector which is the size of the “radius”, shown in the setup GUI, to not exceed the image bounds? So in other words, if I have an image of dimensions x=A, y= B, and z=C, and the “radius” of the analysis is r, then each point that the program picks as the center of the MIL analysis (center of the sphere) must fall in the (x,y,z) space:
0+r<x<A-r 0+r<y<B-r 0+r<z<C-r

In other words, I just want to make sure that the analysis is not sampling empty space outside of the image. After a point is picked, MIL analysis is run. Then, the process is repeated until results converge below some tolerance?

Thanks for the help!

-Steve

1 Like

Yes, they are. All the code is doing is creating a tensor from the ellipsoid coefficients, then eigendecomposing it.

Yes, exactly.

Cowin just says “fabric tensor is generalized to mean any symmetric second rank tensor that
characterizes the local geometric arrangement of solid material or microstructure of a porous material”, so BoneJ’s usage is not wrong.

DA is calculated according to the standard in the bone field (I actually followed the Skyscan CTAn technical manual), as DA = 1 - minEigenvalue / maxEigenvalue.

That’s correct. The bit of code that does the work is not too code-y:

//get the equation coefficients
double[] coEf = (double[]) ellipsoid[3]; 
    //pack the coefficients into a tensor
		double[][] tensor = { { coEf[0], coEf[3], coEf[4] }, 
				{ coEf[3], coEf[1], coEf[5] }, { coEf[4], coEf[5], coEf[2] } };
//turn it into a Matrix object to do matrix-y stuff
		Matrix M = new Matrix(tensor);
//do the eigendecomposition
		EigenvalueDecomposition E = M.eig();
//get the eigenvalues
		Matrix EigenVal = E.getD();
		double[] diag = EigenVal.diag().getColumnPackedCopy();
//calculate degree of anisotropy
		da = 1 - diag[0] / diag[2];

That’s correct. There’s a buffer zone of one sphere radius for sampling sphere centres so that sampling probes don’t go outside the image bounds. Then sampling spheres are seeded at random points, the DA calculated for each seed point and further spheres get seeded until the coefficient of variation in the last 100 DA measurements drops below 0.005. I’m not that happy with this approach, but it was the best I could come up with at the time. The main point is that fabric orientation passes through all points in all directions, but a single sampling sphere samples most points in only one direction.

Thanks for this discussion and I’m more than happy to incorporate alternative approaches to BoneJ’s anisotropy plugin. I’ve never been that convinced by DA / MIL, perhaps due to this particular implementation and difficulties in validating it properly.

1 Like

A couple of comments/questions regarding your response, above.

1.) The BoneJ MIL analysis seems to be more robust than the CTAn implementation for two reasons:
A) In BoneJ, the radius of the sampling spheres is set so that the sphere is completely enclosed within the image stack. In CTAn, I believe the radius is set so that the sphere encloses the entire image stack. So…if your stack is not cubical but is instead something more like a thin rectangular panel, then CTAn’s results will be biased (and misleading).
B.) BoneJ can sample hundreds of spheres from a given stack to give the final results. As far as I am aware, to do this in CTAn, one would have to create each individual spherical VOI and then run the analysis on it.

2.) I found an earlier post (06/19/14) from the Google groups in which you asked:
Would it be helpful to have an option to visualise the fabric tensor alongside the original image in the 3D viewer, something like is done for Moments of Inertia?
To me, this would be extremely helpful. I am trying to visualise the angle between the dominant direction (eigenvector with the smallest eigenvalue) and the Z direction of my stack.

3.) The validity/usefulness of the MIL analysis seems to be dependent on at least two things:
(A) the size of the sampling sphere with respect to the spatial frequency of features. For example, if my sampling sphere is only 50 pixels wide and my trabecular spacing is on average about 30 pixels, then I will probably not have enough “boundary crossings” (which is what MIL measures) for the results to be significant. In other words, the number of trabecular spacings within the sampling sphere seems to be important.
(B) the correlation coefficients for the ellipsoidal fit.

CTAn provides statistics about (A) and (B) in its anisotropy results. For (A), it outputs the average number of intersections. For (B) it does something with the correlation coefficients in the calculation of confidence intervals and t-values for the eigenvalues. (If you have access to CTAn’s method notes, this is described in detail in Method Note 31: Anisotropy, MIL, and Stereology.)

This is currently missing in BoneJ. Any way that some statistics on (A) or (B) could be included in the results? Or do you have other ideas as to how to assess the validity of the MIL analysis?

Again, I appreciate the help!

1 Like

Hi Steve and Michael,

This discussion was helpful for me, but I still have ambiguities in the representation of the fabric tensor.

Michael, as you have mentioned in your post, the output fabric tensor calculated by bonej shows the representation of the ellipsoid, right? (based on (Harrigan and Mann, 1984))
I tried to plot this ellipsoid.

Subsequently, I obtained MIL from CTan software and build the fabric tensor. Based on the following formula:
M = ∑ = m i M i = ∑ = m i (m i ⨂m i )

I plotted the representative ellipsoid related to the fabric tensor M.

Is this approach correct, and if so, why these two ellipsoids are not oriented in the same way?

Thank you very much for your help!
Mohammad

1 Like

The MIL results could be completely different between CTan and BoneJ. In other words, for a given VOI, the two programs could end up showing two completely different orientations. This is, in part, due to the issues I described in my post above. CTan creates a single sphere and uses it for MIL, whereas BoneJ runs a convergence routine- usually going through hundreds of spheres (hundreds of MIL analyses) before outputting results (unless you selected “single sphere” option in BoneJ)?

You would need to share more information on what you ran. The VOI was exactly the same between BoneJ and CTan? What was the dimension of the VOI? Was it spherical, cubical, or something else? Did you just use the default parameters in BoneJ Anisotropy or something custom? What were the Stereology settings in CTan (test line spacing and number of orientations)?

Perhaps the quickest way to see what is going on would be to look at the actual vector results from BoneJ/CTan. If you want to post them here, it would be interesting to see.

-Steve

1 Like

Thanks for your response Steve.

I have used similar VOI and images in both software. I made the binary image in ImageJ and used it in bonej and CTan. The dimension of the VOI were: 10.69 x 10.69 x 18.7 mm^3 with the resolution of 20 μm.
The images are taken from the cylindrical foam aluminum samples. I used default parameters in Bonej Anisotropy as I could not find proper custom ones. Regarding the Stereology setting, I used the default definition and did not change them.

Here are the bonej and CTan vectors:

Bonej

#DA	v1,1	v1,2	v1,3	v2,1	v2,2	v2,3	v3,1	v3,2	v3,3	D1	D2	D3
#0.163	-0.236	0.459	0.857	-0.123	0.86	-0.495	0.964	0.222	0.147	4.465	4.735	5.332

I rebuild the vectors as:

m1 = np.array([-0.236, -0.123,  0.964])
m2 = np.array([ 0.459,  0.86 ,  0.222])
m3 = np.array([ 0.857, -0.495,  0.147])

M1 = 4.465
M2 = 4.735
M3 = 5.332
F = np.outer(m1,m1) * M1 + np.outer(m2,m2) * M2 + np.outer(m3,m3) * M3     #fabric tensor

CTan

Eigenvectors:,,,
E-vector 1:,,0.00559, -0.07066, 0.99748,
E-vector 2:,,0.75220, 0.65757, 0.04237,
E-vector 3:,,-0.65891, 0.75008, 0.05683,
Eigenvalue 1,,0.05839,
Eigenvalue 2,,0.07918,
Eigenvalue 3,,0.08398,

As I wanted to have the longest axis in 3rd direction, I have sorted the result to have M3 > M2 > M1.

m1 = np.array ([-0.65891, 0.75008,  0.05683])
m2 = np.array ([ 0.7522 , 0.65757,  0.04237])
m3 = np.array ([ 0.00559,-0.07066,  0.99748])

M3 = 0.3713972898
M2 = 0.3189203985
M1 = 0.3096823118

F2 = np.outer(m1,m1) * M1 + np.outer(m2,m2) * M2 + np.outer(m3,m3) * M3   #fabric tensor

and I have plotted the ellipsoid fitted to F1 and F2, which you can see in my previous post.

But looking at the results, now I find out what my problem is. The different is coming from the way I have sorted the eigenvectors and eigenvalues if you agree. Therefore, I should change my question. Can I do the similar sorting for the results I obtained from bonej? In other words, from the type of my image, I expect that to have an ellipsoid with longer diameter align to the z-direction. Is such a conclusion right?

Thanks for your help.
Mohammad

1 Like

Two things to point out that I can see right now:

  1. Your VOI is not an isotropic size. In other words, it is not a cube or a sphere. If you presented to CTan that VOI of dimensions 10.69x10.69x18.7, the results are going to be wrong (at least based on my understanding!) CTan will create a sphere big enough to accommodate all 3 dimensions of your VOI, which means that part of the sphere will be empty space outside of your original VOI. Therefore, the CTan results are going to have a pre-bias in one direction (the dimension corresponding to the length of 18.7).
    (see my previous post. This is one of the advantages of BoneJ because it will run MIL on spheres that are smaller than your VOI.)

2.) Looks like you are trying to build the MIL tensor from the eigenvalues and eigenvectors. From what I can tell, looks like you have it right for BoneJ. But for CTan, you are taking the “principal MILs” to be the eigenvalues of the MIL tensor, but they are not. In the CTan output file, the eigenvalues of the MIL tensor are labelled as “Eigenvalue 1”, “Eigenvalue 2”, and “Eigenvalue 3”. (If I remember correctly, the “principal MILs” are actually the inverse square root of the eigenvalues.)

Surprisingly, the results don’t look that much different. CTan says that the dominant direction is roughly the z dimension, since E-vector 1 is roughly (0,0,1). This is about the same for BoneJ. So if you take into account point 2.), you might see the results get closer.

Any word from Michael on these last three posts???

2 Likes

See the second post in this thread. The eigenvectors are reported as v(row, column), so your matrix is not correct.

Yes, a goodness of fit measurement would be helpful. If you have a good reference to use as guidance to make the calculation I’d gladly implement it.

1 Like

Yes that’s important. BoneJ gets around this issue by summing all the boundary crossings between sampling spheres, so the counts in each direction continue to grow as you add more spheres. That’s one reason why the BoneJ DA result stabilises over time: the number of boundary hits on each unit vector continue to grow until a stable result is achieved. A better design would be to do this without sampling spheres and to run grids of vectors through the image, and to increase the density of the vectors / number of directions iteratively until the result converges. This would also avoid issues of uneven sampling around the edges of images (the corners never get sampled in the current scheme).

1 Like

Michael,

Based on Anisotropy manual in Bonej I build the matrix. It says "Eigenvalue and eigenvector matrices as column vectors in the log window."
From your answer to “pra04dpc” I understood that eigenvectors are sorted as m1=[v11,v21,v31] m2=[v12,v22,v32] m3=[v13,v23,v33].
Are you sure with your comment?

Steve,

I edited my post and replaced the eigenvalues instead of the principal MIL values. Thanks for your comment, but I still do not understand your observation. You mentioned: “CTan says that the dominant direction is roughly the z dimension since E-vector 1 is roughly (0,0,1). This is about the same for BoneJ.”

E-vector 1 has the smallest eigenvalue in both software. How can we conclude that it is the dominant direction?
It arises my previous question again, are the output eigenvectors already sorted in CTan and bonej? Because as I have mentioned, I resort the eigenvector as you see the results of CTan, and can I do the same for the results of bonej?
Many Thanks for your kind attention,

You are correct, my misreading of your vectors as untransposed.

In BoneJ they are. The Eigendecomposition is performed by the JAMA library, so the resulting matrices are already sorted so that the columns of the eigenvector and eigenvalue matrices correspond.

1 Like

Thanks, Michael,
This conversation was very helpful for me. I replot the ellipsoid that is fitted on the fabric tensor calculated by bonej. It looks similar to one calculated by CTan, and I edit my first post.

I have last question sorry for the disturbance.
To calculate the angle of the largest eigenvalue with respect to z-axis, we can use acos(m3.e3) where m3 is eigenvector and e3=[0,0,1], right? Based on my results, it gives me the wrong answer. I mean, in my case I have m3 = np.array([ 0.857, -0.495, 0.147]) which gives angle of 81.5.
Given 1,2,3 correspond to x-y-z in the output, I doubt if I am mistaken with this calculation. Could you please help me with that?
Thank you so much!

1 Like

Hi @MoG,

I played around a bit with the BoneJ code to help visualise your issues and come up with the following observations / confirmations:

  1. The MIL point cloud in 3D red/green that BoneJ can show is really well approximated by an ellipsoid (in blue points in the image) with radii calculated as the inverse square root of the eigenvalues and oriented to the eigenvectors of the fabric tensor.

  2. Adding the input image to the same 3D viewer shows that the long axis of the MIL point cloud is oriented parallel to the plates of a piece of trabecular bone, which makes sense because this is the direction that line probes have to travel furthest to hit the next boundary.

As for your last issue, yes the 3rd element of each eigenvector should correspond to the z component of a unit vector, so you can form a triangle with the hypotenuse = 1, z-axis length = 3rd vector element, then calculate the angle between them using the relation cos(theta) = adjacent / hypotenuse rearranged to cos^-1(adjacent / hypotenuse) = acos(3rd element of eigenvector) to work out the angle between the vector and the z axis. What is the ‘correct answer’ you are expecting? Can you get it from the third element of one of the other eigenvectors? Bear in mind that

  1. Eigen matrices are usually sorted in size order - JAMA does it in ascending order so the largest eigenvalue is at (3,3). You are already doing this.

  2. The largest eigenvalue corresponds to the shortest ellipsoid radius (in real units of length), hence highest spatial frequency, because ellipsoid radius = 1/sqrt(eigenvalue)

2 Likes

Many appreciations, Michael,

True, so I should calculate the angle of m1 which corresponds to the largest diameter/ radius of the ellipsoid. With this, the angle would be 15.4 (for bonej) and 4.06 (for CTan), and it is what I would expect by looking at the images.
Based on our previous conversation, I use the results of the bonej.