主动轮廓模型:Snake模型的python实现

本文介绍使用Skimage和Numpy库实现Snake模型进行图像分割的方法。Snake模型是一种活动轮廓模型,适用于2D单通道或多通道图像的分割。通过调整长度、平滑度、亮度吸引和边缘吸引等参数,Snake模型可以捕捉图像特征并拟合轮廓。文中提供了详细的函数参数说明、代码示例及结果展示。

质量声明:原创文章,内容质量问题请评论吐槽。如对您产生干扰,可私信删除。
主要参考:Active Contour Model — skimage v0.16.dev0 docs - scikit-image


skimage实现

函数声明

Active contours by fitting snakes to features of images. Supports single and multichannel 2D images. Snakes can be periodic (for segmentation) or have fixed and/or free ends. The output snake has the same length as the input boundary. As the number of points is constant, make sure that the initial snake has enough points to capture the details of the final contour.

active_contour(image, snake, alpha=0.01, beta=0.1, w_line=0, w_edge=1, gamma=0.01,
        bc='periodic', max_px_move=1.0, max_iterations=2500, convergence=0.1)

    Parameters
    ----------
    image : (N, M) or (N, M, 3) ndarray
        Input image.
    snake : (N, 2) ndarray
        Initial snake coordinates. For periodic boundary conditions, endpoints
        must not be duplicated.
    alpha : float, optional
        Snake length shape parameter. Higher values makes snake contract
        faster.
    beta : float, optional
        Snake smoothness shape parameter. Higher values makes snake smoother.
    w_line : float, optional
        Controls attraction to brightness. Use negative values to attract toward
        dark regions.
    w_edge : float, optional
        Controls attraction to edges. Use negative values to repel snake from
        edges.
    gamma : float, optional
        Explicit time stepping parameter.
    bc : {'periodic', 'free', 'fixed'}, optional
        Boundary conditions for worm. 'periodic' attaches the two ends of the
        snake, 'fixed' holds the end-points in place, and 'free' allows free
        movement of the ends. 'fixed' and 'free' can be combined by parsing
        'fixed-free', 'free-fixed'. Parsing 'fixed-fixed' or 'free-free'
        yields same behaviour as 'fixed' and 'free', respectively.
    max_px_move : float, optional
        Maximum pixel distance to move per iteration.
    max_iterations : int, optional
        Maximum iterations to optimize snake shape.
    convergence: float, optional
        Convergence criteria.

    Returns
    -------
    snake : (N, 2) ndarray
        Optimised snake, same shape as input parameter.

    References
    ----------
    .. [1]  Kass, M.; Witkin, A.; Terzopoulos, D. "Snakes: Active contour
            models". International Journal of Computer Vision 1 (4): 321
            (1988). DOI:`10.1007/BF00133570`

代码示例

import numpy as np
from matplotlib import pyplot as plt
from skimage.color import rgb2gray
from skimage import data
from skimage.filters import gaussian
from skimage.segmentation import active_contour

img = data.astronaut() # 读入图像
img = rgb2gray(img) # 灰度化

# 圆的参数方程:(220, 100) r=100
t = np.linspace(0, 2*np.pi, 400) # 参数t, [0,2π]
x = 220 + 100*np.cos(t)
y = 100 + 100*np.sin(t)

# 构造初始Snake
init = np.array([x, y]).T # shape=(400, 2)

# Snake模型迭代输出
snake = active_contour(gaussian(img,3), snake=init, alpha=0.1, beta=1, gamma=0.01, w_line=0, w_edge=10)

# 绘图显示
plt.figure(figsize=(5, 5))
plt.imshow(img, cmap="gray")
plt.plot(init[:, 0], init[:, 1], '--r', lw=3)
plt.plot(snake[:, 
<think>我们正在处理一个关于“自动结合局部与全局信息活动轮廓模型”的Python实现问题。活动轮廓模型(Active Contour Model)是一种图像分割技术,它结合了局部和全局信息以改进分割效果。常见的模型有经典的Snake模型,以及后来改进的模型如CV模型(Chan-Vese)和局部区域模型。这里我们可能需要实现一个结合了局部和全局信息的模型。 根据问题,我们需要实现一个能够自动结合局部和全局信息的活动轮廓模型。在图像分割领域,Chan-Vese模型是一种全局模型,而局部二元拟合(LBF)模型则是一种局部模型。因此,我们可能会参考一些结合两者的方法,例如Local and Global Intensity Fitting (LGIF)模型。 由于实现一个完整的LGIF模型可能比较复杂,我们将提供一个简化版本的实现思路,并给关键代码步骤。 步骤: 1. 初始化轮廓(通常是一个矩形或圆形)。 2. 定义能量函数,该函数同时包含全局项(如CV模型中的全局区域拟合)和局部项(如LBF模型中的局部区域拟合)。 3. 使用水平集方法(Level Set)来表示轮廓。 4. 通过梯度下降法最小化能量函数,从而演化轮廓。 我们将实现一个结合了全局和局部信息的活动轮廓模型,这里我们参考LGIF模型的能量函数: 能量函数一般形式: $$E = \alpha E_{global} + \beta E_{local} + \mu \cdot \text{正则项}$$ 其中: - $E_{global}$ 是全局能量项(如Chan-Vese的能量项) - $E_{local}$ 是局部能量项(如LBF的能量项) - 正则项用于保持轮廓的光滑性 具体实现中,我们将使用水平集函数$\phi$,并求解偏微分方程(PDE)来更新水平集函数。 由于完整的实现非常复杂,我们将提供一个简化版本,重点展示框架和关键步骤。 注意:我们将使用skimage库中的一些图像处理函数,以及numpy进行数值计算。 实现步骤: 1. 初始化水平集函数(例如,一个矩形区域内部为负,外部为正)。 2. 定义全局能量项(Chan-Vese)和局部能量项(LBF)需要的变量。 3. 迭代更新水平集函数,直到收敛或达到最大迭代次数。 这里,我们参考Chan-Vese模型和LBF模型,结合两者的能量项。 由于代码较长,我们将分步骤编写。 注意:由于这是一个复杂的模型,我们只提供一个基础框架,实际应用可能需要调整参数和细节。 下面我们开始编写代码: </think> ### 自动结合局部与全局信息活动轮廓模型Python实现 活动轮廓模型(Active Contour Model)是一种用于图像分割的经典方法,其中**结合局部与全局信息**的模型能够更有效地处理灰度不均匀的图像。以下是一个基于改进的Chan-Vese模型(全局)和局部二值拟合(LBF)的Python实现框架: ```python import numpy as np import cv2 import matplotlib.pyplot as plt from scipy.ndimage import gaussian_filter class ActiveContourModel: def __init__(self, image, init_contour, lambda1=1, lambda2=1, mu=0.2, nu=0, epsilon=1.0, dt=0.1, max_iter=200): """ 参数说明: image: 输入图像 (灰度图) init_contour: 初始轮廓 (Nx2数组) lambda1, lambda2: 全局能量项权重 mu: 长度正则化权重 nu: 面积正则化权重 epsilon: Heaviside函数平滑参数 dt: 时间步长 max_iter: 最大迭代次数 """ self.image = image.astype(np.float32) self.contour = init_contour.copy() self.lambda1 = lambda1 self.lambda2 = lambda2 self.mu = mu self.nu = nu self.epsilon = epsilon self.dt = dt self.max_iter = max_iter self.height, self.width = image.shape # 创建符号距离函数 self.phi = self._create_sdf(init_contour) def _create_sdf(self, contour): """创建初始符号距离函数(水平集函数)""" phi = np.ones((self.height, self.width)) * -1 cv2.drawContours(phi, [contour], 0, 1, thickness=cv2.FILLED) dist = cv2.distanceTransform( (phi < 0).astype(np.uint8), cv2.DIST_L2, 3 ) return np.where(phi > 0, dist, -dist) def _heaviside(self, x): """平滑的Heaviside函数""" return 0.5 * (1 + (2/np.pi) * np.arctan(x/self.epsilon)) def _compute_local_energy(self): """计算局部二值拟合能量项""" # 使用高斯核计算局部均值 K = self._gaussian_kernel(9, 3.0) f1 = cv2.filter2D(self.image, -1, K) f2 = cv2.filter2D(1 - self.image, -1, K) return f1, f2 def _gaussian_kernel(self, size, sigma): """生成高斯核""" ax = np.linspace(-(size-1)/2, (size-1)/2, size) x, y = np.meshgrid(ax, ax) kernel = np.exp(-(x**2 + y**2)/(2*sigma**2)) return kernel / kernel.sum() def evolve(self): """执行轮廓演化""" for _ in range(self.max_iter): # 计算全局能量项(Chan-Vese) H = self._heaviside(self.phi) c1 = np.sum(self.image * H) / (np.sum(H) + 1e-8) c2 = np.sum(self.image * (1 - H)) / (np.sum(1 - H) + 1e-8) global_force = self.lambda1 * (self.image - c1)**2 - self.lambda2 * (self.image - c2)**2 # 计算局部能量项(LBF) f1, f2 = self._compute_local_energy() local_force = self.lambda1 * (self.image - f1)**2 - self.lambda2 * (self.image - f2)**2 # 组合全局和局部能量 total_force = 0.7*global_force + 0.3*local_force # 计算曲率正则化 dy, dx = np.gradient(self.phi) dxx = np.gradient(dx, axis=1) dyy = np.gradient(dy, axis=0) curvature = (dx**2 * dyy - 2*dx*dy*dxx + dy**2 * dxx) / (dx**2 + dy**2 + 1e-8) # 更新水平集函数 dphi_dt = total_force * self._dirac(self.phi) + self.mu * curvature self.phi += self.dt * dphi_dt # 重新初始化水平集(保持SDF属性) self.phi = self._reinitialize(self.phi) # 提取新轮廓 self.contour = self._extract_contour() return self.contour def _dirac(self, x): """Dirac delta函数""" return self.epsilon / (np.pi * (self.epsilon**2 + x**2)) def _reinitialize(self, phi): """重新初始化水平集函数为符号距离函数""" return cv2.distanceTransform( (phi > 0).astype(np.uint8), cv2.DIST_L2, 3 ) - cv2.distanceTransform( (phi <= 0).astype(np.uint8), cv2.DIST_L2, 3 ) def _extract_contour(self): """从水平集函数提取轮廓""" contours, _ = cv2.findContours( (self.phi > 0).astype(np.uint8), cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE ) return contours[0] if contours else self.contour # 使用示例 if __name__ == "__main__": # 读取图像并预处理 image = cv2.imread('input.png', cv2.IMREAD_GRAYSCALE) image = cv2.GaussianBlur(image, (5,5), 0) # 创建初始轮廓(中心矩形) h, w = image.shape init_contour = np.array([[[50,50]], [[50,h-50]], [[w-50,h-50]], [[w-50,50]]]) # 初始化并演化轮廓 acm = ActiveContourModel(image, init_contour, lambda1=1, lambda2=2, mu=0.04, max_iter=300) final_contour = acm.evolve() # 可视化结果 result = cv2.cvtColor(image, cv2.COLOR_GRAY2BGR) cv2.drawContours(result, [final_contour], -1, (0,255,0), 2) plt.imshow(result), plt.title('分割结果'), plt.show() ``` ### 关键实现要点 1. **水平集方法**:使用水平集函数$\phi$表示轮廓,零水平集对应实际边界 2. **能量函数组合**: - **全局项**:Chan-Vese模型的区域能量$E_{global} = \lambda_1 \int_{\Omega} (I-c_1)^2 dxdy + \lambda_2 \int_{\Omega} (I-c_2)^2 dxdy$ - **局部项**:局部二值拟合能量$E_{local} = \lambda_1 \int_{\Omega} K_\sigma * (I-f_1)^2 dxdy + \lambda_2 \int_{\Omega} K_\sigma * (I-f_2)^2 dxdy$ 3. **正则化项**:曲率项保持轮廓平滑$\mu \cdot \text{div}(\frac{\nabla \phi}{|\nabla \phi|})$ 4. **自适应权重**:通过调整全局和局部项的权重系数平衡两者贡献 ### 参数调优建议 1. **λ1/λ2**:控制前景/背景拟合权重(通常λ2 > λ1) 2. **μ**:长度约束权重(0.01-0.5) 3. **ε**:Heaviside函数平滑参数(1.0-2.0) 4. **dt**:时间步长(0.05-0.5)需满足CFL条件 实际应用中,该方法在医学图像(如MRI、CT)分割中效果显著,特别是对于边界模糊或灰度不均匀的目标[^1]。
评论 13
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值