Level Set图像分割方法及python实例

Minerva

修改于 2020-05-23 10:30:35

4.7K2

修改于 2020-05-23 10:30:35

文章被收录于专栏：Python编程和深度学习Python编程和深度学习

Level Set 模型介绍

Level Set是基于能量的图像分割方法，通过求解最小能量泛函，得到目标轮廓的表达式。其轮廓表示借鉴了水平集的概念，将低维度的曲线嵌入了高维度的曲面中。

一。能量函数的构造

以下图为例，黑色图形为带分割图形，白色圆为轮廓曲线C，常用的能量方程如下所示，

其中c1,c2,分别为圆圈内/外像素点的平均值。u0为图像的像素值。整个能量项的意义即为：圆圈外像素减去他们均值的平方和加上圆圈内像素减去他们均值的平方和，还有一些规则化项，如曲线的长度，曲线围成的面积等项。u,v 大于等于0，lamda1,lamda2 大于0. 一般规定v=0; lamda1=lamda2=1.一般来说，最小化上述能量方程即可求得轮廓C的方程。

二 . 水平集思想

水平集LevelSet需要给定初始的轮廓曲线，然后根据泛函能量最小化，进行曲线演化，由初始曲线演化到目标轮廓曲线。水平集的方法，用的是一种隐式函数的方法，通过高维空间来表示低位的空间。其更新的不是曲线离散点的坐标（snake算法曲线演化），而是更新整张图片像素点到曲线的有向距离场

下面用一个实例展示LevelSet方法，源代码：https://github.com/Ramesh-X/Level-Set 通过运行python -m lv_set.Main 即可实现LevelSet目标分割。下面贴一下主函数：

from scipy.misc import imread
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import scipy.ndimage.filters as filters
from skimage import measure
import lv_set.drlse_algo as drlse
import numpy as np
import cv2

img = np.array(imread('gourd.bmp', True), dtype='float32')
# im_t = img[:, :, 1]
# plt.imshow(img)
# plt.show()
# parameters
timestep = 1        # time step
mu = 0.2/timestep   # coefficient of the distance regularization term R(phi)
iter_inner = 4
iter_outer = 25
lmda = 2            # coefficient of the weighted length term L(phi)
alfa = -9           # coefficient of the weighted area term A(phi)
epsilon = 2.0       # parameter that specifies the width of the DiracDelta function

sigma = 0.8         # scale parameter in Gaussian kernel
img_smooth = filters.gaussian_filter(img, sigma)    # smooth image by Gaussian convolution
[Iy, Ix] = np.gradient(img_smooth)
f = np.square(Ix) + np.square(Iy)
g = 1 / (1+f)    # edge indicator function.

# initialize LSF as binary step function
c0 = 2
initialLSF = c0 * np.ones(img.shape)
cv2.imwrite('initialLSF.png', initialLSF)
# generate the initial region R0 as two rectangles
# initialLSF[24:35, 19:25] = -c0
print('initialLSF.shape', initialLSF.shape)
initialLSF[24:35, 20:26] = -c0
phi = initialLSF.copy()
cv2.imwrite('phi.png', phi)
plt.ion()
fig1 = plt.figure(1)
# plt.imshow(phi)
# plt.show()
# plt.pause(30)
def show_fig1():
    ax1 = fig1.add_subplot(111, projection='3d')
    y, x = phi.shape
    x = np.arange(0, x, 1)
    y = np.arange(0, y, 1)
    X, Y = np.meshgrid(x, y)
    ax1.plot_surface(X, Y, -phi, rstride=2, cstride=2, color='r', linewidth=0, alpha=0.6, antialiased=True)
    ax1.contour(X, Y, phi, 0, colors='g', linewidths=2)

show_fig1()
fig2 = plt.figure(2)


def show_fig2():
    contours = measure.find_contours(phi, 0)
    ax2 = fig2.add_subplot(111)
    ax2.imshow(img, interpolation='nearest', cmap=plt.cm.gray)
    for n, contour in enumerate(contours):
        ax2.plot(contour[:, 1], contour[:, 0], linewidth=2)


show_fig2()
print('show fig 2 first time')

potential = 2
if potential == 1:
    potentialFunction = 'single-well'  # use single well potential p1(s)=0.5*(s-1)^2, which is good for region-based model
elif potential == 2:
    potentialFunction = 'double-well'  # use double-well potential in Eq. (16), which is good for both edge and region based models
else:
    potentialFunction = 'double-well'  # default choice of potential function

# start level set evolution
for n in range(iter_outer):
    phi = drlse.drlse_edge(phi, g, lmda, mu, alfa, epsilon, timestep, iter_inner, potentialFunction)
    if np.mod(n, 2) == 0:
        print('show fig 2 for %i time' % n)
        fig2.clf()
        show_fig2()
        fig1.clf()
        show_fig1()
        plt.pause(10)

# refine the zero level contour by further level set evolution with alfa=0
alfa = 0
iter_refine = 10
phi = drlse.drlse_edge(phi, g, lmda, mu, alfa, epsilon, timestep, iter_refine, potentialFunction)

finalLSF = phi.copy()
print('show final fig 2')
fig2.clf()
show_fig2()
fig1.clf()
show_fig1()
plt.pause(10)
plt.show()

原始图像为