简介

图像分割的核心是从复杂背景中分离出目标区域，但传统边缘-based方法往往受限于模糊边界或弱对比度。2001年，Tony F. Chan和Luminita Vese提出的ACWE模型（Active Contours Without Edges），结合了Mumford-Shah分割理论与水平集方法，无需依赖图像梯度（边缘信息），而是通过区域统计特征引导轮廓演化，特别适合复杂结构（如医学影像、弱边缘目标）的分割。

一、水平集方法：从曲线到高维函数的隐式表达

水平集（Level Set）是一种隐式轮廓跟踪方法，由Sethian和Osher于1988年提出。其核心思想是：将低维空间中的运动曲线/曲面，隐含地表示为高维函数的“零水平集”（即函数值为0的点集）。
举个直观例子：二维平面上的圆 $(x^2+y^2=1$)，可以看作三维空间中函数 $(f(x,y)=x^2+y^2$) 的“1水平集”（所有满足 $(f(x,y)=1$) 的点）。若我们让这个圆随时间演化（比如扩张或收缩），只需修改高维函数 $(f(x,y,t)$)，再从每个时刻的函数中提取“零水平集”，就能得到动态的轮廓。

水平集的优势在于：

1. 自动处理拓扑变化：无需手动调整轮廓的连接/分裂（如分割多个目标时，轮廓可自然分裂）；
2. 数学严谨性：通过偏微分方程（PDE）描述轮廓演化，便于融入正则项（如平滑、长度约束）。

二、Chan-Vese模型：无边缘的区域主动轮廓

Chan-Vese模型（简称CV模型）是水平集方法与Mumford-Shah分割的结合，其能量函数仅依赖区域的统计特征（如灰度均值），而非边缘梯度。

1. 能量函数的构成

对于图像$(u_0(x,y))$，CV模型的能量函数定义为：$$E(\varphi )=\mu L(\varphi )+uA(\varphi )+{\lambda }_1{\int }_{\Omega }(u_0-c_1{)}^{2}H(\varphi )dxdy+{\lambda }_2{\int }_{\Omega }(u_0-c_2{)}^2(1-H(\varphi ))dxdy$$
其中：

• $(L(\varphi ))$：轮廓的长度（正则项，约束轮廓平滑）；
• $(A(\varphi )$)：轮廓包围的面积（可选正则项，控制轮廓大小）；
• $(c_1$)：轮廓内部（$(\varphi >0$)）的灰度均值；
• $(c_2$)：轮廓外部（$(\varphi <0$)）的灰度均值；
• $(H(\varphi )$)：Heaviside函数（区分内外区域）；
• $(\mu ,u,{\lambda }_1,{\lambda }_2$)：权重参数（平衡各能量项）。

2. 平滑的Heaviside与Dirac函数

原始Heaviside函数是阶跃函数（不连续），Dirac函数是冲激函数（不可导），无法直接用于数值计算。因此CV模型使用平滑近似：

• 平滑Heaviside函数：$$H_{ϵ}(\varphi )=\frac{1}{2}\left(1+\frac{2}{\pi }\arctan \left(\frac{\varphi }{ϵ}\right)\right)$$
• 平滑Dirac函数（Heaviside的导数）：$${\delta }_{ϵ}(\varphi )=\frac{ϵ}{\pi ({ϵ}^2+{\varphi }^2)}$$其中 ($epsilon$) 是控制平滑程度的小参数（通常取1）。

3. 水平集的迭代更新

通过变分法最小化能量函数 $(E(\varphi )$)，可得到水平集函数 $(\varphi$) 的演化方程：$$\frac{\partial \varphi }{\partial t}={\delta }_{ϵ}(\varphi )\left(\mu \kappa -u-{\lambda }_1(u_0-c_1{)}^2+{\lambda }_2(u_0-c_2{)}^2\right)$$其中 $(\kappa$) 是轮廓的曲率（控制轮廓平滑）。

迭代过程中，需交替计算：

1. $(H_{ϵ}(\varphi )$) 和 $({\delta }_{ϵ}(\varphi )$)；
2. 内部均值 $(c_1$) 和外部均值 $(c_2$)；
3. 曲率 $(\kappa$)；
4. 更新 $(\varphi$)。

三、python实现

import cv2
import numpy as np
import math


class LevelSet:
    def __init__(self, src):
        """构造函数：处理输入图像"""
        if len(src.shape) == 3:
            self.src_ = cv2.cvtColor(src, cv2.COLOR_BGR2GRAY)
            self.image_ = src.copy()  # 保留彩色图用于可视化
        else:
            self.src_ = src.copy()
            self.image_ = cv2.cvtColor(src, cv2.COLOR_GRAY2BGR)

        self.src_ = self.src_.astype(np.float32)  # 转换为32位浮点型
        self.phi_ = np.zeros(self.src_.shape, dtype=np.float32)
        self.dirac_ = np.zeros(self.src_.shape, dtype=np.float32)
        self.heaviside_ = np.zeros(self.src_.shape, dtype=np.float32)
        self.curvature_ = np.zeros(self.src_.shape, dtype=np.float32)

        # 参数设置
        self.iter_num_ = 100  # 迭代次数
        self.lambda1_ = 1.0  # c1拟合项权重
        self.lambda2_ = 1.0  # c2拟合项权重
        self.mu_ = 0.1 * 255 * 255  # 长度项权重（μ）
        self.nu_ = 0.0  # 面积项权重（ν，默认0）
        self.timestep_ = 0.1  # 时间步长（δt）
        self.epsilon_ = 1.0  # 平滑参数（ε）

        self.c1_ = 0.0  # 内部灰度均值
        self.c2_ = 0.0  # 外部灰度均值

    def initializePhi(self, center, radius):
        """初始化水平集函数（圆）"""
        c = 2.0  # 内外区域的函数值（非零）
        rows, cols = self.src_.shape

        for i in range(rows):
            for j in range(cols):
                dist = math.sqrt((j - center[0]) ** 2 + (i - center[1]) ** 2)
                value = -dist + radius
                if abs(value) < 0.1:
                    self.phi_[i, j] = 0  # 零水平集（初始轮廓）
                elif value >= 0.1:
                    self.phi_[i, j] = c  # 内部
                else:
                    self.phi_[i, j] = -c  # 外部

    def computeGradient(self, src):
        """计算梯度"""
        if src.shape[0] < 2 or src.shape[1] < 2:
            return np.zeros_like(src), np.zeros_like(src)

        grad_x = np.zeros_like(src)
        grad_y = np.zeros_like(src)
        rows, cols = src.shape

        # 计算x方向梯度
        for i in range(rows):
            for j in range(cols):
                if j == 0:
                    grad_x[i, j] = src[i, j + 1] - src[i, j]
                elif j == cols - 1:
                    grad_x[i, j] = src[i, j] - src[i, j - 1]
                else:
                    grad_x[i, j] = (src[i, j + 1] - src[i, j - 1]) / 2.0

        # 计算y方向梯度
        for j in range(cols):
            for i in range(rows):
                if i == 0:
                    grad_y[i, j] = src[i + 1, j] - src[i, j]
                elif i == rows - 1:
                    grad_y[i, j] = src[i, j] - src[i - 1, j]
                else:
                    grad_y[i, j] = (src[i + 1, j] - src[i - 1, j]) / 2.0

        return grad_x, grad_y

    def computeDirac(self):
        """计算平滑Dirac函数"""
        k1 = self.epsilon_ / math.pi
        k2 = self.epsilon_ * self.epsilon_

        for i in range(self.src_.shape[0]):
            for j in range(self.src_.shape[1]):
                phi = self.phi_[i, j]
                self.dirac_[i, j] = k1 / (k2 + phi * phi)

    def computeHeaviside(self):
        """计算平滑Heaviside函数"""
        k = 2.0 / math.pi

        for i in range(self.src_.shape[0]):
            for j in range(self.src_.shape[1]):
                phi = self.phi_[i, j]
                self.heaviside_[i, j] = 0.5 * (1 + k * math.atan(phi / self.epsilon_))

    def computeCurvature(self):
        """计算曲率"""
        grad_x, grad_y = self.computeGradient(self.phi_)

        # 计算梯度的模长
        grad_norm = np.sqrt(grad_x ** 2 + grad_y ** 2)
        grad_norm[grad_norm == 0] = 1e-10  # 避免除零

        # 计算单位梯度
        dx_unit = grad_x / grad_norm
        dy_unit = grad_y / grad_norm

        # 计算单位梯度的二阶导数
        dxx, _ = self.computeGradient(dx_unit)
        _, dyy = self.computeGradient(dy_unit)

        self.curvature_ = dxx + dyy

    def computeC(self):
        """计算区域均值c1、c2"""
        sum_in = 0.0
        sum_out = 0.0
        count_in = 0.0
        count_out = 0.0

        for i in range(self.src_.shape[0]):
            for j in range(self.src_.shape[1]):
                h = self.heaviside_[i, j]
                pix = self.src_[i, j]
                sum_in += h * pix
                count_in += h
                sum_out += (1 - h) * pix
                count_out += (1 - h)

        self.c1_ = sum_in / (count_in + 1e-10)  # 避免除零
        self.c2_ = sum_out / (count_out + 1e-10)

    def showEvolving(self):
        """显示演化过程"""
        show_img = self.image_.copy()
        mask = self.phi_ >= 0  # 内部区域（φ≥0）

        # 膨胀+腐蚀去除噪声，提取轮廓
        kernel = np.ones((3, 3), np.uint8)
        mask = mask.astype(np.uint8)
        mask = cv2.dilate(mask, kernel, iterations=3)
        mask = cv2.erode(mask, kernel, iterations=3)

        # 查找并绘制轮廓
        contours, _ = cv2.findContours(mask, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_NONE)
        cv2.drawContours(show_img, contours, -1, (0, 255, 0), 2)

        cv2.imshow("Evolving Contour", show_img)
        cv2.waitKey(10)  # 延迟10ms，观察演化过程

    def evolving(self):
        """主迭代函数"""
        self.showEvolving()  # 显示初始轮廓

        for iter in range(self.iter_num_):
            self.computeHeaviside()  # 计算H_ε(φ)
            self.computeDirac()  # 计算δ_ε(φ)
            self.computeCurvature()  # 计算曲率κ
            self.computeC()  # 计算c1、c2

            # 更新水平集函数φ
            for i in range(self.src_.shape[0]):
                for j in range(self.src_.shape[1]):
                    phi = self.phi_[i, j]
                    dirac = self.dirac_[i, j]
                    curv = self.curvature_[i, j]
                    pix = self.src_[i, j]

                    # 演化项：长度项 + 面积项 + 拟合项
                    length_term = self.mu_ * dirac * curv
                    area_term = self.nu_ * dirac
                    fit_term = dirac * (-self.lambda1_ * (pix - self.c1_) ** 2
                                        + self.lambda2_ * (pix - self.c2_) ** 2)
                    delta_phi = self.timestep_ * (length_term + area_term + fit_term)

                    self.phi_[i, j] += delta_phi

            self.showEvolving()  # 显示当前演化结果

        cv2.waitKey(0)
        cv2.destroyAllWindows()


# 使用示例
if __name__ == "__main__":
    # 读取图像
    img = cv2.imread("image/11.jpeg")
    if img is None:
        # 如果没有图像，创建一个测试图像
        img = np.ones((300, 300, 3), dtype=np.uint8) * 128
        cv2.circle(img, (150, 150), 80, (255, 255, 255), -1)

    # 创建LevelSet对象
    ls = LevelSet(img)

    # 初始化水平集函数（圆心和半径）
    center = (img.shape[1] // 2, img.shape[0] // 2)
    radius = min(img.shape[0], img.shape[1]) // 4

    ls.initializePhi(center, radius)

    # 开始演化
    ls.evolving()

获取更多资料

我给大家整理了一套全网最全的人工智能学习资料（1.5T），包括：机器学习，深度学习，大模型，CV方向，NLP方向，kaggle大赛，实战项目、自动驾驶，AI就业等免费获取。