Watershed segmentation¶
adaptation for Nissl-stained human cortex¶
Vadim V. Istomin, MD, PhD¶
New York, 2016-2020 All Rights Reserved¶
DESCRIPTION¶
This is a modification of the demo example from MAMBA-Image documentation¶
See mamba-image.org for installation and other instructions¶
Additional comments and task-specific explanations follow below¶
In [109]:
# Importing mamba
from mamba import *
from mambaDisplay.extra import *
import mambaDisplay
import os
import numpy as np
from skimage import exposure
from scipy.optimize import curve_fit
import sklearn.cluster
from statsmodels.nonparametric.kernel_regression import KernelReg
from statsmodels.nonparametric.kde import KDEUnivariate
import matplotlib.pyplot as plt
from matplotlib import cm
%matplotlib inline
from PIL import Image, ImageMath
from PIL import ImageEnhance
from PIL import ImageOps
Image.MAX_IMAGE_PIXELS = 100000*100000
#
# This function can be used to understand how watersheds algorithm works
#
def watershedByGradMinima(imIn, imOut, grid=DEFAULT_GRID):
"""
Computes the watershed line of imIn using its gradient minima as a marker image
"""
grad = imageMb(imIn, 8)
marker = imageMb(imIn, 32)
minbin = imageMb(imIn, 1)
gradient(imIn, grad)
minima(grad, minbin, grid=grid)
label(minbin, marker, grid=grid)
watershedSegment(grad, marker, grid=grid)
copyBytePlane(marker, 3, imOut)
def ShowImage(mbImage, Title="Set Title", figsize=(10, 12)):
pilImage = Mamba2PIL(mbImage)
fig, ax = plt.subplots(figsize=figsize)
plt.title(Title)
ax.imshow(pilImage, cmap=plt.cm.gray)
def ShowImageCrop(mbImage, Title="Set Title", figsize=(10, 12)):
showSize = (512, 512);
imSize = mbImage.getSize()
imDepth = mbImage.getDepth()
mbCrop = imageMb(showSize[0], showSize[1], imDepth)
posXY = ((imSize[0]-512)/2, (imSize[1]-512)/2)
cropCopy(mbImage, posXY, mbCrop, (0,0), showSize)
pilImage = Mamba2PIL(mbCrop)
fig, ax = plt.subplots(figsize=figsize)
plt.title(Title)
ax.imshow(pilImage, cmap=plt.cm.gray)
def ShowImageCrop(mbImage, Title="Set Title", figsize=(10, 12)):
showSize = (512, 512);
imSize = mbImage.getSize()
imDepth = mbImage.getDepth()
mbCrop = imageMb(showSize[0], showSize[1], imDepth)
posXY = ((imSize[0]-512)/2, (imSize[1]-512)/2)
cropCopy(mbImage, posXY, mbCrop, (0,0), showSize)
pilImage = Mamba2PIL(mbCrop)
fig, ax = plt.subplots(figsize=figsize)
plt.title(Title)
ax.imshow(pilImage, cmap=plt.cm.gray)
def ShowPilImage(PilImage, Title="Set Title", figsize=(10, 12)):
fig, ax = plt.subplots(figsize=figsize)
plt.title(Title)
ax.imshow(PilImage, cmap=plt.cm.gray)
def ShowPilImageCrop(PilImage, Title="Set Title", figsize=(10, 12)):
mbImage = imageMb(PilImage.width, PilImage.height, 8)
PIL2Mamba(PilImage, mbImage)
showSize = (512, 512);
mbCrop = imageMb(showSize[0], showSize[1], imDepth)
posXY = ((imSize[0]-512)/2, (imSize[1]-512)/2)
cropCopy(mbImage, posXY, mbCrop, (0,0), showSize)
pilImage = Mamba2PIL(mbCrop)
fig, ax = plt.subplots(figsize=figsize)
plt.title(Title)
ax.imshow(pilImage, cmap=plt.cm.gray)
def ShowImageColor(mbImage, Title="Set Title", figsize=(10, 12), palette = 0):
pilImage = Mamba2PIL(mbImage)
pilImage.putpalette(palette)
fig, ax = plt.subplots(figsize=figsize)
plt.title(Title)
ax.imshow(pilImage, cmap=plt.cm.prism)
def ShowImageCropColor(mbImage, Title="Set Title", figsize=(10, 12), palette = 0):
showSize = (512, 512);
imSize = mbImage.getSize()
imDepth = mbImage.getDepth()
mbCrop = imageMb(showSize[0], showSize[1], imDepth)
posXY = ((imSize[0]-512)/2, (imSize[1]-512)/2)
cropCopy(mbImage, posXY, mbCrop, (0,0), showSize)
pilImage = Mamba2PIL(mbCrop)
pilImage.putpalette(palette)
fig, ax = plt.subplots(figsize=figsize)
plt.title(Title)
ax.imshow(pilImage, cmap='hot')
def ShowImageCropColorXY(mbImage, X, Y, Title="Set Title", figsize=(10, 12), palette = 0):
showSize = (512, 512);
imSize = mbImage.getSize()
imDepth = mbImage.getDepth()
mbCrop = imageMb(showSize[0], showSize[1], imDepth)
posXY = ((imSize[0]-512)/2, (imSize[1]-512)/2)
cropCopy(mbImage, (X, Y), mbCrop, (0,0), showSize)
pilImage = Mamba2PIL(mbCrop)
pilImage.putpalette(palette)
fig, ax = plt.subplots(figsize=figsize)
plt.title(Title)
ax.imshow(pilImage, cmap='hot')
In [88]:
def kde_fit(x, x_grid, bandwidth):
"""Kernel Density Estimation with Scipy"""
# Note that scipy weights its bandwidth by the covariance of the
# input data. To make the results comparable to the other methods,
# we divide the bandwidth by the sample standard deviation here.
#kde = gaussian_kde(x, bw_method=bandwidth/x.std(ddof=1.0))
#kde = gaussian_kde(x, bw_method=50.0)
#kde = gaussian_kde(x, bw_method = 'scott')
kde = KDEUnivariate(x)
kde.fit(bw=51.2)
return kde.evaluate(x_grid)
#return kde.density
def histogram(imIn, Title):
maxRange = computeMaxRange(imIn)
xarray = np.linspace(-maxRange[1]/2, maxRange[1]/2, maxRange[1]+1)
#imWrk = imageMb(imIn)
#imOut = imageMb(imIn)
#imBin = imageMb(imIn, 1)
##threshold(imIn, imBin, 1, 230)
##convert(imBin, imWrk)
##logic(imIn, imWrk, imOut, "inf")
histArray = getHistogram(imIn)
## y_axis values must be normalised
histArray[0] = 0
sum_ys = sum(histArray)
print "Sum_Y=", sum_ys
histArrayNorm = [1.0*_/sum_ys for _ in histArray]
kr = KernelReg(histArrayNorm,xarray,'c')
histSmoothed, histStd = kr.fit(xarray)
y = np.concatenate([histArrayNorm])
#print(y)
histKDE = kde_fit(y, xarray, 5.0)
# Prints numerical values of the histogram
#print(histSmoothed)
plt.plot(xarray+maxRange[1]/2, histArrayNorm, label='Raw')
plt.plot(xarray+maxRange[1]/2, histSmoothed, label='Smoothed')
plt.plot(xarray+maxRange[1]/2, histKDE, label='KDE Fit')
plt.xlabel('Grey Levels')
plt.ylabel('Frequency')
#plt.title(r'$\mathrm{Histogram\ of\ Image}\ \mu=XX,\ \sigma=XX$')
plt.title(Title)
plt.axis([maxRange[0], maxRange[1], 0.0, max(histArrayNorm)])
plt.grid(True)
plt.show()
return histSmoothed
Open files and create images¶
In [89]:
#sourceFolder1 = "C:\\Image_Processing_Mamba\\A19_Reg_No_Repeats\\Croped_1024_1024\\"
sourceFolder1 = "D:\\Columns\\Selected_Area_9m_Fragments\\Not_Yet_Used_Images\\"
#sourceFolder1 = "C:\\Image_Processing_Mamba\\A19_Reg_No_Repeats\\original_frag\images_1\\"
#sourceFolder1 = "N:\\Mouse_cshl\\170613_RMC_MsNT_NoMeOH\\Ch1_512x512\\"
# No background image for human samples
#sourceFolder2 = "N:\\Mouse_cshl\\170613_RMC_MsNT_NoMeOH\\Ch2_512x512\\"
#destFolder = "C:\\Image_Processing_Mamba\\A19_Reg_No_Repeats\\Results_mouse_version\\"
sourceFolder1 = "D:\\Columns\\Selected_Area_10l_Fragments\\Not_Yet_Used_Images\\"
destFolder = "D:\\Columns\\Selected_Area_10l_Fragments\\Segmented_Equalised\\"
destFolderDemo = "D:\\Columns\\Selected_Area_10l_Fragments\\Segmented_Demo\\"
destFolderNormalized = "D:\\Columns\\Selected_Area_10l_Fragments\\Images_Normalized\\"
#destFolder = "N:\\Mouse_cshl\\170613_RMC_MsNT_NoMeOH\\Segmented\\"
palette=mambaDisplay.getPalette("rainbow")
redPalette = mambaDisplay.tagOneColorPalette(255, (255,0,0))
allFiles1 = os.listdir(sourceFolder1)
#allFiles2 = os.listdir(sourceFolder2)
In [90]:
# Set to True to save intermediate steps images
Save_Demo = False
Demo = True
totalNumber = 0
#for i in range(0, len(allFiles1)):
for i in range(2, 3):
#im1 = imageMb(sourceFolder2+allFiles2[0])
im0 = imageMb(sourceFolder1+allFiles1[i])
imDepth = im0.getDepth()
imSize = im0.getSize()
print "File number: ", i, "File Name: ", sourceFolder1+allFiles1[i]
pilImage = Mamba2PIL(im0)
pilImage2 = Mamba2PIL(im0)
#factor = 2.0
#enhancer = ImageEnhance.Sharpness(pilImage)
#enhancer.enhance(factor).save("temp.png")
#pilImage = Image.open("temp.png")
#factor = 3.0
#enhancer = ImageEnhance.Contrast(pilImage)
#enhancer.enhance(factor).save("temp.png")
imIn = imageMb(imSize[0], imSize[1], 8)
if Demo: ShowImage(im0, "Initial Image", (15,15))
### Normalization of initial image
### Mandatory step: adaptation to staining variability
if (True):
img = np.array(pilImage)
#img = io.imread(sourceFolder1+allFiles1[i])
# Contrast stretching
p2, p98 = np.percentile(img, (2, 98))
img_rescale = exposure.rescale_intensity(img, in_range=(p2, p98))
#pilImage = Image.fromarray(img_rescale).convert('L')
pilImage2 = Image.fromarray(img_rescale)
PIL2Mamba(pilImage2, imIn)
if Demo: ShowImage(imIn, "Initial Image Normalized", (15,15))
if Save_Demo: imIn.save(destFolderNormalized+allFiles1[i]+"_initial_normalized.png")
In [91]:
if (Demo): hist1 = histogram(im0, "Initial")
In [92]:
hist1 = histogram(imIn, "Normalized")
bkgrThres = 255.0 - (255.0 - np.argmax(hist1))*2.0
print bkgrThres
Less efficient method of adaptation: histogram equalization¶
In [ ]:
if (False):
#imIn = imageMb("temp.png")
pilImage2 = ImageOps.equalize(pilImage)
PIL2Mamba(pilImage2, imIn)
if Demo: ShowImage(imIn, "Initial Equalized Image", (15,15))
if Save_Demo: imIn.save(destFolderDemo+allFiles1[i]+"_initial.png")
In [93]:
#im1 = imageMb(sourceFolder2+allFiles2[i])
#print(imDepth, imSize)
fileName, ext = os.path.splitext(destFolder+allFiles1[i])
In [94]:
# Defining working images.
#imDepth = 8
imWrk1 = imageMb(imSize[0], imSize[1], imDepth)
imWrk2 = imageMb(imSize[0], imSize[1], imDepth)
imWrk32 = imageMb(imSize[0], imSize[1], 32)
im1 = imageMb(imSize[0], imSize[1], 8)
imWrk1_8 = imageMb(imSize[0], imSize[1], 8)
imWrk2_8 = imageMb(imSize[0], imSize[1], 8)
imBckgrSub = imageMb(imSize[0], imSize[1], 8)
imWrk1_1 = imageMb(imSize[0], imSize[1], 1)
imWrk2_1 = imageMb(imSize[0], imSize[1], 1)
blobsMarkers = imageMb(imSize[0], imSize[1], 1)
backgroundMarker= imageMb(imSize[0], imSize[1], 1)
blobsContours = imageMb(imSize[0], imSize[1], 1)
finalMarkers = imageMb(imSize[0], imSize[1], 1)
In [38]:
# Some simple threshold to avoid processing empty (blank) image
#imVolume = computeVolume(im0)
'''
# Another way to eliminate the backgtound
if Demo: im0.save(fileName+"_initial.png", palette=mambaDisplay.getPalette(name))
strongLevelling(im0, imWrk1, 10, False)
if Demo: imWrk1.save(fileName+"_levelling.png", palette=mambaDisplay.getPalette(name))
floorSub(im0, imWrk1, imWrk2)
if Demo: imWrk2.save(fileName+"_backgr_substr.png", palette=mambaDisplay.getPalette(name))
'''
# Substract background from signal channel (im0, im1)
#sub(im0, im1, im0)
#if Save_Demo: imIn.save(fileName+"_initial.png", palette=redPalette)
#if Demo: ShowImage(imIn, "Initial Image", (15,15))
Out[38]:
In [95]:
# As an alternative, to supress the background, the initial image can be filtered with an alternate filter
negate(imIn, im1)
if Demo: ShowImageCrop(im1, "Initial Negated", (20,15) )
In [98]:
b_Subtracted = False
if (b_Subtracted):
# TRUE
strongLevelling(im1, imWrk2, 12, False)
#largeDodecagonalAlternateFilter(im1, imWrk2, 6, 60, 10, True)
#fullAlternateFilter(im1, imWrk2, 10, True)
if Demo: ShowImageCrop(imWrk2, "Smoothed Background", (20,15) )
sub(im1, imWrk2, imWrk1_8)
if Demo: ShowImageCrop(imWrk1_8, "Smoothed, Background Subtracted", (30,15) )
if Save_Demo: imWrk1_8.save(destFolderDemo+allFiles1[i]+"_bckgr_substr.png")
copy(imWrk1_8, imBckgrSub)
else:
# FALSE
#sub(im1, imWrk2, imWrk1_8)
copy(im1, imWrk1_8)
if Save_Demo: imWrk1_8.save(destFolderDemo+allFiles1[i]+"_bckgr_not_substr.png")
copy(imWrk1_8, imBckgrSub)
if Demo: ShowImageCrop(imWrk1_8, "Smoothed, Background Not Subtracted", (30,15) )
In [99]:
if (Demo):
if (b_Subtracted): hist = histogram(imWrk2, "Smoother")
if (b_Subtracted): hist = histogram(imBckgrSub, "Background Subtracted")
else: hist = histogram(imBckgrSub, "Background Not Subtracted")
In [102]:
# The minima of this filtered image can be used as markers of the blobs.
strongLevelling(imWrk1_8, imWrk1, 5, False)
#fullAlternateFilter(imWrk1_8, imWrk1, 2, False)
if Demo:
ShowImageCrop(imWrk1, "Neurons smoothed by Strong Levelling", (20,20))
In [103]:
if Demo:
hist = histogram(imWrk1, "Neurons smoothed by Strong Levelling")
In [111]:
if (True):
if (b_Subtracted): constSub = 0
else: constSub = 0;
subConst(imWrk1_8, constSub, imWrk2)
highMaxima(imWrk2, blobsMarkers, 50)
#if Demo: ShowImage(blobsMarkers, "High Maxima", (20,20))
convert(blobsMarkers, imWrk2_8)
if Save_Demo: imWrk2_8.save(destFolderDemo+allFiles1[i]+"_markers.png")
copy(imIn, im1)
multiSuperpose(im1, blobsMarkers)
if Demo: ShowImageCropColorXY(im1, 100, 100, "Initial image + Markers", (20,20), palette= redPalette)
if Demo: ShowImageCropColorXY(im1, 100, 1200, "Initial image + Markers", (20,20), palette= redPalette)
In [ ]:
# Another version
if (False):
maxDynamics(imWrk1_8, blobsMarkers, 20)
convert(blobsMarkers, imWrk2_8)
if Demo: ShowImageCrop(imWrk2_8, "Max Dynamics", (20,20))
#if Save_Demo: imWrk1_8.save(destFolderDemo+allFiles1[i]+"_markers.png", palette=mambaDisplay.getPalette("rainbow"))
if Save_Demo: imWrk2_8.save(destFolderDemo+allFiles1[i]+"_markers.png")
copy(imIn, im1)
multiSuperpose(im1, blobsMarkers)
if Demo: ShowImageCropColor(im1, "Initial image + Markers", (20,20), palette= redPalette)
# Yet another
if (False):
threshold(imWrk1, imWrk1_1, 100,255)
if Demo: ShowImage(imWrk1_1, "Threshold", (20,20))
ultimateErosion(imWrk1_1, blobsMarkers, imWrk32)
dilate(blobsMarkers, blobsMarkers, 5)
if Demo: ShowImage(blobsMarkers, "Maxima by ultimate erosion", (20,20))
In [75]:
# Background marker
fastSKIZ(blobsMarkers, backgroundMarker)
negate(backgroundMarker, backgroundMarker)
convert(im0, imIn)
multiSuperpose(imIn, backgroundMarker)
if Save_Demo:
imIn.save(destFolderDemo+allFiles1[i]+"_backgr_markers.png")
if Demo: ShowImageCropColor(imIn, "SKIZ", (30,15), palette=palette)
In [76]:
# The two sets of markers are merged. We must however make sure that they are
# separated by at least one pixel.
dilate(backgroundMarker, blobsContours)
diff(blobsMarkers, blobsContours, blobsMarkers)
logic(blobsMarkers, backgroundMarker, finalMarkers, "sup")
convert(im0, imIn)
multiSuperpose(imIn, finalMarkers)
if Save_Demo:
imIn.save(destFolderDemo+allFiles1[i]+"_all_markers.png")
if Demo: ShowImageCropColor(imIn, "Blobs and Background Markers", (30,15), palette=redPalette)
In [77]:
#strongLevelling(imBckgrSub, imWrk1, 3, True)
#gradient(imWrk1, imWrk2, 1)
negate(im0, imIn)
gradient(imIn, imWrk2, 1)
#halfGradient(imWrk1, imWrk2, 1)
#regularisedGradient(imBckgrSub, imWrk1, 1)
if Save_Demo:
#copy(im0, imIn)
#threshold(imWrk1, blobsContours, 100, 255)
#multiSuperpose(imIn, blobsContours)
imWrk2.save(destFolderDemo+allFiles1[i]+"_gradient.png", palette=redPalette)
if Demo: ShowImageCrop(imWrk2, "Gradient", (30,15))
In [78]:
imWrk32 = imageMb(imSize[0], imSize[1], 32)
#imDepth1 = finalMarkers.getDepth()
#imDepth2 = imWrk2.getDepth()
#print imDepth1, imDepth2
label(finalMarkers, imWrk32)
# imWrk32 is output image
watershedSegment(imWrk2, imWrk32)
# byte plane 3 is used for contours
copyBytePlane(imWrk32, 3, imWrk1_8)
if Demo: ShowImageCrop(imWrk1_8, "Segmentation contours", (20,20))
In [79]:
threshold(imWrk1_8, finalMarkers, 1, 255)
closeHoles(finalMarkers, imWrk2_1)
if (False):
erode(imWrk2_1, imWrk1_1, 3)
build(imWrk2_1, imWrk1_1)
if (True):
copy(imWrk2_1, imWrk1_1)
imWrk1_1.save(fileName+"_segmented.png")
if Save_Demo: imWrk1_1.save(destFolderDemo+allFiles1[i]+"_segmented.png")
# Superposing the result to the original image and saving the result.
convert(im0, imIn)
multiSuperpose(imIn, imWrk1_1)
#if Save_Demo:
if Save_Demo: imIn.save(destFolderDemo+allFiles1[i]+"_result_overlay.png", palette=redPalette)
if Demo: ShowImageCropColor(imIn, "Segmentation Overlay", (20,20), palette = palette)
In [80]:
# We build the labelled catchment basins
# grey contours negated
negate(imWrk1_8, imWrk1_8)
# first byte plane with grey markers
copyBytePlane(imWrk32, 0, imWrk2_8)
if Demo: ShowImageCrop(imWrk2_8, "Bite Plane", (30,15)) #, mambaDisplay.getPalette("patchwork"))
logic(imWrk1_8, imWrk2_8, imWrk1_8, "inf")
#if Demo: ShowImageCrop(imWrk1_8, "All bassins collected", (30,15)) # mambaDisplay.getPalette("patchwork"))
In [82]:
# Then, we obtain the final (and better) result. Each grain is "area labelled"
threshold(imWrk32, finalMarkers, 3, 65535)
areaLabelling(finalMarkers, imWrk32)
copyBytePlane(imWrk32, 0, imWrk2_8)
if Demo: ShowImageColor(imWrk2_8, "Bassins area labeled", (30,15), mambaDisplay.getPalette("patchwork"))
if Demo: ShowImageCropColor(imWrk2_8, "Bassins area labeled", (30,15), mambaDisplay.getPalette("patchwork"))
In [61]:
build(finalMarkers, blobsMarkers)
convert(finalMarkers, imWrk2_8)
logic(imWrk1_8, imWrk2_8, imWrk1_8, "inf") # Reconstruct blob markers to eliminate segmented particles that belong to background (incorrect SKIZ)
if Demo: ShowImageCropColor(imWrk1_8, "Final", (30,15), mambaDisplay.getPalette("patchwork"))
In [62]:
if Demo: ShowImageColor(imWrk1_8, "Final", (30,15), mambaDisplay.getPalette("patchwork"))
if Save_Demo: imWrk1_8.save(destFolderDemo+allFiles1[i]+"_final_result_color.png", mambaDisplay.getPalette("patchwork"))
if Save_Demo: imWrk1_8.save(destFolderDemo+allFiles1[i]+"_final_result_grey.png")
In [ ]:
print "Done"