NeuroXplorer

Introductory Python Watershed Process Documentation

Augusto Ramirez, Thomas Keady, Albert Lee

Part 1 : Data collection

Data collection tutorials can be found at https://github.com/openconnectome/ndio-demos.git. Detailed documentation also exists at https://github.com/openconnectome/ndio. Data collected from ndio and saved as a numpy file with the name 'nXp_data.npy' can be run through our watershed process when placed in the same folder .

The data used for our watershed process can be found at https://www.dropbox.com/s/95rsvfdhkiodq1x/nXp_data.npy?dl=0. The rest of the necessary programs are located on our github repo at https://github.com/Connectomics-Classes/NeuroXplorer to compare results.

Part 2 : Watershed

The complete watershed file can be found in our github repo in the folder 'testing' under the name 'watershed_test.py'. As long as the 'nxp_data.np' file is in the same folder, the program can be immediately run from terminal using the command 'python watershed_test.py'. The following steps deconstruct the program on a very basic level.

Step 1:

First launch python and import the following packages. If an error states that a certain module does not exist simply use 'pip install insert package name here' in order to download the package. Detailed instructions on using pip can be found at https://pip.pypa.io/en/stable/installing/.


In [ ]:
### Watershed segmentation imports
import numpy as np
import matplotlib.pyplot as plt
from scipy import ndimage as ndi

from skimage.morphology import watershed
from skimage.feature import peak_local_max
from skimage.morphology import closing
from skimage.morphology import binary_dilation
from skimage.morphology import binary_erosion
from skimage.filters.rank import threshold
from skimage.exposure import rescale_intensity
from skimage.filters import threshold_otsu
###

print "Done importing packages"

Step 2

Next load the necessary data from our 'nxp_data.npy file'. If you are using your own data, make sure you have already replaced our npy file with your own. In any case, please make sure the file you are using is named 'nxp_data.np'.


In [ ]:
# Gather images from nXp_data.npy file
print "Getting images"
membrane_images = np.load('nXp_data.npy')

# Assigning one image from array to variable pre_image
pre_image = membrane_images[0]

Step 3

Several different image processing tools are necessary to increase precision.


In [ ]:
########
# Binary closing on image
########
print "Running closing"

close_image = closing(pre_image, selem=np.ones((4,4)))

print "Running median filter"
# Remove white speckles of noise
img_med = ndi.median_filter(close_image, size=4)

print "Running threshold"
#binary_img = threshold(close_image, selem=np.ones((200,200)))
threshold = threshold_otsu(img_med)
binary_img = img_med > threshold

image = binary_img
foreground = 1 - image

# image = foreground

Step 4

The below code is for the actual watershed segmentation. The distance parameters can be modified to produce different results.


In [ ]:
########
# Watershed segmentation
########
print "Running watershed"

# 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(foreground) # parameter: sampling=[100,100])
distance_float = rescale_intensity(distance, in_range='image', out_range='float')

local_maxi = peak_local_max(distance, indices=False, footprint=np.ones((10, 10)),
                            labels=foreground)
markers = ndi.label(local_maxi)[0] #[0]

#markers = distance > 5
#print markers
#markers = threshold(distance_float, selem=np.ones((10,10)))
labels =  watershed(foreground, markers, mask=foreground)
#labels =  watershed(distance, markers, mask=image)

fig, axes = plt.subplots(ncols=3, figsize=(10, 10), sharex=True, sharey=True, subplot_kw={'adjustable':'box-forced'})
ax0, ax1, ax2 = axes

ax0.imshow(image, cmap=plt.cm.gray, interpolation='nearest')
ax0.set_title('Overlapping objects')
ax1.imshow(distance, cmap=plt.cm.jet, interpolation='nearest')
ax1.set_title('Distances')
ax2.imshow(labels, cmap=plt.cm.spectral, interpolation='nearest')
ax2.set_title('Separated objects')

for ax in axes:
    ax.axis('off')

Step 5

Finally a plot of the annotated image should pop up.


In [ ]:
fig.subplots_adjust(hspace=0.01, wspace=0.01, top=0.9, bottom=0, left=0,
                    right=1)

plt.show()

Step 6

Lastly, evaluate the data to test your parameters.


In [ ]:
########
# Evaluation
########
print "Running evaluation"

print "...Split VI Evaluation" # Want closest to [0,0] for (undersegmentation (false_merges), oversegmentation (false_splits))
split_eval = ev.split_vi(ws, membrane_ground[0])
print split_eval

print "...Adjusted Rand Index Evaluation" # Want the closest to 1.0
adjusted_rand = ev.adj_rand_index(ws, membrane_ground[0])
print adjusted_rand

print "...Fowlkes-Mallows Index Evaluation" # Want closest to 1
fm = ev.fm_index(ws, membrane_ground[0])
print fm