CLIJ_label map to ROI fast version

Hi
@haesleinhuepf
I have a suggestion for the labelmap to ROI manager plugin in CLIJ2/X. The current version is quite slow on large images. I have written some code that works fast and thought I’d get input about this.
It essentially works by getting the centroid of each label, creating a selection at each label centroid using the wand tool and adding this selection to the ROI Manager. Let me know if there are any bugs.

roiManager("reset");
//assumes active image is a label image
label_image=getTitle();
run("Clear Results");

// Init GPU
run("CLIJ2 Macro Extensions", "cl_device=");
Ext.CLIJ2_push(label_image);
//statistics of labelled pixels
Ext.CLIJ2_statisticsOfLabelledPixels(label_image, label_image);
Ext.CLIJ2_release(label_image);
Ext.CLIJ2_clear();

//get centroid of each label
selectWindow("Results");
x=Table.getColumn("CENTROID_X");
y=Table.getColumn("CENTROID_Y");

//use wand tool to create selection at each label centroid and add the selection to ROI manager
//will not add it if there is no selection or if the background is somehow selected
selectWindow(label_image);
for(i=0;i<x.length;i++)
{
	//use wand tool; quicker than the threshold and selection method
	doWand(x[i], y[i]);	
	intensity=getValue(x[i], y[i]);
	//if there is a selection and if intensity >0 (not background), add ROI
	if(selectionType()>0 && intensity>0) roiManager("add");
}
run("Clear Results");

Cheers
Pradeep

3 Likes

Hey @pr4deepr,

wow, what a great idea! I’m sure also @bramvdbroek will be happy to see this. I’ll translate it to Java as soon as I can (bit busy these days) unless anybody else wants to adapt the current code? :smirk:

Cheers,
Robert

2 Likes

Wait a second @pr4deepr, what if the centroid is out of the label?

2 Likes

oooo, good point!!!
Assuming this would be the case for Centre of mass as well, right???
I was thinking we could use reduce labels to labelled spots, but the final labelled spot is also on the centroid.

How about:
If it returns False for the following

if(selectionType()>0 && intensity>0)

It searches within the bounding box coordinates till the condition is True? and then adds it to ROI?
Its better than searching the whole image. Only issue is if there is another object within the bounding box?

2 Likes

This seems to work:

	label_image=getTitle();
	run("CLIJ2 Macro Extensions", "cl_device=");
	Ext.CLIJ2_push(label_image);
	//reindex the labels so 
	Ext.CLIJ2_closeIndexGapsInLabelMap(label_image, reindex);
		//statistics of labelled pixels
	Ext.CLIJ2_statisticsOfLabelledPixels(reindex, reindex);
	Ext.CLIJ2_release(label_image);
	Ext.CLIJ2_pull(reindex);
	Ext.CLIJ2_clear();


	//get centroid of each label
	selectWindow("Results");
	x=Table.getColumn("CENTROID_X");
	y=Table.getColumn("CENTROID_Y");

	//getting the identifiers as the values correspond to the label values
	identifier=Table.getColumn("IDENTIFIER");
	
	x1=Table.getColumn("BOUNDING_BOX_X");
	y1=Table.getColumn("BOUNDING_BOX_Y");
	x2=Table.getColumn("BOUNDING_BOX_END_X");
	y2=Table.getColumn("BOUNDING_BOX_END_Y");

	//use wand tool to create selection at each label centroid and add the selection to ROI manager
	//will not add it if there is no selection or if the background is somehow selected
	selectWindow(reindex);
	for(i=0;i<x.length;i++)
	{
		//use wand tool; quicker than the threshold and selection method
		doWand(x[i], y[i]);	
		intensity=getValue(x[i], y[i]);
		//if there is a selection and if intensity >0 (not background), add ROI
		if(selectionType()>0 && intensity>0) roiManager("add");
		//if there is no intensity value at the centroid, its probably coz the object is not circular
		// and centroid is not in the object
		else{
			//get the length of the bounding box
			x_b=x2[i]-x1[i];
			//get y coordinate
			y_temp=y1[i];
			for(j=x1[i];j<x_b;j+=2)
			{
				//print(j,y_temp);
				//wand tool at the initial coordinate
				doWand(j, y_temp);	
				intensity=getValue(j, y_temp);
				//only if intensity is greater than zero and if it matches the identifier value from stats table
				//identifier value ensures that it only finds object of that particular value
				if(intensity>0 && intensity==identifier[i])
				{
					//print("FOUND, intensity matches: "+intensity);
					roiManager("add");
					j=x_b+1000;
				}
				y_temp+=2;
			}
		}
	}
2 Likes

I’ve used the identifier within the stats table to ensure that the correct object is identified…
Let me know if this makes sense…

1 Like

That’s a rather suboptimal solution for parallel computing (such as on GPUs). I assume that would be very slow.

Technially it might be fine, but any kind of for-loop over X and Y (and Z) make things slow. OpenCL and CLIJ are fast because we define operations in a way that they can be approached in parallel. Therefor you need to unwrap your algorithm. When parallelizing this operation, ask yourself: What condition must any pixel fullfill to make a condition true? And how can one identify a pixel with a condition in parallel. That can be very tricky (such as in this case).

I would be for making a new pullRoundishLabelsToROIManager (better name suggestions welcome), which mentions in the documentation that it’s fast, but has an assumption: Centroids should lie within the object. Furthermore, split labels (two regions with the same label) would also not work I guess.

2 Likes

Agree… It may not scale with large datasets. It seems to work fast for me, but not a good solution for parallel computing…

Thanks for the feedback.
Just spitballing here, but what if you convert the labels to label edges? and then get coordinates for the label edges? Wouldn’t that be the outlines? Can this be used to create selections?
Again, just a thought…

1 Like

pullStarConvexLabelsToROIManager maybe?

1 Like

That’s a great idea. I’m just afraid that this increases the coding challenge - for me at least.

Great thoughts! I went down the same road in my mind at some point when thinking about edge-to-edge distances of labels. Also going along edges is not super easy to parallelize. Where would you start to walk along the edge and how can you prevent other walkers to walk along the same edge?

Btw. if you’re note aware yet of the visualizeOutlinesOnOriginal method. In combination with the “Hi” lookup table (a sibling of the HiLo lookup table) , you can easily visualize oulines. It’s made for the assistant but maybe also helpful for visualization in different contexts. It’s super fast and also works in 2D and 3D. It just doesn’t use ROIs :wink:

Together with MorpholibJs measurement operations, which also work in 3D and take label images and no ROIs, one could actually get rid of ROIs completely. :smiley:

image image

3 Likes

With the walkers, I don’t know how it works, but can each walker have a different or unique intensity value?

Thanks for this Robert.
The visualise outlines is extremely helpful… makes a big difference when fine tuning parameters…

1 Like

Nice workaround @pr4deepr, even though some issues remain with non-round labeling.
Amazing to see that the magic wand is soo much quicker than the current method, where converting threshold to selection is the actual slow step. (@haesleinhuepf This is also carried out on the CPU, right?)

I tested it on an image of 3000x3000 pixels with 1104 ROIs. Processing is >100-fold faster.
However, the new workaround still misses 9 out of 1104 (semi-roundish) ROIs in this case.

Ah, found it: your line for(j=x1[i];j<x_b;j+=2) should be for(j=x1[i];j<x2[i];j+=2).
Is there a specific reason why you step with 2 pixels?
Also, doWand(j, y_temp); can be moved inside the if(intensity>0 && intensity==identifier[i]) {... statement., to save some extra time.

To further speed up things, after the centroid, but before the real brute force search method, one could first look for nonzero pixels in a spiral fashion, e.g. like these:
image image
I created the first spiral with:

newImage("spiral", "16-bit black", 200, 200, 1);
getDimensions(width, height, channels, slices, frames);

pitch = 1;
omega = 15 * (PI/180);	//Angle, converted to radians

x = 0; y = 0; i = 0;
while (x <= width/2 || y <= height/2) {
	x = i*pitch*omega*sin(omega*i);
	y = i*pitch*omega*cos(omega*i);
	setPixel(x + width/2, y + height/2, i);
	i++;
}
resetMinAndMax;
run("Fire");

And the second (Archimedean) spiral using:

newImage("spiral", "16-bit black", 200, 200, 1);
getDimensions(width, height, channels, slices, frames);
size = 100;
pitch = 4;
angle = 0; r = 0; i = 0;

while(r <= width/2 || r <= height/2) {
    r = sqrt(i)*pitch;
    angle += atan(1/r);
    x = (r)*cos(angle*pitch);
    y = (r)*sin(angle*pitch);
    setPixel(x + width/2, y + height/2, i);
    i++;
}
resetMinAndMax;
run("Fire");

Best regards,
Bram

3 Likes

Thanks for the feedback and picking up my mistake @bramvdbroek . Appreciate it as its a learning experience for me.

That was just for speed really. My idea was to search every second pixel than every pixel. Its a bit arbitrary; could be higher if needed.

Good point!

Just so I understand this, this search strategy is for within the bounding box, right? The idea is that this search is more efficient than a brute-force search every pixel in bounding box, thus reducing time taken to find a nonzero pixel within the box?
Any preference over which search method?
Thanks for including the code…

BTW, I found that @romainGuiet has included the doWand option for label to ROI in his plugin LaRoMe, specifically here.

2 Likes

Don’t know if my 2 cents are helpful here and if that is implementable in CLIJ2 and makes sense for GPU processing. So,…
Some time ago, I was using @pr4deepr’s method with the centroid as well until I came across the problem of a center outside. The solution in a standard ImageJ Macro is to use the X_Start and Y_Start points for the wand selection, since those are always part of the object, independent of its shape.
That works perfectly. I just don’t know if CLIJ2 is reading those start coordinate in the same way ImageJ’s “record start” does in the Analyze Particle… function. Having the X-Start and Y-Start is generally a good backup to have an inside point at hand.
I can further elaborate on that if needed.

2 Likes

Thanks @biovoxxel
Had no idea about this.
Found this documentation:

The Record Starts option allows plugins and macros to recreate particle outlines using the doWand(x,y) function. The CircularParticles macro demonstrates how to use this feature.

https://imagej.nih.gov/ij/docs/menus/analyze.html

In the example macro above, if you run this command (with display results ticked):

run(“Analyze Particles…”, “minimum=1 maximum=999999 display clear record”);

It will show results table with these options:
image

where XStart and YStart only appear if “Record starts” is ticked in “Analyze Particles”. These correspond

I had a look and could not find this option in CLIJ. @haesleinhuepf , is there a similar functionality? If not, is it easy to add this (when you have time)? Based on what @biovoxxel is saying, then it wouldn’t matter what shape the object is… you get the XStart and YStart and run doWand on each x,y value.

Pradeep

2 Likes

Something like this doesn’t exist yet in CLIJ. Does anybody know how XStart and YStart are mathematically defined? Side note: This will only work for connected labels, but not for labels in general as mentioned above.

CLIJ in very general hardly supports ImageJ ROIs. I prefer working with label images because these are platform independent (in python for example) and also work in 3D (and in napari for example :wink: ). So better support for ROIs is not on the roadmap.

However, I’m happy to support anybody who plans to work on it. Just posting the link to the code again:

And a link to the clij plugin template:

And a link to the developer documentation:
https://clij.github.io/clij2-docs/development

Let me know how I can help.

Cheers,
Robert

2 Likes

That is true and it will not work in a 3D case (at least not in the sense of the original idea

3 Likes

Btw, the ImageJ 3D Suite by @ThomasBoudier can turn label images into ROIs as well. Did anybody try how fast it performs?

1 Like

As far as I tested that, it was pretty fast. However, I think it depends strongly on the amount of labels and stack thickness.

Here a small scale testing:

run("Particles");
part = getTitle();
run("Invert LUT");

for (i = 0; i < 10; i++) {
	run("Duplicate...", " ");
}
run("Images to Stack", "name=Stack title=particles use");
stack = getTitle();

run("3D Manager");
selectWindow(stack);
start = getTime();
Ext.Manager3D_Segment(128, 255);
afterLabel = getTime();
Ext.Manager3D_AddImage();
afterRoiCreation = getTime();
labelingTime = afterLabel - start;
roiCreationTime = afterRoiCreation - afterLabel;

print("labeling time = " + (afterLabel - start) + " ms");
print("Roi from label time = " + (afterRoiCreation - afterLabel) + " ms");

For me it took 600ms to get 5097 ROIs from the labels created on a 10 slice stack, while in a 30 slice stack with the same number of objects it took already 2220ms. For 50 slices it was 3325ms.
Looks pretty much linear to me following the slice number. Might be similar for object numbers (?)

Still a pretty short time for a coffee break :wink:

3 Likes

Hi Pradeep @pr4deepr and others,

Indeed, I thought of the spiral probing as intermediate search possibility within the ROI bounds, in case the centroid pixel is zero. If no pixel belonging to the correct label is found after doing the spiral you could still probe every pixel in the bounding box.
I had fun making a macro for the spirals yesterday, but didn’t get to implement it in your code.
Anyway, thanks to Jan’s @biovoxxel suggestion using X_start and Y_start that may not be necessary any more.

The ImageJ 3D Suite option indeed seems to be a great and fast choice as well (although it sends the ROIs to the 3D Manager of course).

Generally, Robert @haesleinhuepf has me convinced that working with label maps is harder, better, faster and stronger than working with ROIs. Until recently I was missing some visualization functionality, but that has at least partially been solved, see for instance this tutorial macro to generate outlines and ROI numbers.

3 Likes