notebook.community

Edit and run

Chan-Vese segmentation algorithm demonstration

Explore the use of Chan-Vese active contours for image segmentation on a small sample of natural and medical images

Author: Jacob Reinhold

Setup notebook





In [108]:



    


import os
import sys

from IPython.display import HTML
import matplotlib.pyplot as plt
from matplotlib import animation, rc
import numpy as np
from PIL import Image
from scipy.ndimage import distance_transform_edt as distance
from skimage import data, img_as_float
from skimage.color import rgb2gray
from skimage.draw import circle_perimeter
from skimage.filters import gaussian
from skimage.segmentation import chan_vese, active_contour

Support in-notebook plotting





In [2]:



    


%matplotlib inline

Report versions





In [3]:



    


print('numpy version: {}'.format(np.__version__))
from matplotlib import __version__ as mplver
print('matplotlib version: {}'.format(mplver))



    


















    






numpy version: 1.15.4
matplotlib version: 3.0.2

















In [4]:



    


pv = sys.version_info
print('python version: {}.{}.{}'.format(pv.major, pv.minor, pv.micro))



    


















    






python version: 3.7.1

Reload packages where content for package development





In [5]:



    


%load_ext autoreload
%autoreload 2

Define helper functions





In [51]:



    


def plot_results(img, cv, rotn=3, figsize=(8,8)):
    fig, axes = plt.subplots(2, 2, figsize=figsize)
    ax = axes.flatten()

    ax[0].imshow(np.rot90(img,rotn), cmap="gray")
    ax[0].set_axis_off()
    ax[0].set_title("Original Image", fontsize=12)

    ax[1].imshow(np.rot90(cv[0],rotn), cmap="gray")
    ax[1].set_axis_off()
    title = "Chan-Vese segmentation - {} iterations".format(len(cv[2]))
    ax[1].set_title(title, fontsize=12)

    ax[2].imshow(np.rot90(cv[1],rotn), cmap="gray")
    ax[2].set_axis_off()
    ax[2].set_title("Final Level Set", fontsize=12)

    ax[3].plot(np.array(cv[2])/1000)
    ax[3].set_title("Evolution of energy over iterations", fontsize=12)

    fig.tight_layout()

    
def animate(image, segs, interval=20):
    fig, ax = plt.subplots(figsize=(8, 8))
    ax.imshow(image, cmap='gray');
    im = ax.imshow(segs[0], alpha=0.5, cmap='inferno');
    ax.axis('off')

    def init():
        im.set_data(segs[0])
        return [im]

    def animate(i):
        im.set_array(segs[i])
        return [im]

    anim = animation.FuncAnimation(fig, animate, init_func=init,
                                   frames=len(segs), interval=20, blit=True);
    return anim

Compare and contrast edge-based active contour vs. Chan-Vese

On Gaussian image





In [118]:



    


img = np.zeros((100, 100))
rr, cc = circle_perimeter(35, 45, 25)
img[rr, cc] = 1
img = gaussian(img, 10)

Edge-based

Kass, M.; Witkin, A.; Terzopoulos, D. “Snakes: Active contour models”. International Journal of Computer Vision 1 (4): 321 (1988). DOI:10.1007/BF00133570





In [119]:



    


s = np.linspace(0, 2*np.pi, 100)
init = 50 * np.array([np.cos(s), np.sin(s)]).T + 50
snake = active_contour(img, init, w_edge=0, w_line=1)



    











In [120]:



    


fig, ax = plt.subplots(figsize=(7, 7))
ax.imshow(img, cmap=plt.cm.gray)
ax.plot(init[:, 0], init[:, 1], '--r', lw=3)
ax.plot(snake[:, 0], snake[:, 1], '-b', lw=3)
ax.set_xticks([]), ax.set_yticks([])
ax.axis([0, img.shape[1], img.shape[0], 0]);

Chan-Vese





In [121]:



    


init = small_disk(img.shape, (49, 49), 40)
plt.imshow(init,cmap='gray'); plt.axis('off'); plt.title('Initial level set');



    


















    
























In [122]:



    


cv = chan_vese(img, mu=0.07, lambda1=1, lambda2=1, tol=1e-3, max_iter=500,
               dt=0.5, init_level_set=init, extended_output=True)



    











In [123]:



    


plot_results(img, cv, 0, (10,10))

On Gaussian image with uniform noise in [0,1]





In [130]:



    


img = np.zeros((100, 100))
rr, cc = circle_perimeter(35, 45, 25)
img[rr, cc] = 1
img = gaussian(img, 10)
img += np.random.rand(*img.shape) * 0.01

Edge-based





In [131]:



    


s = np.linspace(0, 2*np.pi, 100)
init = 50 * np.array([np.cos(s), np.sin(s)]).T + 50
snake = active_contour(img, init, w_edge=0, w_line=1)



    











In [132]:



    


fig, ax = plt.subplots(figsize=(7, 7))
ax.imshow(img, cmap=plt.cm.gray)
ax.plot(init[:, 0], init[:, 1], '--r', lw=3)
ax.plot(snake[:, 0], snake[:, 1], '-b', lw=3)
ax.set_xticks([]), ax.set_yticks([])
ax.axis([0, img.shape[1], img.shape[0], 0]);

Chan-Vese





In [133]:



    


init = small_disk(img.shape, (49, 49), 40)
plt.imshow(init,cmap='gray'); plt.axis('off'); plt.title('Initial level set');



    


















    
























In [136]:



    


cv = chan_vese(img, mu=0.07, lambda1=1, lambda2=1, tol=1e-3, max_iter=5000,
               dt=0.5, init_level_set=init, extended_output=True)



    











In [137]:



    


plot_results(img, cv, 0, (10,10))

Test edge-based active contour on natural images





In [7]:



    


img = data.astronaut()
img = rgb2gray(img)

s = np.linspace(0, 2*np.pi, 400)
x = 220 + 100*np.cos(s)
y = 100 + 100*np.sin(s)
init = np.array([x, y]).T

snake = active_contour(gaussian(img, 3),
                       init, alpha=0.015, beta=10, gamma=0.001)



    











In [9]:



    


fig, ax = plt.subplots(figsize=(7, 7))
ax.imshow(img, cmap=plt.cm.gray)
ax.plot(init[:, 0], init[:, 1], '--r', lw=3)
ax.plot(snake[:, 0], snake[:, 1], '-b', lw=3)
ax.set_xticks([]), ax.set_yticks([])
ax.axis([0, img.shape[1], img.shape[0], 0]);

Test Chan-Vese on natural images

Astronaut

Circle initialization (same as in the snake)





In [16]:



    


def small_disk(image_size, center, radius):
    centerX, centerY = center
    res = np.ones(image_size)
    res[centerY, centerX] = 0.
    return (radius-distance(res)) / (radius*3)



    











In [44]:



    


original = data.astronaut()
astro = rgb2gray(original)
image = img_as_float(astro)
init = small_disk(image.shape, (220, 100), 100)
cv = chan_vese(image, mu=0.07, lambda1=1, lambda2=1, tol=1e-3, max_iter=1,
               dt=0.5, init_level_set=init, extended_output=True)



    











In [47]:



    


plt.figure(figsize=(6,6)); plt.imshow(cv[1],cmap='gray'); plt.title('Initial level set'); plt.axis('off');



    


















    
























In [30]:



    


original = data.astronaut()
astro = rgb2gray(original)
image = img_as_float(astro)
init = small_disk(image.shape, (220, 100), 100)
cv = chan_vese(image, mu=0.07, lambda1=1, lambda2=1, tol=1e-3, max_iter=5000,
               dt=0.5, init_level_set=init, extended_output=True)



    











In [31]:



    


plot_results(image, cv, 0, (10,10))
## Create slides for animation original = data.astronaut() astro = rgb2gray(original) init = small_disk(image.shape, (220, 100), 100) image = img_as_float(astro) segs = [] for i in range(1,500): segs.append(chan_vese(image, mu=0.07, lambda1=1, lambda2=1, tol=1e-3, max_iter=i, dt=0.5, init_level_set=init)) np.save('circastro.npy', np.array(segs)) 




In [ ]:



    


segs = np.load('circastro.npy')
anim = animate(image, segs);



    











In [58]:



    


HTML(anim.to_html5_video())



    


















    
Out[58]:

Checkerboard initialization





In [36]:



    


original = data.astronaut()
astro = rgb2gray(original)
image = img_as_float(astro)
cv = chan_vese(image, mu=0.07, lambda1=1, lambda2=1, tol=1e-3, max_iter=1,
               dt=0.5, init_level_set="checkerboard", extended_output=True)



    











In [43]:



    


plt.figure(figsize=(6,6)); plt.imshow(cv[1],cmap='gray'); plt.title('Initial level set'); plt.axis('off');



    


















    
























In [32]:



    


original = data.astronaut()
astro = rgb2gray(original)
image = img_as_float(astro)
cv = chan_vese(image, mu=0.07, lambda1=1, lambda2=1, tol=1e-3, max_iter=5000,
               dt=0.5, init_level_set="checkerboard", extended_output=True)



    











In [33]:



    


plot_results(image, cv, 0, (10,10))
## Create slides for animation original = data.astronaut() astro = rgb2gray(original) image = img_as_float(astro) segs = [] for i in range(1,500): segs.append(chan_vese(image, mu=0.25, lambda1=1, lambda2=1, tol=1e-3, max_iter=i, dt=0.5, init_level_set=init)) np.save('astrosegs.npy', np.array(segs)) 




In [ ]:



    


segs = np.load('astrosegs.npy')
anim = animate(image, segs);



    











In [55]:



    


HTML(anim.to_html5_video())



    


















    
Out[55]:

Cameraman





In [11]:



    


image = img_as_float(data.camera())
cv = chan_vese(image, mu=0.25, lambda1=1, lambda2=1, tol=1e-3, max_iter=200,
               dt=0.5, init_level_set="checkerboard", extended_output=True)



    











In [12]:



    


plot_results(image, cv, 0)

Text on paper





In [22]:



    


image = img_as_float(data.text())
cv = chan_vese(image, mu=0.01, lambda1=1, lambda2=1, tol=1e-3, max_iter=5000,
               dt=0.5, init_level_set="checkerboard", extended_output=True)



    











In [23]:



    


plot_results(image, cv, 0, (18,8))

Test on medical image

MR T1-w brain image





In [27]:



    


from glob import glob
fns = sorted(glob('imgs/*.tif'))
imgs = [Image.open(f) for f in fns]



    











In [28]:



    


def show(x,ax,t=''): ax.imshow(np.rot90(x,3),cmap='gray');ax.axis('off');ax.set_title(t);



    











In [29]:



    


f, ax1 = plt.subplots(1,1,figsize=(12,8))
show(imgs[0],ax1,'Img 1');



    


















    
























In [30]:



    


img1, img2 = [np.asarray(im) for im in imgs]



    











In [31]:



    


cv = chan_vese(img1, mu=0.01, lambda1=1, lambda2=1, tol=1e-3, max_iter=5000,
               dt=0.5, init_level_set="checkerboard", extended_output=True)



    











In [32]:



    


plot_results(img1, cv, 3, (8,10))

Content source: jcreinhold/jcreinhold.github.io

Similar notebooks:

notebook.community | gallery | about