Exercise Solutions¶

1. Creating Synthetic Image¶

Exercise 1: Different slices can be visualised by changing the parameter plot_slice_index in Cell 5, Line 2. You can use this parameter to step-wise move through the image, visualising the edge and whole cross section of spheres, or the zone two spheres overlap.

Exercise 2: Multiple slices can be visualised next to each other using the plot_panels functions with a number of inputs n greater than one. The arguments need to be extend into lists with n entries in accordance. The suggested three slices can be replacing Cell 5, line 3-10 with:

plot_panels(n=3, data_list=[image3d]*3,plot_func=plot_2d_slice_with_values, plot_kwargs_list=[{"slice_index": i} for i in [4, 7, 10]])

Ecercise 3: A sphere can be added to the synthetic image by adding another code block of the same syntax as used for the first two spheres.

# --- 3. Sphere parameter definition ---
centre_s3 = (3, 3, 3)      # (cz, cy, cx)
radius_s3 = 2                 # in voxels
intensity_s3 = 140         # uint8 grayscale

# Add the sphere by calling the function
sphere3 = create_sphere(image3d.shape, centre_s3, radius_s3, intensity_s3)
# Merge with max again
image3d = np.maximum(image3d, sphere3)

The overlap rule is set through the line that merges the sphere into the existing background image image3d = np.maximum(image3d, sphere2). The current default is to prioritise the sphere with the highest intensity value. Another option is to overwrite the background image, i.e. give priority to the sphere added the latest. This can be achieved by replacing the line with

np.copyto(image3d, sphere2, where=(sphere2 != 0)

2. Thresholding¶

Exercise 1: The sphere intensity can be changed through adapting the intensity_s1 variable. Set this to 180 and 190 for the first and second sphere respectively. You will see that this has no significant effect on the Otsu threshold separating background and foreground.

A different scenario emerges, if one of the sphere intensities is moved closer to the background intensity, e.g. if intensity_s1 = 50. Now you can still visually distinguish them, however the found Otsu threshold is 57 and as a result the sphere is not present in the binary image. This is a direct result of the Otsu method maximising the in-between variance of the separated regions. This explains also, why the Otsu threshold produces a good result again once both sphere intensities are moved closer to the background, e.g. intensity_s2 = 50.

Exercise 2: The region of interest can be reduced through cropping the synthetic image/array accordingly. This can be achieved by adding the line in Cell 2.

image3d = image3d[3:11,1:9,0:8] # This crops the image to a region of interest (RoI) containing the first sphere only.

This reduction of the RoI focuses on the second sphere, on bigger images this can be an effective way to reduce computational requirements, and to improve accuracy in finding the optimal separation threshold, as it removes noise and other objects with deviating voxel intensities.

3. Distance Transform¶

Exercise 1: The weights can be changed by adjusting Cell 3, Line weights.append(1.0) for each of the three levels faces, edges, corners, accordingly to 3,4,5.

You can observe that the qualitative pattern of the distance transform are contained. However, the distances are now not corresponding to the absolute distance with respect to image(voxel count) dimension anymore. While this does not have an immediate impact in can lead to problems with the segmentation in the downstream processing steps.

Exercise 2: A histogram showing only the foreground distance values can be plotted through adding a cell with the code:

# Removing all background voxels (value "0")
vals = dt[binary_image != 0].ravel()

# Plot the histogram
plt.figure(figsize=(5,3)); plt.hist(vals, bins=50, color="steelblue"); plt.title("Foreground DT histogram")
plt.xlabel("Distance to background"); plt.ylabel("Count"); plt.tight_layout(); plt.show()

You can see that the highest values do indeed correspond to the set radii of the spheres, and that there are more voxels with lower distances as would be exacted for a spherical shape.

4. Finding seedpoints¶

Exercise 1: You will find that while only two spheres are added to the synthetic image, there are three local extrema identified. In a case where spheres overlap, as it is the case here, the distance map can form plateaus or elongated regions, creating extra extrema that don’t correspond to actual sphere centers. This happens because overlapping areas increase the distance to the background, producing unexpected peaks or valleys. Visualizing slices of the distance transform helps reveal these artifacts. It can for example bee seen that the voxel right at the border/overlap point of the two spheres has a distance value of 2.4, while neighboring voxels have lower distances.

5. Watershed algorithm¶

Exercise 1: The weights of the segmentation voting scheme can be changed by adjusting Cell 5, Line 70 as well as Cell 6, line 27 for each of the three levels faces, edges, corners, accordingly to 5,3,1.

 # Assign larger weights to more direct neighbors
 ...
 w = 5 if nzc == 1 else (3 if nzc == 2 else 1)

or alternatively for equal weights

 w = 5 if nzc == 1 else (3 if nzc == 2 else 1)

You can observe that using weights of 5,3,1 does not alter the segmentation result because the clear step distribution remains intact (similar to 3,2,1). Similarly, using the default synthetic data (two spheres overlapping in one voxel), the segmentation produces the same results when equal weights are applied. However, a look at the distance transform and segmented image reveals this works only in such a special case, as the number of each face, edge, and corner neighbours with one of the two spheres is bigger than the other. In such a case the weighting does not have any effect on the label assignment.

On the other side if we apply the watershed to the example case with a 3rd sphere 'location = 3,3,3, radius = 2, intensity = 140' as build in Exercise 3 of the data creation step, we see that equal weight voting significantly impacts the segmentation. Here the segmentation works correctly when having step weights, however with equal weights, half of the smaller 3rd sphere is grouped together with the bigger neighbouring sphere, it is essentially broken in two pieces. This happens because in one of the overlapping voxels, the number of corner and edge neighbours with the larger sphere is bigger than the number of face neighbours with the added 3rd sphere. As a result this voxel is wrongly assigned to the large sphere, and the incorrect labelling propagates from there into the added sphere.

Exercise 2: The connectivity of voxels can be changed in cell 4 of the step notebook. We first create the full 26-voxel neighbourhood but then remove either remove all but face or all but face and edge neighbours.

For only a 6-voxel (only face neighbours) we add:

offsets = np.array(np.meshgrid([-1, 0, 1], [-1, 0, 1], [-1, 0, 1])).T.reshape(-1, 3)

# Keep only face neighbors (sum of absolute values == 1)
offsets = offsets[np.sum(np.abs(offsets), axis=1) == 1]

Visually inspecting the segmentation, we observe that a restriction to only face neighbours distorts the algorithm. This shows as gaps at the sphere boundary (watershed line), and one sphere growing into the other one. While a general segmentation is still obtained this can be understood as a considerable decrease in segmentation quality.

For an 18-voxel (face and edge neighbours) we add:

# keep only face + edge neighbors (sum of abs values 1 or 2)
face_edge_offsets = offsets[np.sum(np.abs(offsets), axis=1) <= 2]

Taking into account the 18-voxel neighbourhood, we can see that the segmentation result is much more accurate again. However, at closer inspection we find that boundary (watershed line) is offset by one voxel when compared to accounting for the full 26-voxel neighbourhood. In this simple example, this would be still understood as a accurate segmentation. However, when moving to more complex images with less smooth contact/overlap regions, this can introduce inaccuracies that then further propagate along the watershed line.

6. Quantitative Analysis¶

The output table with the quantitative data we obtain is

Label | Volume (voxels) | COM (z, y, x)           | Eq. Diameter (voxels)
----------------------------------------------------------------------
    0 |            3725 | (8, 8, 8)                 | 19.23
    1 |             253 | (7, 5, 4)                 | 7.85
    2 |             118 | (8, 9, 8)                 | 6.09

We can identify region 0 as the background, as it has by far the biggest voxel count, and the CoM is exactly at the centre of the 16x16x16 voxel image. Accordingly, region 1 and 2 must be the spheres. The calculated CoMs correspond exactly with the sphere centres defined in Section 01. Similarly, the calculated equivalent diameters are in good agreement with the defined radii of 4 and 3 voxels. A slight deviation is present, which is a result of the overlapping voxels, which can only be assigned to one of the spheres.