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")
File number:  2 File Name:  D:\Columns\Selected_Area_10l_Fragments\Not_Yet_Used_Images\Image_14_Plate_3760000240.jpg_17013.37_10459.53_348.00.jpg