Finding the Interface of two regions of a segmented image-CodePudding

I have a segmented (by watershed) image of two regions that share one boundary. How do I easily find the position of the pixels on the interface? I tried using hints from this answer but could not get it working. Here is my example code:

import matplotlib.pyplot as plt
from scipy import ndimage as ndi

from skimage.segmentation import watershed
from skimage.feature import peak_local_max

from skimage import future
from skimage.measure import label, regionprops, regionprops_table


# Generate an initial image with two overlapping circles
x, y = np.indices((80, 80))
x1, y1, x2, y2 = 28, 28, 44, 52
r1, r2 = 16, 20
mask_circle1 = (x - x1)**2   (y - y1)**2 < r1**2
mask_circle2 = (x - x2)**2   (y - y2)**2 < r2**2
image = np.logical_or(mask_circle1, mask_circle2)

# Now we want to separate the two objects in image
# Generate the markers as local maxima of the distance to the background
distance = ndi.distance_transform_edt(image)
coords = peak_local_max(distance, footprint=np.ones((3, 3)), labels=image)
mask = np.zeros(distance.shape, dtype=bool)
mask[tuple(coords.T)] = True
markers, _ = ndi.label(mask)
labels = watershed(-distance, markers, mask=image)

fig, axes = plt.subplots(ncols=3, figsize=(9, 3), sharex=True, sharey=True)
ax = axes.ravel()

ax[0].imshow(image, cmap=plt.cm.gray)
ax[0].set_title('Overlapping objects')
ax[1].imshow(-distance, cmap=plt.cm.gray)
ax[1].set_title('Distances')
ax[2].imshow(labels, cmap=plt.cm.nipy_spectral)
ax[2].set_title('Separated objects')

for a in ax:
    a.set_axis_off()

fig.tight_layout()
plt.show()

#---------------- find the interface pixels (either of the two interfaces) of these two objects -----------

rag = future.graph.RAG(labels) 
rag.remove_node(0)
for region in regionprops(labels):
    nlist=list(rag.neighbors(region.label))
print(nlist)

The nlist seems to be just a list containing one element 1: [1]. I was expecting position of pixels. I do not have much experience in using the graph and RAG. It seems that rag creates a graph/network of the regions and has the information of which region is next to which one but I cannot extract that information in the form of the interface pixels. Thanks for any help.

CodePudding user response：

Currently the RAG object doesn't keep track of all the regions and boundaries, though we hope to support that in the future. What you found is just the list of adjacent regions.

For now, if you only have two regions, it's not too expensive to do this manually:

from skimage.morphology import dilation

label1 = labels == 1
label2 = labels == 2

boundary = dilation(label1) & dilation(label2)