Hi @haesleinhuepf,
I’m trying to go through the Tribolium embryo morphometry sample macro but the commands that start with show
are not recognised.
CLIJ and CLIJ2 are installed and up to date.
Am I missing something?
Thank you!
Hi @haesleinhuepf,
I’m trying to go through the Tribolium embryo morphometry sample macro but the commands that start with show
are not recognised.
CLIJ and CLIJ2 are installed and up to date.
Am I missing something?
Thank you!
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
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!
@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
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
Cheers,
Robert
Addendum: In daily routine it might be easier to replace NaNs with zeros before doing the maximum projection.
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
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.
How do I incorporate my voxel dimension in the distance calculations?
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
?
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
If you have measurements of any kind you can print them out or pull them to the results table
You’re most welcome
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,
I did attempt to use the smallerConstant
function but sadly also failed.
Would thresholding the distance matrix be an alternative to your suggestion?
Without thresholding the distance matrix we have this
I did mange to pull the averageDistanceOfTouchingNeighbors
to the results table so
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.
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
Keep pushing pixels, I think you’re almost there!
Cheers,
Robert
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
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 . The distance map looks fine.
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:
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
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);
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
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
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.