Measuring Skeletal Length

I created a macro that creates a skeleton of an object and then I want to measure the curved length because feret and bounding box lengths are inappropriate.

Using the segmented line tool my object has an approximate length of 1045 pixels.

Previously I read on the email list that an easy way to get skeletal length is to take the perimeter of a selection of the skeleton and divide by two. I followed this advice, but a colleague pointed out this morning that the discrepancy between the length measured in this way and the length measured using a segmented line were sometimes as much as 30% different.

Can anyone help me understand the discrepancy or offer up a better solution to measuring skeletal length.

Using the attached example image and the following macro code:

//open example.tif
open();

// segmented line ~1045 px

setAutoThreshold("Default");
setOption("BlackBackground", false);
run("Convert to Mask");

// fix the f/ing LUT
invertedLUT = is("Inverting LUT");
if (invertedLUT == 1) {
  run("Invert LUT");
  run("Invert");
}

run("Create Selection");
run("Measure");

// perim/2 == 2296/2 = 1148

// how many black pixels are there in the image? 941
getHistogram(values, counts, 256);  
LENG=counts[0]; 
print(LENG);

I’m using ImageJ 2.0.0-rc-43/1.50h35; Java 1.8.0_66 64 bit on Mac OSX 10.11.4

Hello @Brandon_Hurr,

What about you use the Analyze Skeleton plugin? It should be right what you’re looking for.

ignacio

3 Likes

@iarganda,

I believe that maximum branch length is what I’m looking for. Can you help me understand what the source code is doing to measure length?

It’s clearly longer than the number of pixels observed so there is some interpolation going on but I don’t understand the code enough to know what’s going on.

Thanks,
Brandon

1 Like

Hi @Brandon_Hurr ,

I don’t know how exactly the Analyze skeleton algorithm calculates the length but in general you would not need to go via the perimeter, since the skeleton is a one-pixel line. This means its area = length. This assumption is incorrect, see comment from @anon96376101 [quote=“Herbie, post:5, topic:1262”]
Think of diagonal pixels, i.e. sqrt(2)
[/quote]

.
The discrepancy you (@Brandon_Hurr) descibe is coming from the way how the skeleton is created which always results in a shortening at both ends.

You might wanna give the “Particle Length (via Skeleton)” from the BioVoxxel toolbox update site. This algorithm takes a binary image, creates the skeleton and compensates the missing parts at the ends to overcome those discrepancy.

But if you have branches this simple algorithm measures the complete length of all branches in one tree.

Here is the extracted IJ macro code:

    setOption("BlackBackground", true);
    	original = getTitle();
    	if (is("binary")==false) {
    		exit("8-bit binary image necessary"); 
    	}
    	//setup dialog
    	Dialog.create("Setup");
    			Dialog.addString("Size:", "0-Infinity");
    			Dialog.addString("Circularity:", "0.00-1.00");
    			Dialog.addChoice("Show", newArray("Nothing", "Outlines", "Bare Outlines", "Ellipses", "Masks", "Count Masks", "Overlay Outlines", "Overlay Masks"), "Nothing");
    			Dialog.addCheckbox("clear results", false);
    			Dialog.addCheckbox("include holes", false);
    			Dialog.addCheckbox("exclude edges", false);
    			Dialog.addCheckbox("mininmal results", true);
    			Dialog.show();
    			size = Dialog.getString();
    			circ = Dialog.getString();
    			show = Dialog.getChoice();
    			clear = Dialog.getCheckbox();
    			holes = Dialog.getCheckbox();
    			edges = Dialog.getCheckbox();
    			display = Dialog.getCheckbox();
    			if(clear==true) {
    				clear = "clear ";
    			} else {
    				clear = "";
    			}
    			if(holes==true) {
    				hole = "include ";
    			} else {
    				hole = "";
    			}
    			if(edges==true) {
    				edge = "exclude ";
    			} else {
    				edge = "";
    			}
    	
    	setBatchMode(true);
    	
    	//create distance map of particles
    	selectWindow(original);
    	run("Duplicate...", "title=distance");
    	distance = getTitle();
    	run("Distance Map");
    	
    	//create particle skeleton
    	selectWindow(original);
    	run("Duplicate...", "title=skeleton");
    	skeleton = getTitle();
    	run("Skeletonize");
    	
    	//create endpoint erosded skeleton
    	run("Duplicate...", "title=eroded");
    	eroded = getTitle();
    	run("Options...", "iterations=1 count=7 black edm=Overwrite do=Erode");
    	
    	//create binary endpoint image
    	run("Calculator Plus", "i1="+skeleton+" i2="+eroded+" operation=[Subtract: i2 = (i1-i2) x k1 + k2] k1=1 k2=0 create");
    	endpoints = getTitle();
    	
    	//create image with endpoint intensities 
    	run("Select All");
    	run("Copy");
    	selectWindow(distance);
    	run("Select All");
    	setPasteMode("Transparent-white");
    	run("Paste");
    	run("Select None");
    	
    	correction = getTitle();
    	
    	
    	//results display mode and redirect to the endpoint intensity image 
    	if(display==true) {
    		run("Set Measurements...", "area centroid integrated redirect=["+correction+"] decimal=3");
    	} else {
    		run("Set Measurements...", "area mean standard modal min centroid center perimeter bounding fit shape feret's integrated median skewness kurtosis area_fraction redirect=["+correction+"] decimal=3");
    	}
    	
    	//analysis
    	selectWindow(skeleton);
    	run("Analyze Particles...", "size="+size+" circularity="+circ+" show=["+show+"] display " + edge + clear + hole +" record");
    	setBatchMode("show");
    	//correct the length by addition of skeleton area to endpoint intensities
    	for(n=0; n<nResults;n++) {
    		setResult("Corr. Length", n, ((getResult("Area", n))+(getResult("RawIntDen", n))));
    	}
    	//close intermediate images
    	close(distance);
    	close(skeleton);
    	close(eroded);
    	close(endpoints);
    	close(correction);
    	setBatchMode(false);

Edit: remark: the macro presented above does nott take the longer diagonal length of the digital binary pixel into account!

1 Like

This is must be doubted. Think of diagonal pixels, i.e. sqrt(2). See the more recent (26. Feb.) discussion on the IJ-list with topic “automate distance measure of uneven surface”

Best

Herbie

1 Like

@biovoxxel

When I create the skeleton, I compensate for the missing ends of my object by extending a line along the slope of the last 40 pixels of the skeleton to the tip of the object. It works fairly well, but is not perfect. Once I do that I have a full skeleton like in my example.

That being the case, your suggestion of area would be the same as counting the pixels from the histogram. I do some smoothing and other things to ensure I only have a single branch (mostly works).

@anon96376101,
True! My mistake. However, to correct this in the macro would most likely be too slow calculation wise.
Thus, I keep the macro code here if someone is interested in but with the remark of the non-corrected diagonal pixel size.

1 Like

@anon96376101 @biovoxxel @iarganda

Ok, so I altered the code found in that thread and followed it.

/// different method using "automate distance measure of uneven surface" Feb 26 ImageJ list

//open example.png
open();

// segmented line ~1045 px

setAutoThreshold("Default");
setOption("BlackBackground", false);
run("Convert to Mask");

// fix the f/ing LUT
invertedLUT = is("Inverting LUT");
if (invertedLUT == 1) {
  run("Invert LUT");
  run("Invert");
}

doWand(559, 81, 0, "Legacy");  //assumes (0,0) is air

getSelectionCoordinates(x,y);
Array.print(x);
Array.print(y);

makeSelection("Polyline", x, y);
run("Set Measurements...", "perimeter redirect=None decimal=3");
run("Measure");

Now the code produces a perimeter of 2665. I assume you divide that by 2, which would be 1332.5.

So I have:
941 px (area and histogram counts)
973 px Maximum branch length using Analyze Skeleton
1045 px if using segmented line tool (by hand, not automatable)
1148 px if creating a selection perim/2
1332 px if using wand tool perim/2

I’m not sure the “most correct” solution is any more clear to me. :confused:

Good day Brandon,

sorry, but I’ve nothing to do with this code or any code others have written for skeleton analyses. I’ve only remarked that the assumption “This means its area = length.” must be doubted.

Best

Herbie

@iarganda @biovoxxel
Ok, to follow up on @anon96376101 's point about sqrt(2), I calculated the distance between every point along the skeleton and added those up. I also realized that my image isn’t actually a skeleton, there were a few points that caused issues. Added skeletonization step.

I get 1099 px length, which is very close to the 1045 that was done by hand and selection perim/2.

What I don’t get is how this differs from @iarganda’s analyze skeleton method. As far as I can tell, it’s using pythag’s theorem to calculate distance between points in the same way.

/// different method using pythag to calculate distance between points

//open example.png
open();



// make binary image
setAutoThreshold("Default");
setOption("BlackBackground", false);
run("Convert to Mask");

// fix the f/ing LUT
invertedLUT = is("Inverting LUT");
if (invertedLUT == 1) {
  run("Invert LUT");
  run("Invert");
}

run("Invert"); // skeletons are white on black invert to reflect this

run("Skeletonize (2D/3D)");

// find end points
rename("Image1"); 
run("Duplicate...", "title=[Image2]");
run("Duplicate...", "title=finalline");
selectWindow("Image1"); 
run("Options...", "iterations=1 count=7 black edm=Overwrite do=Erode"); 
imageCalculator("Subtract create", "Image2","Image1"); 
selectWindow("Result of Image2"); 


// Get "Head" and "Tail" points 
run("Find Maxima...", "noise=10 output=[Point Selection]"); 
getSelectionCoordinates(x, y); 	

// head (blossom end due to image orientation)
// tail (stem end due to image orientation)

HEADX=x[1]; 
HEADY=y[1]; 

TAILX=x[0]; // needed only for angle measurement
TAILY=y[0]; // needed only for angle measurement

// close unneccesary windows
selectWindow("Image1");
run("Close");
selectWindow("Image2");
run("Close");
selectWindow("Result of Image2");
run("Close");


// This block of code goes along the skeleton and finds every pixel position
// This allows you to find the width of the object at any position along the skeleton 
selectWindow("finalline"); 
getHistogram(values, counts, 256);  
LENG=counts[255]; 

COORDX=newArray(LENG-1); 
COORDY=newArray(LENG-1); 

Array.print(COORDX);
Array.print(COORDY);
        
// To scan the neighboors 
selectWindow("finalline");
run("Duplicate...", "title=finalline2");
selectWindow("finalline");
seqalongline(HEADX, HEADY, LENG, COORDX, COORDY); // function crawls along line reports back coords x/y

Array.print(COORDX);
Array.print(COORDY);

selectWindow("finalline");


distArray = newArray(LENG-2); 
distValue = 0;

for (i=0; i<LENG-2; i++)  { 
  distArray[i] =  pythag(COORDX[i], COORDY[i], COORDX[i+1], COORDY[i+1]);
  distValue = distValue + distArray[i];
}

Array.print(distArray);
print(distValue);


/////////////////////////////////////////////////////////////////////////
function pythag (x1,y1, x2,y2) {
  return sqrt( ((x2-x1)*(x2-x1)+ (y2-y1) * (y2-y1)));
}
/////////////////////////////////////////////////////////////////////////


/////////////////////////////////////////////////////////////////////////
  // This block of code goes along the skeleton and finds every pixel position
  // This allows you to find the width of the object at any position along the skeleton 
  // required input:
  //	- 8-bit *B&W* image with skeleton to follow in foreground
  //	- start position (startX, startY)
  //	- COORDX & COORDY array of length LENG

function seqalongline(startX, startY, LENG, COORDX, COORDY) {
  //reset HEADX and HEADY
  HEADX = startX;
  HEADY = startY;

  X_OFFSET = newArray(  0,  1,  1,  1,  0, -1, -1, -1 ); 
  Y_OFFSET = newArray( -1, -1,  0,  1,  1,  1,  0, -1 ); 
  
  for (j=0; j<LENG-1; j++)  { 
    print(HEADX);
    print(HEADY);
    
    for (i=0; i<8; i++) { 
    
      neighbor = getPixel(HEADX+X_OFFSET[i], HEADY+Y_OFFSET[i]); 
        
      if (neighbor == 255) { 
        COORDX[j]=HEADX; 
        COORDY[j]=HEADY; 
        setPixel(HEADX, HEADY, 100); 
        HEADX=HEADX+X_OFFSET[i]; 
        HEADY=HEADY+Y_OFFSET[i]; 
        i=8; // break (quit the loop) 

      } 
    } 
  }
}
/////////////////////////////////////////////////////////////////////////
1 Like

Hello @Brandon_Hurr,

Sorry to jump late into the conversation.

That’s very important. Analyze Skeleton will only work if the input image is a skeleton following the description of the Skeletonize3D plugin page.

I’m having trouble reproducing your result using the macro you posted. My guess is the coordinates are not exactly in the center positions of the skeleton voxels.

Good point @iarganda . I reran the analyze skeleton and it came up with 1100. So it’s:

  • 941 px (area and histogram counts)
  • 1045 px if using segmented line tool (by hand, not automatable)
  • 1099 px Macro pythag distance
  • 1100 px Maximum branch length using Analyze Skeleton
  • 1148 px if creating a selection perim/2
  • 1332 px if using wand tool perim/2

There is good agreement between my macro version of your plugin and your plugin. No sense in reinventing the wheel so I’ll just use your analyzer’s maximum branch.

Thanks for helping me understand. The perimeter method is too unstable for accurate results.

3 Likes

As Herbie pointed out, area is a bad estimator of line length.
I wrote some time ago a line version of Particles8 (called Lines8, included in the Morphology zip file at my page) which measures the length of a line in 2 different ways. One is based on the length of the perimeter making sure that you do not count the same pixel twice (what you call perimeter/2). This works fine if there are no ‘closed loops’ in the line/skeleton, like in the image above. The other method is based on connectivity of a skeleton so you measure the distance between connected pixels (1, sqrt(2)). Depending on what you are doing, you might want one or the other.
On the example image above, I get the same results for both methods:
Perimeter/2: 1100.785889
Skeleton: 1100.785889
Pixels: 939
Area: 0

I guess that the difference to your figures is because of they way the perimeter lengths are computed. Lines8 uses the Freeman algorithm to compute these. IJ, I think, uses something different.
Also what you call “Area” is in reality ‘number of pixels’. If you consider pixels as sampled points, and the perimeter is the path across the boundary pixels, then it should follow that the ‘area’ is strictly the ‘area confined by polygon made of the boundary pixels’, not the number of pixels or samples. They are close, but they are not the same. Lines and points should have no area, although they can be made of many pixels.
Particles8 reports both, area and pixels. The particle analyzer in IJ reports pixels as area.

Strange (i.e. unintuitive) things happen in discrete lattices, e.g. two digital lines can cross without overlapping/sharing any pixels!

8 Likes

Hello everyone,
This in general is a statistics problem, and current research uses the factor of dividing by Pi.
Bob