CLIJ2 Tribolium embryo morphometry sample macro

Hi @haesleinhuepf,
I’m trying to go through the Tribolium embryo morphometry sample macro but the commands that start with show are not recognised.
image

image

CLIJ and CLIJ2 are installed and up to date.
image

Am I missing something?

Thank you!

1 Like

Hey @LPUoO,

the ‘show’ function is at the end of the macro you were looking at.

I put it there to not confuse the reader at the beginning.

Let us know if this solves your issue!

Cheers,
Robert

1 Like

Ahh, yes thank you,
I was attempting to run small individual fragments of the macro so I could see exactly the effect of the individual commands!

Problem solved!

1 Like

@haesleinhuepf,
I am having some trouble running the end of the " Quantitative analysis of distance between neighbors" section.

For some reason all the outputs are blank (pixel values = NaN). I’m unsure why? The only thing I did was to disable line 150 (exit()).

Also, reading the description I am slightly confused by the diference between the distance map and the neighbor mean distance map.
The distance map is the the averge distance between a node and of all its neighbors while the neighbor mean distance map would be the mean of the averge distance between a node and of all its neighbors?

Thank you

1 Like

Hey @LPUoO,

this bug was discussed earlier here in the forum. The issue was that the parametric 3D stacks contained NaN pixels and a maximum-z-projection of some pixels were any are NaN, is NaN. It was fixed by introducing some lines in the example code:

Apparently, the website wasn’t update after that. Apologies. I’m going to update that asap.

Yes, you got it right. In the first place, you measure distance between neighbors generating a distance matrix. You then average the distance for each cell and visualise it in a map. Afterwards, you compute the mean distance of all neighbors of the neighbors of each cell. If you think about that further, you are averaging distance of first and second order while taking first order into account twice. Thus, this something like a Gaussian Blur or local distances.

Let me know if you want to learn more! This neighbor-part is the interesting part, also for me :wink:

Cheers,
Robert

1 Like

Addendum: In daily routine it might be easier to replace NaNs with zeros before doing the maximum projection.

1 Like

Hey @LPUoO,

I just updated the website. While doing so, I realised that the bug you observed was still there, when running the script on an NVidia GPU, but not on Intel. It’s fixed now on the website and in the example code as well. I need to dig deeper regaring the difference between results on NVidia and Intel. Thanks for pointing me to this issue!

Cheers,
Robert

1 Like

Yes, I am very interested in learning more about this.

Also. a few additional question. I am now trying to apply your macro to my data.

  1. How do I incorporate my voxel dimension in the distance calculations?
    image

  2. In the distance map is it possible to exclude maximas where the average distance between a node and of all its neighbors are superior to a given value X?

  3. How can I export the average distance between a node and of all its neighbors as a simple list?

Thank you very much!

Hey @LPUoO

that’s great questions!

Multiply the point-list with the voxel size as shown here

Yes! After measuring the distance, you have a vector of measurements in a 1d-image. Apply a threshold to that vector, e.g. using the smallerConstant method and you receive a binary vector, which specifies which rows/columns to keep in the distance matrix. Depending on if you want to eliminate rows or columns, you might want to transpose the vector before multiplying the vector with the distance matrix. Use multiplyImages for this operation and not matrix multiplication :wink:

If you have measurements of any kind you can print them out or pull them to the results table

You’re most welcome :slight_smile:

1 Like

Thank you @haesleinhuepf

I’m having some trouble incorporating my voxel dimension in the distance calculations.
I multiply the point_list with the voxel_size matrix and I then use the transformed_pointlist calculate the touch matrix mesh,


but my distance matrix, distance mesh and distance map are blank, so I am not sure where is my mistake? :thinking:



I did attempt to use the smallerConstant function but sadly also failed.
Would thresholding the distance matrix be an alternative to your suggestion?
image

Without thresholding the distance matrix we have this


but this after thresholding, so I assume this is working okay?


I did mange to pull the averageDistanceOfTouchingNeighbors to the results table so :+1:

Thank you very much!

In case touchMatrixToMesh doesn’t work, try printing out the parameters you send to it. I think you should keep the pointlist when calling touchMatrixToMesh because the mesh should be drawn in pixel coordinates, right? Assume your image is 100 pixels wide and your pixel size 0.01: The mesh would then be drawn inside one pixel in the top left corner. :wink:

Can you crop out the top left 10x10 pixels of the distance matrix and print it out? You will then see what numbers are in there and can think about a proper thresholding strategy. Furthermore, if code using smallerConstant doesn’t work, please post it here so that we can take a look :wink:

Keep pushing pixels, I think you’re almost there! :slight_smile:

Cheers,
Robert

1 Like

Thank you @haesleinhuepf,
I made some progress

My issue was that when I did:

Ext.CLIJ2_pushArray(voxel_size, newArray(0.2, 0.2, 0.5), 1, 3, 1);
Ext.CLIJ2_multiplyImages(pointlist, voxel_size, pointlist_calibrated);

Then pointlist_calibrated was just blank. I’m not sure why.

As of now, my ridiculous solution is :face_with_hand_over_mouth:

Ext.CLIJ2_getDimensions(pointlist, Number_width, Number_height, Number_depth);
newImage("Image voxel", "32-bit black", Number_width, Number_height, Number_depth);

for(x=0; x<Number_width; x++){
	setPixel(x,0,2.7676);
	setPixel(x,1,2.7676);
	setPixel(x,2,5);
}
Ext.CLIJ2_push("Image voxel");
Ext.CLIJ2_multiplyImages(pointlist, "Image voxel", pointlist_calibrated);

But I’m pretty sure this is not the solution you had in mind.
My only issue is that the distance matrix “doesn’t fit” any more in the image :sleepy:. The distance map looks fine. :+1:


Is there any major advantage in incorporating the voxel dimension in the distance calculations? Because it looks much simpler to make the voxels isotropic by re-slicing the stack then if necessary multiply the results by voxel size.


I think I managed to do this with help from the superpixel_segmentation tutorial

run("Clear Results");
Ext.CLIJ2_statisticsOfBackgroundAndLabelledPixels(distance_map, labels);
Ext.CLIJ2_pushResultsTableColumn(MEAN_INTENSITY, "MEAN_INTENSITY");
//"MEAN_INTENSITY" is actually the distance extracted from the distance map.
run("Clear Results");

Ext.CLIJ2_replaceIntensities(labels, MEAN_INTENSITY, MEAN_INTENSITY_map);
MEAN_INTENSITY_threshold = 20;

threshold_vector_and_visualise(MEAN_INTENSITY, labels, MEAN_INTENSITY_threshold);

// This function takes a vector, binarizes it by using a threshold and 
// visualizes the results as parametric image by the given label map:
function threshold_vector_and_visualise(vector, labelmap, threshold) {

	// threshold the vector in two vectors:
	Ext.CLIJ2_smallerConstant(vector, small_objects, threshold);
	Ext.CLIJ2_greaterOrEqualConstant(vector, large_objects, threshold);
	
	// visualise resulting binary images
	Ext.CLIJ2_replaceIntensities(labelmap, small_objects, small_objects_map);
	show(small_objects_map, "below threshold");	
}
Ext.CLIJ2_push("below threshold");
Ext.CLIJ2_multiplyImages("below threshold", distance_map, filtered_distance_map);
Ext.CLIJ2_pull(filtered_distance_map);
run("glasbey_on_dark");

Thank you very much for all your help

Hm. If I execute this:

run("CLIJ2 Macro Extensions", "cl_device=");

Ext.CLIJ2_pushArray(pointlist, newArray(
	1, 2, 3, 4,
	1, 2, 3, 4,
	1, 2, 3, 5
), 4, 3, 1);

Ext.CLIJ2_pushArray(voxel_size, newArray(0.2, 0.2, 0.5), 1, 3, 1);
Ext.CLIJ2_multiplyImages(pointlist, voxel_size, pointlist_calibrated);

Ext.CLIJ2_print(pointlist_calibrated);

It outputs that:
image

Which appears correct to me. Now the big question is: What’s different in your code? Feel free to provide a full executable example macro

Yeah, your solution is indeed ridiculous. Feel free to try also this ridiculous solution :slight_smile:

run("CLIJ2 Macro Extensions", "cl_device=");

width = 10;
height = 3; 
bit_depth = 32;
image = "Image voxel";
Ext.CLIJ2_create2D(image, width, height, bit_depth);

Ext.CLIJ2_setRow(image, 0, 2.7676);
Ext.CLIJ2_setRow(image, 1, 2.7676);
Ext.CLIJ2_setRow(image, 2, 5);

Ext.CLIJ2_print(image);

image

For drawing the mesh, you should take the distance matrix as you got it from generateDistanceMatrix and for this, you should take the pointlist without multiplying voxel sizes. Meshes are drawn in pixel coordinate system. CLIJ doesn’t have any operation taking physical units into account.

Hm, not sure. Resampling brings issues. But I also prefer making every dataset isotropic first and then have an easier life in the rest of the workflow. In my data (looking at nuclei and their positions), often voxel size 1x1x1 um works nicely and makes interpretation of distance measurements super easy.

Let me know if this helps!

Cheers,
Robert

Thank you @haesleinhuepf

The distance mesh is now sorted out :+1:

Below is my macro.
Make sure that the path in line 7 is correct.
The issue seems to be between lines 88 and 108.
A dialog will let you run your method (that works) or mine (that doesn’t).

run("CLIJ2 Macro Extensions", "cl_device=[GeForce RTX 2060]");
Ext.CLIJ2_clear();

run("Close All");

// Load a data set
path = "C:/structure/teaching/clij2_example_data/";
open(path + "lund1051_resampled.tif");
input = getTitle();


Ext.CLIJ2_push(input);
//run("Close All");

// visualize the dataset
//show(input, "input");
/*
## Spot detection
After some noise removal/smoothing, we perform a local maximum detection:
*/
// gaussian blur
sigma = 2;
Ext.CLIJ2_gaussianBlur3D(input, blurred, sigma, sigma, sigma);

// detect maxima
radius = 2.0;
Ext.CLIJ2_detectMaximaBox(blurred, detected_maxima, radius);
show_spots(detected_maxima, "detected maxima");
/*
## Spot curation
Now, we remove spots with values below a certain intensity and label the remaining spots.
*/
// threshold
threshold = 300.0;
Ext.CLIJ2_threshold(blurred, thresholded, threshold);

// mask
Ext.CLIJ2_mask(detected_maxima, thresholded, masked_spots);

// label spots
Ext.CLIJ2_labelSpots(masked_spots, labelled_spots);
show_spots(labelled_spots, "selected, labelled spots");
run("glasbey_on_dark");
/*
Let's see how many spots are left:
*/
Ext.CLIJ2_getMaximumOfAllPixels(labelled_spots, number_of_spots);
print("Number of detected spots: " + number_of_spots);
/*
## Expanding labelled spots
Next, we spatially extend the labelled spots by applying a maximum filter.
*/
// label map closing
number_of_dilations = 10;
number_of_erosions = 4;
Ext.CLIJ2_copy(labelled_spots, flip);
for (i = 0; i < number_of_dilations; i++) {
	Ext.CLIJ2_onlyzeroOverwriteMaximumBox(flip, flop);
	Ext.CLIJ2_onlyzeroOverwriteMaximumDiamond(flop, flip);
	if (i % 2 == 0) {
		//show(flip, "Extended spots after " + (i * 2) + " dilations");
		run("glasbey_on_dark");
	}
}
/*
Afterwards, we erode all labels in the map and get a final result of cell segementation.
*/
Ext.CLIJ2_threshold(flip, flap, 1);
for (i = 0; i < number_of_erosions; i++) {
	Ext.CLIJ2_erodeBox(flap, flop);
	Ext.CLIJ2_erodeBox(flop, flap);
}
Ext.CLIJ2_mask(flip, flap, labels);
show(labels, "cell segmentation");
run("glasbey_on_dark");
/*
We also save all labels to disc to use them as starting point in other notebooks, later.
*/
Ext.CLIJ2_pull(labels);
saveAs("TIF", path + "lund1051_labelled.tif");
close();
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
/*
## Draw connectivity of the cells as a mesh
We then read out all current positions of detected nuclei as a pointlist to generate 
a distance matrix of all nuclei towards each other:
*/
Ext.CLIJ2_labelledSpotsToPointList(labelled_spots, pointlist);
Ext.CLIJ2_pull(pointlist);
rename("pointlist");
Ext.CLIJ2_print(pointlist);
#@ String(choices={"haesleinhuepf's method (works)", "LPUoO's method (doesn't work)"}, style="radioButtonHorizontal") method

Ext.CLIJ2_getDimensions(pointlist, Number_width, Number_height, Number_depth);
voxel_size = "Image voxel";
	if (method == "haesleinhuepf's method (works)") {
		Ext.CLIJ2_create2D(voxel_size, Number_width, Number_height, Number_depth);
		Ext.CLIJ2_setRow(voxel_size, 0, 2.7676);
		Ext.CLIJ2_setRow(voxel_size, 1, 2.7676);
		Ext.CLIJ2_setRow(voxel_size, 2, 5);
		}
	else {
		Ext.CLIJ2_pushArray(voxel_size, newArray(2.7676, 2.7676, 5), 1, 3, 1);
	}

Ext.CLIJ2_multiplyImages(pointlist, voxel_size, pointlist_calibrated);
Ext.CLIJ2_pull(pointlist_calibrated);
rename("pointlist_calibrated");


Ext.CLIJ2_generateDistanceMatrix(pointlist_calibrated, pointlist_calibrated, distance_matrix);
show(distance_matrix, "distance matrix");


/*
Starting from the label map of segmented cells, we generate a touch matrix:
*/
Ext.CLIJ2_generateTouchMatrix(labels, touch_matrix);
show_spots(touch_matrix, "touch matrix");
//Ext.CLIJ2_print(touch_matrix);


/*
Using element by element multiplication of a distance matrix and a touch matrix, we calculate the length of 
each edge. We use this result to draw a mesh with a color gradient of distance (between 0 and 50 micron):
*/
Ext.CLIJ2_multiplyImages(touch_matrix, distance_matrix, touch_matrix_with_distances);
Ext.CLIJ2_getDimensions(input, width, height, depth);
Ext.CLIJ2_create3D(mesh, width, height, depth, 32);
Ext.CLIJ2_touchMatrixToMesh(pointlist, touch_matrix_with_distances, mesh);//////////////////////////////////////////////////////////////////////////////////////
show(mesh, "distance mesh");
run("Green Fire Blue");
setMinAndMax(0, 50);
/*
## Quantitative analysis of distance between neighbors
Next, we determine the averge distance between a node and of all its neighbors. The resulting 
vector has as many entries as nodes in the graph. We use this vector to color-code the 
label map of segmented cells. This means, label 1 gets replaced by the average distance to 
node 1, label 2 by the average distance to node 2, et cetera.
*/
Ext.CLIJ2_setColumn(touch_matrix, 0, 0);


Ext.CLIJ2_averageDistanceOfTouchingNeighbors(distance_matrix, touch_matrix, distances_vector);
// set the first column to zero to ignore all object touching the background
//Ext.CLIJ2_setColumn(distances_vector, 0, 0);
//Ext.CLIJ2_setColumn(touch_matrix, 0, 0);



Ext.CLIJ2_undefinedToZero(distances_vector, distances_vector1);
Ext.CLIJ2_replaceIntensities(labels, distances_vector1, distance_map);
show(distance_map, "distance map");
run("Fire");
setMinAndMax(0, 50);

//Ext.CLIJ2_pullToResultsTable(distances_vector);/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

/*
Also, let's see how much of GPU memory got used by this workflow. At the end, cleaning up remains important.
*/
//Ext.CLIJ2_reportMemory();

// finally, clean up
Ext.CLIJ2_clear();

/*
Following methods are convenient for a proper visualization in a notebook:
*/
function show(input, text) {
	Ext.CLIJ2_maximumZProjection(input, max_projection);
	Ext.CLIJ2_pull(max_projection);
	setColor(100000);
	drawString(text, 20, 20);
	Ext.CLIJ2_release(max_projection);
}

function show_spots(input, text) {
	Ext.CLIJ2_maximum3DBox(input, extended, 1, 1, 0);
	Ext.CLIJ2_maximumZProjection(extended, max_projection);
	Ext.CLIJ2_pull(max_projection);
	setColor(100000);
	drawString(text, 20, 20);
	Ext.CLIJ2_release(extended);
	Ext.CLIJ2_release(max_projection);
}
run("Tile");
print("FINISH");

Apologies the macro is a little bit messy.

What issues do you have in mind?

Thank you very much for your help :relaxed:

1 Like

Just the obvious ones:
If you sample down, you may loose resolution/information. You may miss nuclei then for example.
If you sample up, your processing performance may become worse. It takes longer.

1 Like