Intersection of longest principal axis and object surface

Hello there,

I am trying to automatically extract bone length of segmented skeletal elements, without user intervention. While max. Feret’s diameter from the Particle analyser is a good approximation, it is not really a straight line along the center of the bone, which is what I ideally want to extract from the image. However, the longest principal axis that is used for the moments of inertia looks like a better approximation. Is it possible to get the coordinates of the two points at the intersections between the principal axis and the surface of the object? It would be where the red line touches the bone in the attached picture.

Thank you, I hope the question is clear.

Screen Shot 2021-02-19 at 17.06.42

The Moments plugin does something like this internally at least while calculating a minimal aligned bounding box to set up the rotated image stack.

Note though it’s not the intersections of the principal axes with the surface of the bone, but a caliper radius of the bone in the 3 orthogonal eigenvectors. Really irregular bones will give results you might not expect. However, you could get most of the way to the answer you want just by aligning your bone with Moments and then extracting the maximum distances from the centre of the stack in ± x, y, and z.

Hello Michael,

Thank you for your prompt reply. I understand what you suggest, which is a great idea. The real issue is that we are trying to develop a pipeline to extract the length of multiple bones (each with a different orientation) from a whole-body scan. Is there any way to recursively align every segmented object with Moments and extracting the maximum distances from the centre of the substack that contains it? Please see image of our ideal collection of objects.
Thanks again!
190426_DTA12_40um 3d view.tiff (678.3 KB)

Thanks for the image. In this case Particle Analyser is probably the best bet, augmented with a bit of code borrowed / inspired by the method posted above from Moments. It would use the inertia tensor to find a minimal aligned bounding box for each particle (but with a centre free to move from the geometric centroid). Would that work for you?

Yes! I think that’s exactly what I need. Not sure how to do it, though. I have pasted below what I think is the code to get the size of the bounding box, but it would need to start from the objects already aligned to the inertia tensors. Can this be done directly from the results of the Particle analyser? Apologies if this is obvious, I don’t have much experience with coding in ImageJ

/**
	 * Find side lengths in pixels of the smallest stack to fit the aligned image
	 *
	 * @param E Rotation matrix
	 * @param imp Source image
	 * @param centroid 3D centroid in 3-element array {x,y,z}
	 * @param startSlice first slice of source image
	 * @param endSlice last slice of source image
	 * @param min minimum threshold
	 * @param max maximum threshold
	 * @return Width, height and depth of a stack that will 'just fit' the aligned image
	 */
	private static int[] getRotatedSize(final Matrix E, final ImagePlus imp,
		final double[] centroid, final int startSlice, final int endSlice,
		final double min, final double max)
	{
		final ImageStack stack = imp.getImageStack();
		final Calibration cal = imp.getCalibration();
		final double xC = centroid[0];
		final double yC = centroid[1];
		final double zC = centroid[2];
		final Rectangle r = imp.getProcessor().getRoi();
		final int rW = r.x + r.width;
		final int rH = r.y + r.height;
		final int rX = r.x;
		final int rY = r.y;

		final double vW = cal.pixelWidth;
		final double vH = cal.pixelHeight;
		final double vD = cal.pixelDepth;

		final double[][] v = E.getArrayCopy();
		final double v00 = v[0][0];
		final double v10 = v[1][0];
		final double v20 = v[2][0];
		final double v01 = v[0][1];
		final double v11 = v[1][1];
		final double v21 = v[2][1];
		final double v02 = v[0][2];
		final double v12 = v[1][2];
		final double v22 = v[2][2];

		final int d = imp.getStackSize();
		double[] sliceXTmax = new double[d + 1];
		double[] sliceYTmax = new double[d + 1];
		double[] sliceZTmax = new double[d + 1];

		final AtomicInteger ai = new AtomicInteger(startSlice);
		final Thread[] threads = Multithreader.newThreads();
		for (int thread = 0; thread < threads.length; thread++) {
			threads[thread] = new Thread(() -> {
				for (int z = ai.getAndIncrement(); z <= endSlice; z = ai.getAndIncrement()) {
					IJ.showStatus("Getting aligned stack dimensions...");
					final ImageProcessor ip = stack.getProcessor(z);
					final double zCz = z * vD - zC;
					final double zCzv20 = zCz * v20;
					final double zCzv21 = zCz * v21;
					final double zCzv22 = zCz * v22;
					double xTmax = 0;
					double yTmax = 0;
					double zTmax = 0;
					for (int y = rY; y < rH; y++) {
						final double yCy = y * vH - yC;
						final double yCyv10 = yCy * v10;
						final double yCyv11 = yCy * v11;
						final double yCyv12 = yCy * v12;
						for (int x = rX; x < rW; x++) {
							final double pixel = ip.get(x, y);
							if (pixel < min || pixel > max) continue;

							// distance from centroid in
							// original coordinate system
							// xCx, yCx, zCx
							final double xCx = x * vW - xC;
							// now transform each coordinate
							// transformed coordinate is dot product of original
							// coordinates
							// and eigenvectors
							final double xT = xCx * v00 + yCyv10 + zCzv20;
							final double yT = xCx * v01 + yCyv11 + zCzv21;
							final double zT = xCx * v02 + yCyv12 + zCzv22;
							// keep the biggest value to find the greatest distance
							// in x, y and z
							xTmax = Math.max(xTmax, Math.abs(xT));
							yTmax = Math.max(yTmax, Math.abs(yT));
							zTmax = Math.max(zTmax, Math.abs(zT));
						}
					}
					sliceXTmax[z] = xTmax;
					sliceYTmax[z] = yTmax;
					sliceZTmax[z] = zTmax;
				}
			});
		}

The best bet is to add a new method or two to ParticleAnalysis, reusing where possible code that already exists and results that have already been calculated.

Keep an eye on this issue:

And this development branch:

@arosello
I have something like this working - please take a look at these images and let me know if it’s what you have in mind.

This is a synthetic image of 703 ellipsoids, all 40 × 16 × 8 pixels long×wide×deep, at various rotations. The yellow boxes are the axis-aligned measuring boxes:

The box length, width & depth all agree with the expected dimensions, with a bit of error due to discretisation artefact. Centroid is in total agreement with expectation. This is the test image stack: Tilted_Ellipsoids.zip (398.7 KB)

Numerical results look like this:

Label ID Vol. (pixels³) x Cent (pixels) y Cent (pixels) z Cent (pixels) Box x (pixels) Box y (pixels) Box z (pixels) Box l0 (pixels) Box l1 (pixels) Box l2 (pixels)
Tilted Ellipsoids.tif 1 2621 22 22 23 22.000 22.000 23 40 16.000 8.000
Tilted Ellipsoids.tif 2 2671 66 22 23 66.000 22.000 23 40 14.495 7.866
Tilted Ellipsoids.tif 3 2689 110 22 23 110.000 22.000 23 40 15.218 7.740
Tilted Ellipsoids.tif 4 2683 154 22 23 154.000 22.000 23 40 15.126 7.461
Tilted Ellipsoids.tif 5 2661 198 22 23 198.000 22.000 23 40 15.620 7.459
Tilted Ellipsoids.tif 6 2661 242 22 23 242.000 22.000 23 40 15.620 7.459
Tilted Ellipsoids.tif 7 2683 286 22 23 286.000 22.000 23 40 15.126 7.461
Tilted Ellipsoids.tif 8 2689 330 22 23 330.000 22.000 23 40 15.218 7.740
Tilted Ellipsoids.tif 9 2671 374 22 23 374.000 22.000 23 40 14.495 7.866
Tilted Ellipsoids.tif 10 2621 418 22 23 418.000 22.000 23 40 16.000 8.000
Tilted Ellipsoids.tif 11 2671 462 22 23 462.000 22.000 23 40 14.495 7.866
Tilted Ellipsoids.tif 12 2689 506 22 23 506.000 22.000 23 40 15.218 7.740
Tilted Ellipsoids.tif 13 2683 550 22 23 550.000 22.000 23 40 15.126 7.461
Tilted Ellipsoids.tif 14 2661 594 22 23 594.000 22.000 23 40 15.620 7.459
Tilted Ellipsoids.tif 15 2661 638 22 23 638.000 22.000 23 40 15.620 7.459
Tilted Ellipsoids.tif 16 2679 682 22 23 682.000 22.000 23 40 15.124 7.466
Tilted Ellipsoids.tif 17 2689 726 22 23 726.000 22.000 23 40 15.218 7.740
Tilted Ellipsoids.tif 18 2671 770 22 23 770.000 22.000 23 40 14.495 7.866
Tilted Ellipsoids.tif 19 2621 814 22 23 814.000 22.000 23 40 16.000 8.000
Tilted Ellipsoids.tif 20 2671 858 22 23 858.000 22.000 23 40 14.495 7.866
Tilted Ellipsoids.tif 21 2689 902 22 23 902.000 22.000 23 40 15.218 7.740
Tilted Ellipsoids.tif 22 2683 946 22 23 946.000 22.000 23 40 15.126 7.461
Tilted Ellipsoids.tif 23 2661 990 22 23 990.000 22.000 23 40 15.620 7.459
Tilted Ellipsoids.tif 24 2661 1034 22 23 1034.000 22.000 23 40 15.620 7.459
Tilted Ellipsoids.tif 25 2679 1078 22 23 1078.000 22.000 23 40 15.124 7.466
Tilted Ellipsoids.tif 26 2689 1122 22 23 1122.000 22.000 23 40 15.218 7.740
Tilted Ellipsoids.tif 27 2671 1166 22 23 1166.000 22.000 23 40 14.495 7.866
Tilted Ellipsoids.tif 28 2621 1210 22 23 1210.000 22.000 23 40 16.000 8.000
Tilted Ellipsoids.tif 29 2671 1254 22 23 1254.000 22.000 23 40 14.495 7.866
Tilted Ellipsoids.tif 30 2689 1298 22 23 1298.000 22.000 23 40 15.218 7.740
Tilted Ellipsoids.tif 31 2683 1342 22 23 1342.000 22.000 23 40 15.126 7.461
Tilted Ellipsoids.tif 32 2689 1518 22 771 1518.000 22.000 771.000 38.965 15.690 7.693
Tilted Ellipsoids.tif 33 2701 770 22 67 770.000 22.000 67 38.854 15.855 7.882
Tilted Ellipsoids.tif 34 2671 1078 22 507 1078.000 22.000 507 39.289 15.628 7.927

These blobs are osteocyte lacunae:

Image stack is here:
small_osteocyte_lacunae_binary.tif.zip (34.8 KB)

Numerical results look like this:

Label ID Vol. (mm³) x Cent (mm) y Cent (mm) z Cent (mm) Box x (mm) Box y (mm) Box z (mm) Box l0 (mm) Box l1 (mm) Box l2 (mm)
small_osteocyte_lacunae_binary.tif 1 1.263E-7 0.010 0.019 0.003 0.010 0.019 0.003 0.010 0.008 0.003
small_osteocyte_lacunae_binary.tif 2 3.753E-8 0.066 0.019 7.311E-4 0.067 0.019 8.560E-4 0.007 0.003 0.003
small_osteocyte_lacunae_binary.tif 3 2.664E-8 0.051 0.024 9.293E-4 0.051 0.024 9.659E-4 0.005 0.003 0.003
small_osteocyte_lacunae_binary.tif 4 3.628E-7 0.046 0.037 0.005 0.046 0.037 0.004 0.014 0.011 0.005
small_osteocyte_lacunae_binary.tif 5 6.256E-8 0.032 0.048 0.001 0.033 0.048 0.002 0.008 0.005 0.003
small_osteocyte_lacunae_binary.tif 6 2.087E-7 0.006 0.051 0.005 0.006 0.051 0.006 0.013 0.008 0.003
small_osteocyte_lacunae_binary.tif 7 2.623E-8 0.053 0.050 6.935E-4 0.053 0.050 5.426E-4 0.005 0.004 0.002
small_osteocyte_lacunae_binary.tif 8 2.058E-8 0.064 0.064 2.318E-4 0.064 0.064 3.375E-4 0.007 0.003 9.439E-4
small_osteocyte_lacunae_binary.tif 9 1.278E-6 0.025 0.060 0.017 0.025 0.061 0.016 0.036 0.011 0.007
small_osteocyte_lacunae_binary.tif 10 1.715E-7 0.083 0.076 0.007 0.084 0.075 0.008 0.018 0.005 0.004
small_osteocyte_lacunae_binary.tif 11 4.278E-8 0.014 0.080 0.001 0.014 0.079 0.001 0.006 0.004 0.003
small_osteocyte_lacunae_binary.tif 12 7.103E-8 0.073 0.087 0.004 0.073 0.087 0.004 0.010 0.005 0.002
small_osteocyte_lacunae_binary.tif 13 2.357E-7 0.050 0.081 0.008 0.050 0.080 0.009 0.017 0.008 0.003
small_osteocyte_lacunae_binary.tif 14 1.848E-7 0.084 0.060 0.011 0.084 0.060 0.011 0.019 0.005 0.003
small_osteocyte_lacunae_binary.tif 15 1.909E-7 0.070 0.048 0.011 0.070 0.048 0.010 0.016 0.005 0.004
small_osteocyte_lacunae_binary.tif 16 2.991E-7 0.025 0.083 0.009 0.024 0.084 0.009 0.015 0.010 0.006
small_osteocyte_lacunae_binary.tif 17 3.265E-7 0.032 0.005 0.012 0.032 0.005 0.013 0.019 0.010 0.003
small_osteocyte_lacunae_binary.tif 18 4.803E-8 9.315E-4 0.033 0.008 0.001 0.033 0.008 0.007 0.004 0.002
small_osteocyte_lacunae_binary.tif 19 2.991E-7 0.063 0.003 0.014 0.064 0.002 0.014 0.016 0.011 0.003

It looks great, Michael, very close to what I had in mind. Thank you very much for the effort.
I have been trying to implement it, but I have some problems with cloning it from github. I will keep trying and let you know how that goes.
Best,

A minimal axis-aligned bounding box feature is now in general release, in BoneJ (styloid-r11) and will be added to user installations after running the ImageJ updater. User docs are here.