In [1]:
""" TITLE #######################################################################
# Counting cell in mouse brain
# VVI, New York, 2016-2018 All rights reserved
## DESCRIPTION ##################################################################
## This is a modification od the demo example from MAMBA-Image documentation
## See mamba-image.org for installation and other instructions
#################################################################################
"""
# Importing mamba
from mamba import *
from mambaDisplay.extra import *
import mambaDisplay
import os
import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt
from PIL import Image, ImageMath
Image.MAX_IMAGE_PIXELS = 100000*100000
#
# This function cab 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, figsize=(10, 10)):
pilImage = Mamba2PIL(mbImage)
pilImage = pilImage.convert('P')
fig, ax = plt.subplots(figsize=figsize)
ax.imshow(pilImage)
def ShowImageColor(mbImage, palette, figsize=(10, 10)):
pilImage = Mamba2PIL(mbImage)
pilImage.putpalette(palette)
pilImage = pilImage.convert('P', dither=Image.NONE)
fig, ax = plt.subplots(figsize=figsize)
ax.imshow(pilImage)
In [2]:
# Open files and create images
sourceFolder1 = "N:\\Mouse_cshl\\170613_RMC_MsNT_NoMeOH\\Ch1_512x512\\"
sourceFolder2 = "N:\\Mouse_cshl\\170613_RMC_MsNT_NoMeOH\\Ch2_512x512\\"
destFolder = "N:\\Mouse_cshl\\170613_RMC_MsNT_NoMeOH\\Segmented\\"
palette=mambaDisplay.getPalette("rainbow")
name = mambaDisplay.tagOneColorPalette(255, (255,0,0))
allFiles1 = os.listdir(sourceFolder1)
allFiles2 = os.listdir(sourceFolder2)
im0 = imageMb(sourceFolder1+allFiles1[0])
im1 = imageMb(sourceFolder2+allFiles2[0])
imDepth = im0.getDepth()
imSize = im0.getSize()
print(imDepth, imSize)
In [3]:
# Defining working images.
#imDepth = 8
imWrk1 = imageMb(imSize[0], imSize[1], imDepth)
imWrk2 = imageMb(imSize[0], imSize[1], imDepth)
imIn = imageMb(imSize[0], imSize[1], 8)
imWrk1_8 = imageMb(imSize[0], imSize[1], 8)
imWrk2_8 = imageMb(imSize[0], imSize[1], 8)
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 [4]:
# Set to True to save intermediate steps images
Save_Demo = False
totalNumber = 0
#for i in range(0, len(allFiles2)):
for i in range(0, 1):
if (i > 0):
im0 = imageMb(sourceFolder1+allFiles1[i])
im1 = imageMb(sourceFolder2+allFiles2[i])
fileName, ext = os.path.splitext(destFolder+allFiles2[i])
# Some simple threshold to avoid processing empty (blank) image
#imVolume = computeVolume(im0)
imVolume = 3500
if imVolume > 3000:
'''
# 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: im0.save(fileName+"_initial.png", palette=mambaDisplay.getPalette(name))
convert(im0, imWrk1_8)
ShowImage(imWrk1_8, (5,5))
In [5]:
# Tis background image is to bright to be useful
if Save_Demo: im1.save(fileName+"_backgr.png", palette=mambaDisplay.getPalette(name))
convert(im1, imWrk1_8)
ShowImage(imWrk1_8, (5,5))
In [6]:
# As an alternative, to supress the background, the initial image can be filtered with an alternate filter of size 1.
alternateFilter(im0, imWrk1, 1, False)
if Save_Demo: imWrk1.save(fileName+"_alternate_filter.png", palette=mambaDisplay.getPalette(name))
convert(imWrk1, imWrk1_8)
ShowImage(imWrk1_8, (5,5))
In [7]:
# The minima of this filtered image can be used as markers of the blobs.
convert(imWrk1, imWrk1_8)
highMaxima(imWrk1_8, blobsMarkers, 2)
#ultimateErosion(finalMarkers, blobsMarkers, imWrk3)
#dilate(blobsMarkers, blobsMarkers, 1)
# To give more space between centers of markers and blobs edges, which will make gradients
#erode(blobsMarkers, blobsMarkers, 1)
convert(im0, imIn)
multiSuperpose(imIn, blobsMarkers)
if Save_Demo:
imIn.save(fileName+"_blob_markers.png", palette=mambaDisplay.getPalette(name))
ShowImage(imIn, (5,5))
In [8]:
# Then, we must generate a background marker. This marker can be obtained by
# a marker-controlled watershed of the initial image.
#label(blobsMarkers, imWrk2)
#watershedSegment(imWrk1, imWrk2)
#copyBytePlane(imWrk2, 3, imWrk2_8)
# Another way is to do SKIZ, wich is faster
fastSKIZ(blobsMarkers, backgroundMarker)
negate(backgroundMarker, backgroundMarker)
convert(im0, imIn)
multiSuperpose(imIn, backgroundMarker)
if Save_Demo:
imIn.save(fileName+"_backgr_markers.png", palette=mambaDisplay.getPalette(name))
ShowImageColor(imIn, name, (5,5))
In [9]:
# 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(fileName+"_all_markers.png", palette=mambaDisplay.getPalette(name))
ShowImage(imIn, (5,5))
In [10]:
# The contours of the blobs are obtained by a marker-controlled watershed of
# the gradient image.
#gradient(imWrk2, imWrk1)
#convert(imWrk2, imWrk2_8)
#halfGradient(imWrk2_8, imWrk1_8, "extern")
strongLevelling(imWrk1, imWrk2, 3, True)
convert(imWrk2, imWrk2_8)
ShowImage(imWrk2_8, (5,5))
halfGradient(imWrk2, imWrk1, "extern")
if Save_Demo:
#copy(im0, imIn)
#threshold(imWrk1, blobsContours, 100, 255)
#multiSuperpose(imIn, blobsContours)
imWrk1.save(fileName+"_gradient.png", palette=mambaDisplay.getPalette(name))
convert(imWrk1, imWrk1_8)
ShowImage(imWrk1_8, (5,5))
In [11]:
label(finalMarkers, imWrk2)
watershedSegment(imWrk1, imWrk2)
copyBytePlane(imWrk2, 3, imWrk1_8)
ShowImage(imWrk1_8, (5,5))
threshold(imWrk1_8, finalMarkers, 1, 255)
# Superposing the result to the original image and saving the result.
convert(im0, imIn)
multiSuperpose(imIn, finalMarkers)
if Save_Demo: imIn.save(fileName+"_result_overlay.png", palette=mambaDisplay.getPalette(name))
ShowImage(imIn, (5,5))
In [12]:
# We build the labelled catchment basins
negate(imWrk1_8, imWrk1_8)
copyBytePlane(imWrk2, 0, imWrk2_8)
logic(imWrk1_8, imWrk2_8, imWrk1_8, "inf")
# Then, we obtain the final (and better) result. Each grain is labelled
threshold(imWrk2, finalMarkers, 3, 65535)
convert(finalMarkers, imWrk2_8)
logic(imWrk1_8, imWrk2_8, imWrk1_8, "inf")
if Save_Demo: imWrk1_8.save(fileName+"_result.png", palette=mambaDisplay.getPalette("patchwork"))
ShowImageColor(imWrk1_8, mambaDisplay.getPalette("patchwork"), (5,5), )
