CLIJ2

Logo

GPU accelerated image processing for everyone

CLIJ2 home

Tribolium embryo morphometry

Authors: Robert Haase, Daniela Vorkel, April 2020

Source

This script is an example of heavy GPU-accelerated processing. It is recommended to use a dedicated graphics card with at least 8 GB of GDDR6 memory. Otherwise, it may be quite slow.

Let’s initialize that graphics card and mesure the start time.

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

run("Close All");
time = getTime();
Ext.CLIJ2_startTimeTracing();

Load a data set

The dataset is available online. It shows a Tribolium castaneum embryo, imaged by a custom light sheet microscope, at a wavelength of 488nm (Imaging credits: Daniela Vorkel, Myers lab, MPI CBG). The data set has been resampled to a voxel size of 1x1x1 microns. The embryo expresses nuclei-GFP. We will use the dataset to detect nuclei and to generate an estimated cell-segmentation.

All processing steps are performed in 3D space. For visualization purpose, we are using the maximum intensity projection in Z:

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

print("Loading took " + (getTime() - time) + " msec");

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

// visualize the dataset
show(input, "input");

> Loading took 265 msec

CLIJ2_maximumZProjection_result97

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");

CLIJ2_maximumZProjection_result101

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");

CLIJ2_maximumZProjection_result106

Let’s see how many spots are left:

Ext.CLIJ2_getMaximumOfAllPixels(labelled_spots, number_of_spots);
print("Number of detected spots: " + number_of_spots);

> Number of detected spots: 1131

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");
	}
}

CLIJ2_maximumZProjection_result109 CLIJ2_maximumZProjection_result110 CLIJ2_maximumZProjection_result111 CLIJ2_maximumZProjection_result112 CLIJ2_maximumZProjection_result113

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");

CLIJ2_maximumZProjection_result116

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_generateDistanceMatrix(pointlist, pointlist, distance_matrix);
show(distance_matrix, "distance matrix");

CLIJ2_maximumZProjection_result119

Starting from the label map of segmented cells, we generate a touch matrix:

Ext.CLIJ2_generateTouchMatrix(labels, touch_matrix);

// touch matrix:
show_spots(touch_matrix, "touch matrix");

CLIJ2_maximumZProjection_result122

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);

CLIJ2_maximumZProjection_result125

Quantitative analysis of distance between neighbors

Next, we determine the average 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);
// we replace NaN values with zeros so that later maximum-projections work.
Ext.CLIJ2_copy(distances_vector, distances_vector1);
Ext.CLIJ2_undefinedToZero(distances_vector1, distances_vector);

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

CLIJ2_maximumZProjection_result129

Now, we measure the mean between neighbors and visualize it as above.

Ext.CLIJ2_meanOfTouchingNeighbors(distances_vector, touch_matrix, local_mean_distances_vector);
Ext.CLIJ2_replaceIntensities(labels, local_mean_distances_vector, local_mean_pixel_count_map);
show(local_mean_pixel_count_map, "neighbor mean distance map");
run("Fire");
setMinAndMax(0, 50);

CLIJ2_maximumZProjection_result132

We can also use the minimum, median and maximum to measure distances:

Ext.CLIJ2_minimumOfTouchingNeighbors(distances_vector, touch_matrix, local_minimum_distances_vector);
Ext.CLIJ2_replaceIntensities(labels, local_minimum_distances_vector, local_minimum_pixel_count_map);
show(local_minimum_pixel_count_map, "neighbor minimum distance map");
run("Fire");
setMinAndMax(0, 50);

Ext.CLIJ2_medianOfTouchingNeighbors(distances_vector, touch_matrix, local_median_distances_vector);
Ext.CLIJ2_replaceIntensities(labels, local_median_distances_vector, local_median_pixel_count_map);
show(local_median_pixel_count_map, "neighbor median distance map");
run("Fire");
setMinAndMax(0, 50);

Ext.CLIJ2_maximumOfTouchingNeighbors(distances_vector, touch_matrix, local_maximum_distances_vector);
Ext.CLIJ2_replaceIntensities(labels, local_maximum_distances_vector, local_maximum_pixel_count_map);
show(local_maximum_pixel_count_map, "neighbor maximum distance map");
run("Fire");
setMinAndMax(0, 50);

Ext.CLIJ2_standardDeviationOfTouchingNeighbors(distances_vector, touch_matrix, local_stddev_distances_vector);
Ext.CLIJ2_replaceIntensities(labels, local_stddev_distances_vector, local_stddev_pixel_count_map);
show(local_stddev_pixel_count_map, "neighbor standard deviation distance map");
run("Fire");
setMinAndMax(0, 50);

CLIJ2_maximumZProjection_result135 CLIJ2_maximumZProjection_result138 CLIJ2_maximumZProjection_result141 CLIJ2_maximumZProjection_result144

Performance evaluation

Finally, a time measurement. Note that performing this workflow in ImageJ macro markdown is slower, because intermediate results are saved to disc.

print("The whole workflow took " + (getTime() - time) + " msec");


> The whole workflow took 4245 msec

Detailed time tracing for all operations

Ext.CLIJ2_stopTimeTracing();
Ext.CLIJ2_getTimeTracing(time_traces);
print(time_traces);


> > timeTracing
>  > MaximumZProjection
>  < MaximumZProjection           2.9759 ms
>  > Copy
>  < Copy                         9.0152 ms
>  > GaussianBlur3D
>  < GaussianBlur3D               44.7836 ms
>  > DetectMaximaBox
>   > Mean3DBox
>   < Mean3DBox                   37.0359 ms
>  < DetectMaximaBox              53.6302 ms
>  > Maximum3DBox
>   > Copy
>   < Copy                        8.4579 ms
>  < Maximum3DBox                 35.7375 ms
>  > MaximumZProjection
>  < MaximumZProjection           1.5405 ms
>  > Threshold
>   > GreaterOrEqualConstant
>   < GreaterOrEqualConstant      7.437 ms
>  < Threshold                    7.449 ms
>  > Mask
>  < Mask                         7.8011 ms
>  > LabelSpots
>   > Set
>   < Set                         7.2574 ms
>   > SumXProjection
>   < SumXProjection              4.1411 ms
>   > SumYProjection
>   < SumYProjection              0.6554 ms
>  < LabelSpots                   24.7585 ms
>  > Maximum3DBox
>   > Copy
>   < Copy                        7.7649 ms
>  < Maximum3DBox                 32.5717 ms
>  > MaximumZProjection
>  < MaximumZProjection           1.173 ms
>  > GetMaximumOfAllPixels
>   > MaximumOfAllPixels
>    > MaximumZProjection
>    < MaximumZProjection         1.283 ms
>    > MaximumYProjection
>    < MaximumYProjection         0.3715 ms
>    > MaximumXProjection
>    < MaximumXProjection         0.3312 ms
>   < MaximumOfAllPixels          2.9054 ms
>  < GetMaximumOfAllPixels        2.9123 ms
>  > Copy
>  < Copy                         7.2455 ms
>  > OnlyzeroOverwriteMaximumBox
>  < OnlyzeroOverwriteMaximumBox  11.9568 ms
>  > OnlyzeroOverwriteMaximumDiamond
>  < OnlyzeroOverwriteMaximumDiamond2.5826 ms
>  > MaximumZProjection
>  < MaximumZProjection           1.253 ms
>  > OnlyzeroOverwriteMaximumBox
>  < OnlyzeroOverwriteMaximumBox  5.0572 ms
>  > OnlyzeroOverwriteMaximumDiamond
>  < OnlyzeroOverwriteMaximumDiamond2.5588 ms
>  > OnlyzeroOverwriteMaximumBox
>  < OnlyzeroOverwriteMaximumBox  4.7444 ms
>  > OnlyzeroOverwriteMaximumDiamond
>  < OnlyzeroOverwriteMaximumDiamond2.6358 ms
>  > MaximumZProjection
>  < MaximumZProjection           1.5018 ms
>  > OnlyzeroOverwriteMaximumBox
>  < OnlyzeroOverwriteMaximumBox  5.0227 ms
>  > OnlyzeroOverwriteMaximumDiamond
>  < OnlyzeroOverwriteMaximumDiamond2.5737 ms
>  > OnlyzeroOverwriteMaximumBox
>  < OnlyzeroOverwriteMaximumBox  4.7476 ms
>  > OnlyzeroOverwriteMaximumDiamond
>  < OnlyzeroOverwriteMaximumDiamond2.5303 ms
>  > MaximumZProjection
>  < MaximumZProjection           1.2085 ms
>  > OnlyzeroOverwriteMaximumBox
>  < OnlyzeroOverwriteMaximumBox  4.7258 ms
>  > OnlyzeroOverwriteMaximumDiamond
>  < OnlyzeroOverwriteMaximumDiamond2.478 ms
>  > OnlyzeroOverwriteMaximumBox
>  < OnlyzeroOverwriteMaximumBox  4.663 ms
>  > OnlyzeroOverwriteMaximumDiamond
>  < OnlyzeroOverwriteMaximumDiamond2.5056 ms
>  > MaximumZProjection
>  < MaximumZProjection           1.3752 ms
>  > OnlyzeroOverwriteMaximumBox
>  < OnlyzeroOverwriteMaximumBox  4.6049 ms
>  > OnlyzeroOverwriteMaximumDiamond
>  < OnlyzeroOverwriteMaximumDiamond2.5639 ms
>  > OnlyzeroOverwriteMaximumBox
>  < OnlyzeroOverwriteMaximumBox  4.6176 ms
>  > OnlyzeroOverwriteMaximumDiamond
>  < OnlyzeroOverwriteMaximumDiamond2.5136 ms
>  > MaximumZProjection
>  < MaximumZProjection           1.1664 ms
>  > OnlyzeroOverwriteMaximumBox
>  < OnlyzeroOverwriteMaximumBox  4.4245 ms
>  > OnlyzeroOverwriteMaximumDiamond
>  < OnlyzeroOverwriteMaximumDiamond2.5655 ms
>  > Threshold
>   > GreaterOrEqualConstant
>   < GreaterOrEqualConstant      7.9888 ms
>  < Threshold                    8.0028 ms
>  > ErodeBox
>  < ErodeBox                     3.1491 ms
>  > ErodeBox
>  < ErodeBox                     2.9736 ms
>  > ErodeBox
>  < ErodeBox                     2.8717 ms
>  > ErodeBox
>  < ErodeBox                     2.7934 ms
>  > ErodeBox
>  < ErodeBox                     2.7161 ms
>  > ErodeBox
>  < ErodeBox                     2.7864 ms
>  > ErodeBox
>  < ErodeBox                     2.593 ms
>  > ErodeBox
>  < ErodeBox                     2.5627 ms
>  > Mask
>  < Mask                         7.6373 ms
>  > MaximumZProjection
>  < MaximumZProjection           1.6464 ms
>  > LabelledSpotsToPointList
>  < LabelledSpotsToPointList     1.3793 ms
>  > GenerateDistanceMatrix
>  < GenerateDistanceMatrix       1.9646 ms
>  > MaximumZProjection
>  < MaximumZProjection           1.1825 ms
>  > GenerateTouchMatrix
>   > Set
>   < Set                         1.0726 ms
>  < GenerateTouchMatrix          2.9299 ms
>  > Maximum3DBox
>  < Maximum3DBox                 2.2232 ms
>  > MaximumZProjection
>  < MaximumZProjection           0.8365 ms
>  > MultiplyImages
>  < MultiplyImages               1.0607 ms
>  > GetDimensions
>  < GetDimensions                0.0028 ms
>  > TouchMatrixToMesh
>  < TouchMatrixToMesh            7.4039 ms
>  > MaximumZProjection
>  < MaximumZProjection           0.7299 ms
>  > SetColumn
>  < SetColumn                    0.4095 ms
>  > AverageDistanceOfTouchingNeighbors
>  < AverageDistanceOfTouchingNeighbors0.759 ms
>  > Copy
>  < Copy                         0.2599 ms
>  > UndefinedToZero
>  < UndefinedToZero              0.2419 ms
>  > ReplaceIntensities
>  < ReplaceIntensities           7.5108 ms
>  > MaximumZProjection
>  < MaximumZProjection           0.6875 ms
>  > MeanOfTouchingNeighbors
>  < MeanOfTouchingNeighbors      0.8507 ms
>  > ReplaceIntensities
>  < ReplaceIntensities           7.505 ms
>  > MaximumZProjection
>  < MaximumZProjection           0.6716 ms
>  > MinimumOfTouchingNeighbors
>  < MinimumOfTouchingNeighbors   1.1877 ms
>  > ReplaceIntensities
>  < ReplaceIntensities           7.9356 ms
>  > MaximumZProjection
>  < MaximumZProjection           0.8287 ms
>  > MedianOfTouchingNeighbors
>  < MedianOfTouchingNeighbors    1.0459 ms
>  > ReplaceIntensities
>  < ReplaceIntensities           7.8399 ms
>  > MaximumZProjection
>  < MaximumZProjection           0.7293 ms
>  > MaximumOfTouchingNeighbors
>  < MaximumOfTouchingNeighbors   0.9384 ms
>  > ReplaceIntensities
>  < ReplaceIntensities           7.7835 ms
>  > MaximumZProjection
>  < MaximumZProjection           0.7976 ms
>  > ReplaceIntensities
>  < ReplaceIntensities           7.8999 ms
>  > MaximumZProjection
>  < MaximumZProjection           0.6905 ms
> < timeTracing                   4249.9027 ms
>

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();


> GPU contains 28 images.
> - CLIJ2_threshold_result102[net.haesleinhuepf.clij.clearcl.ClearCLPeerPointer@34952cde] 204.8 Mb
> - CLIJ2_maximumOfTouchingNeighbors_result139[net.haesleinhuepf.clij.clearcl.ClearCLPeerPointer@553a28fd] 4.4 kb
> - CLIJ2_generateTouchMatrix_result120[net.haesleinhuepf.clij.clearcl.ClearCLPeerPointer@68258ac7] 4.9 Mb
> - CLIJ2_mask_result103[net.haesleinhuepf.clij.clearcl.ClearCLPeerPointer@22f19cba] 204.8 Mb
> - CLIJ2_minimumOfTouchingNeighbors_result133[net.haesleinhuepf.clij.clearcl.ClearCLPeerPointer@3f9fa7e6] 4.4 kb
> - CLIJ2_onlyzeroOverwriteMaximumBox_result108[net.haesleinhuepf.clij.clearcl.ClearCLPeerPointer@3be61d28] 204.8 Mb
> - CLIJ2_standardDeviationOfTouchingNeighbors_result142[net.haesleinhuepf.clij.clearcl.ClearCLPeerPointer@7fd1db21] 4.4 kb
> - CLIJ2_gaussianBlur3D_result98[net.haesleinhuepf.clij.clearcl.ClearCLPeerPointer@51eb8c2d] 204.8 Mb
> - CLIJ2_replaceIntensities_result137[net.haesleinhuepf.clij.clearcl.ClearCLPeerPointer@1389b161] 204.8 Mb
> - CLIJ2_detectMaximaBox_result99[net.haesleinhuepf.clij.clearcl.ClearCLPeerPointer@188df165] 204.8 Mb
> - CLIJ2_replaceIntensities_result134[net.haesleinhuepf.clij.clearcl.ClearCLPeerPointer@69743986] 204.8 Mb
> - CLIJ2_copy_result107[net.haesleinhuepf.clij.clearcl.ClearCLPeerPointer@64f35a61] 204.8 Mb
> - CLIJ2_copy_result127[net.haesleinhuepf.clij.clearcl.ClearCLPeerPointer@3a2fed27] 4.4 kb
> - CLIJ2_replaceIntensities_result131[net.haesleinhuepf.clij.clearcl.ClearCLPeerPointer@4dde25de] 204.8 Mb
> - CLIJ2_threshold_result114[net.haesleinhuepf.clij.clearcl.ClearCLPeerPointer@5ae9dea] 204.8 Mb
> - CLIJ2_medianOfTouchingNeighbors_result136[net.haesleinhuepf.clij.clearcl.ClearCLPeerPointer@7b70ff6c] 4.4 kb
> - CLIJ2_create3D_result124[net.haesleinhuepf.clij.clearcl.ClearCLPeerPointer@3ef363f0] 204.8 Mb
> - lund1051_resampled.tif[net.haesleinhuepf.clij.clearcl.ClearCLPeerPointer@2bf5125] 204.8 Mb
> - CLIJ2_multiplyImages_result123[net.haesleinhuepf.clij.clearcl.ClearCLPeerPointer@730fdd55] 4.9 Mb
> - CLIJ2_labelSpots_result104[net.haesleinhuepf.clij.clearcl.ClearCLPeerPointer@f78d7c7] 204.8 Mb
> - CLIJ2_mask_result115[net.haesleinhuepf.clij.clearcl.ClearCLPeerPointer@2ab1d93b] 204.8 Mb
> - CLIJ2_averageDistanceOfTouchingNeighbors_result126[net.haesleinhuepf.clij.clearcl.ClearCLPeerPointer@527f9331] 4.4 kb
> - CLIJ2_meanOfTouchingNeighbors_result130[net.haesleinhuepf.clij.clearcl.ClearCLPeerPointer@5f17195e] 4.4 kb
> - CLIJ2_replaceIntensities_result128[net.haesleinhuepf.clij.clearcl.ClearCLPeerPointer@16183b94] 204.8 Mb
> - CLIJ2_labelledSpotsToPointList_result117[net.haesleinhuepf.clij.clearcl.ClearCLPeerPointer@32c4506d] 13.3 kb
> - CLIJ2_replaceIntensities_result143[net.haesleinhuepf.clij.clearcl.ClearCLPeerPointer@1c0cab1f] 204.8 Mb
> - CLIJ2_replaceIntensities_result140[net.haesleinhuepf.clij.clearcl.ClearCLPeerPointer@7b369c29] 204.8 Mb
> - CLIJ2_generateDistanceMatrix_result118[net.haesleinhuepf.clij.clearcl.ClearCLPeerPointer@2ca5dc40] 4.9 Mb
> = 3.4 Gb
>

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);
}

</pre>