简介这份资源提供Chan-Vese图像分割算法的Python实现代码面向正在学习图像处理、计算机视觉或活动轮廓模型的学生与开发者。Chan-Vese模型是首个完全脱离图像梯度信息、依靠区域信息构建能量函数的活动轮廓模型为后续分割方法开辟了新方向适合希望深入理解其数学原理并动手复现的读者。压缩包为7z格式共2个文件包含1个py脚本与1个bmp示例图像整体仅6KB轻量易取脚本内附有注释便于理解能量函数迭代、曲线演化与水平集更新的实现逻辑也方便按需修改参数或迁移到自己的项目中。目前已有1301人学习下载说明该实现具备一定的参考价值。读者可借助示例图直接运行验证分割效果并对照代码梳理CV模型的区域拟合思路为图像分割实验与算法改进提供可复用的基础代码。1. 从灰度不均的医学图像说起Chan-Vese 算法到底解决什么问题做过图像分割的人大概都遇到过这种场景一张工业检测的零件图背景亮度从左上到右下缓慢变化目标区域内部灰度也不均匀用阈值法切出来的结果要么缺一块要么多一块边缘还毛糙。这时候有人会推荐你试试 Chan-Vese 算法。它属于活动轮廓模型里的区域型方法核心思路不是靠梯度找边界而是靠“区域内部灰度尽量一致、区域外部灰度尽量一致”这个假设去驱动一条曲线演化最终停在目标和背景的分界上。和传统基于边缘的 Snake 模型相比Chan-Vese 对弱边缘、模糊边界、噪声的容忍度明显更高代价是计算量偏大、对初始化和参数敏感。这篇笔记面向想用 Python 复现 Chan-Vese 的从业者从数学骨架讲到可运行代码再到参数怎么调、坑在哪尽量让你照着敲一遍就能跑出结果。2. Chan-Vese 的数学骨架与水平集实现路径2.1 能量泛函到底在最小化什么Chan-Vese 模型把图像分成两个区域曲线内部和曲线外部。它假设每个区域的灰度近似为常数于是定义能量F(c1, c2, C) μ·Length(C) ν·Area(inside(C)) λ1·∫|u0 - c1|² dx λ2·∫|u0 - c2|² dx其中 u0 是输入图像C 是演化曲线c1 和 c2 分别是曲线内、外的平均灰度。前两项控制曲线的长度和面积防止曲线无限膨胀或收缩后两项是拟合项衡量用两个常数去近似整幅图像的误差。最小化这个能量就是让曲线在“平滑”和“拟合”之间找平衡。直接对曲线 C 求导很麻烦所以标准做法是引入水平集函数 φ(x,y)用 φ 的零水平集表示曲线。φ 0 是内部φ 0 是外部。再引入 Heaviside 函数 H(φ) 和 Dirac 函数 δ(φ)把能量改写成关于 φ 的泛函。对 φ 求变分得到演化方程∂φ/∂t δ(φ) [ μ·div(∇φ/|∇φ|) - ν - λ1(u0 - c1)² λ2(u0 - c2)² ]这个方程就是代码里要迭代更新的核心。c1 和 c2 每轮根据当前 φ 重新计算c1 ∫ u0·H(φ) dx / ∫ H(φ) dx c2 ∫ u0·(1 - H(φ)) dx / ∫ (1 - H(φ)) dx实际实现时H 和 δ 要用正则化版本否则数值不稳定。常见做法是用 Hε(z) 0.5(1 2/π·arctan(z/ε))δ 就是它的导数。ε 一般取 1 左右。2.2 为什么用水平集而不是直接参数化曲线参数化曲线比如用一组控制点表示 Snake在演化过程中容易发生自交要额外做重采样拓扑变化一条曲线分裂成两条几乎没法处理。水平集把曲线隐式地表示为高维函数的零水平集拓扑变化天然支持曲线分裂、合并都不用特殊处理。代价是计算维度升高整幅图像每个像素都要更新 φ内存和耗时都上去了。对于 512×512 的图单次迭代在纯 NumPy 下大概几十毫秒几百次迭代就是十几秒这是后面要重点优化的地方。2.3 用 NumPy 写一版可读的 Chan-Vese下面这版代码追求可读性注释写清楚每一步对应公式的哪一项。输入图像先转灰度、归一化到 [0,1]然后初始化 φ 为内部正、外部负的符号距离函数。import numpy as np from scipy.ndimage import gaussian_filter def chan_vese(image, num_iter300, mu0.2, nu0.0, lambda11.0, lambda21.0, epsilon1.0, dt0.5): image: 2D float array, 归一化到 [0,1] num_iter: 迭代次数 mu: 长度项权重越大曲线越平滑 nu: 面积项权重正数收缩曲线负数扩张 lambda1/lambda2: 内外拟合项权重 epsilon: 正则化参数控制 H 和 delta 的平滑程度 dt: 时间步长太大不稳定太小收敛慢 # 转灰度并归一化 if image.ndim 3: image image.mean(axis2) img image.astype(np.float64) img (img - img.min()) / (img.max() - img.min() 1e-8) # 初始化 phi中心圆内为正圆外为负 h, w img.shape y, x np.mgrid[0:h, 0:w] cx, cy, r w / 2.0, h / 2.0, min(h, w) / 4.0 phi np.sqrt((x - cx)**2 (y - cy)**2) - r phi phi / (np.abs(phi).max() 1e-8) # 归一化数值更稳 for i in range(num_iter): # 正则化 Heaviside 和 Dirac H 0.5 * (1 (2 / np.pi) * np.arctan(phi / epsilon)) delta (1 / np.pi) * (epsilon / (epsilon**2 phi**2)) # 计算曲线内外的平均灰度 c1, c2 inside H.sum() outside (1 - H).sum() c1 (img * H).sum() / (inside 1e-8) c2 (img * (1 - H)).sum() / (outside 1e-8) # 计算曲率项 div(grad(phi)/|grad(phi)|) phi_y, phi_x np.gradient(phi) grad_norm np.sqrt(phi_x**2 phi_y**2) 1e-8 nx phi_x / grad_norm ny phi_y / grad_norm # 对法向量求散度 n_y, _ np.gradient(ny) _, n_x np.gradient(nx) curvature n_y n_x # 演化方程各项 fit -lambda1 * (img - c1)**2 lambda2 * (img - c2)**2 phi phi dt * delta * (mu * curvature - nu fit) # 每 50 轮做一次轻微平滑抑制数值噪声 if i % 50 0: phi gaussian_filter(phi, sigma0.5) return phi 0 # 返回二值掩膜这段代码里几个关键点值得展开。H和delta用的是 arctan 正则化比直接用阶跃函数稳定得多epsilon控制过渡带宽度取 1.0 是常见默认值。c1和c2每轮重算这是区域型模型的核心也是它比边缘型模型抗噪的原因。曲率项用np.gradient两次求导近似虽然精度一般但对 300 轮迭代来说够用。dt取 0.5 是经验值超过 1.0 很容易发散表现为phi出现 NaN 或者整幅图全白全黑。每 50 轮做一次高斯平滑是血泪经验不加的话迭代到 200 轮以后会出现棋盘状伪影。2.4 跑通最小示例与结果检查拿一张简单的合成图验证import numpy as np from PIL import Image # 生成一张背景渐变、目标内部灰度不均的测试图 h, w 256, 256 y, x np.mgrid[0:h, 0:w] bg 0.3 0.4 * (x / w) # 背景从左到右变亮 target ((x - 128)**2 (y - 128)**2) 60**2 img bg.copy() img[target] 0.6 0.2 * np.sin(x[target] / 20) # 目标内部有波动 img np.random.normal(0, 0.02, img.shape) # 加噪声 mask chan_vese(img, num_iter300, mu0.2, nu0.0) print(前景像素占比:, mask.mean()) Image.fromarray((mask * 255).astype(np.uint8)).save(cv_result.png)跑完看cv_result.png如果前景占比在 0.15 到 0.25 之间说明大致收敛到目标区域了。如果占比接近 0 或 1多半是nu或dt设得不对。检查中间过程可以每 50 轮存一张图观察曲线是从初始圆逐渐贴到目标边界还是一直在震荡。3. 参数怎么调mu、nu、lambda 与迭代次数的实际影响3.1 mu 和 nu 的配合关系mu是长度项权重直接决定曲线平滑程度。mu太小曲线会钻进目标内部的灰度波动里把一块区域切成好几块mu太大曲线会缩成一个圆忽略真实边界。经验上对于 256×256 的图mu在 0.1 到 0.5 之间比较稳图像越大mu可以适当调大。nu是面积项权重控制曲线整体膨胀还是收缩。nu为正时曲线倾向收缩为负时倾向扩张。很多开源实现默认nu0让长度项和拟合项自己平衡。但如果你的初始圆明显小于目标可以设nu-0.01到-0.05帮助扩张初始圆大于目标就设正值。注意nu的量级要比mu小一个数量级否则会压过曲率项导致曲线失控。3.2 lambda1 和 lambda2 什么时候要改默认lambda1 lambda2 1.0表示内外拟合项同等重要。如果目标区域比背景小很多可以适当增大lambda1让模型更关注内部拟合反之增大lambda2。但这两个参数改动超过 5 倍就容易出问题比如lambda110时曲线会疯狂扩张去覆盖更多像素以降低内部方差。我一般只在目标占比低于 5% 时才动它们且每次只调 0.5。3.3 迭代次数与收敛判断固定迭代次数是最省事的做法但不同图像收敛速度差很多。更稳的方式是监控phi的变化量phi_old phi.copy() # ... 更新 phi ... change np.abs(phi - phi_old).mean() if change 1e-4: print(f第 {i} 轮收敛) break1e-4这个阈值对归一化后的phi比较合适。如果跑满 300 轮还没收敛要么是dt太小要么是参数组合让曲线在两组解之间来回跳。这时候先把dt提到 0.8 试试还不行就检查mu是不是太小。3.4 不同参数组合的对比参数组合munudt典型表现默认0.20.00.5适合大多数 256 图收敛约 200 轮强平滑0.50.00.5边界圆滑细小凹陷被忽略弱平滑0.050.00.3边界贴合紧但易受噪声影响扩张型0.2-0.020.5初始圆偏小时用注意别过冲收缩型0.20.020.5初始圆偏大时用可能缩过头这张表是多次试出来的经验区间不是理论最优。实际图像还得自己微调尤其是mu和dt的乘积它俩共同决定每次迭代曲线移动的幅度。4. 避坑与排查Chan-Vese 落地时最容易翻车的五个地方4.1 现象迭代几轮后 phi 全变 NaN原因dt太大或者epsilon太小导致delta在phi接近 0 时数值爆炸。delta epsilon / (pi * (epsilon^2 phi^2))当epsilon0.1且phi0时delta能到 3 以上乘上dt1.0直接让phi跳出发散。解决epsilon不要低于 0.5dt不要超过 0.8。如果必须用大dt加速先把epsilon提到 2.0过渡带变宽但稳定。另外在更新后加一句phi np.clip(phi, -10, 10)兜底。4.2 现象曲线停在初始位置几乎不动原因mu太大或者delta在初始phi处接近 0。初始phi是符号距离函数归一化后很多像素的phi绝对值大于 3delta已经衰减到很小演化速度极慢。解决初始化后不要做过度归一化或者把初始phi乘一个系数让它集中在 [-2, 2]。更直接的办法是前 50 轮把dt临时调到 1.0等曲线开始移动后再降回 0.5。4.3 现象分割结果把目标切碎成多块原因mu太小长度项约束不足曲线钻进目标内部的灰度波动里。区域型模型假设每个区域灰度近似常数目标内部如果有明显纹理模型会认为那是多个区域。解决先把图像做一次高斯平滑sigma1.0到2.0抹掉高频纹理再送入算法。如果还不能解决把mu提高 2 到 3 倍。实在不行就换模型Chan-Vese 本身不适合纹理丰富的目标。4.4 现象边界明显偏移曲线停在目标内部或外部原因nu设得不对或者lambda1和lambda2失衡。nu为正时曲线有收缩倾向如果目标本身比初始圆大曲线就停在目标内部。解决先设nu0跑一遍看曲线往哪个方向偏。往内偏就设负nu往外偏就设正nu每次调 0.005。同时检查c1和c2的值如果两者差距很小比如都接近 0.5说明图像对比度太低模型本身分辨不了得先做对比度增强。4.5 现象跑一次要几十秒调参效率极低原因纯 NumPy 逐像素更新Python 循环 300 次每次都在整幅图上算np.gradient和求和。解决三个方向。第一把图像降采样到 128×128 调参参数确定后再上原图跑。第二用scipy.ndimage的uniform_filter替代部分求和操作。第三如果长期用这个算法把核心迭代用 Cython 或 Numba 重写njit装饰器加在迭代函数上通常能快 20 到 50 倍。下面是一个 Numba 加速的骨架from numba import njit njit(cacheTrue) def evolve(phi, img, mu, nu, lambda1, lambda2, epsilon, dt): h, w phi.shape for _ in range(1): # 外层循环在 Python 里控制 # 这里写逐像素更新避免创建大临时数组 # 具体实现略思路是用双重 for 循环替代向量化 pass return phiNumba 版本要注意不能用np.gradient得手写差分。第一次编译有开销之后每次调用就快了。5. 进阶技巧多相 Chan-Vese 与初始化策略单相 Chan-Vese 只能分前景背景两类遇到多目标场景就力不从心。多相版本用多个水平集函数n 个函数可以表示 2^n 个区域。比如两个水平集 φ1 和 φ2四个区域分别是 {φ10, φ20}、{φ10, φ20}、{φ10, φ20}、{φ10, φ20}每个区域有自己的平均灰度 c11、c10、c01、c00。演化方程变成两个耦合的 PDE实现时把单相的c1、c2换成四个均值拟合项按区域归属分别计算。# 两相版本的核心改动 H1 0.5 * (1 (2/np.pi) * np.arctan(phi1/epsilon)) H2 0.5 * (1 (2/np.pi) * np.arctan(phi2/epsilon)) # 四个区域指示函数 r11 H1 * H2 r10 H1 * (1 - H2) r01 (1 - H1) * H2 r00 (1 - H1) * (1 - H2) c11 (img * r11).sum() / (r11.sum() 1e-8) c10 (img * r10).sum() / (r10.sum() 1e-8) c01 (img * r01).sum() / (r01.sum() 1e-8) c00 (img * r00).sum() / (r00.sum() 1e-8) # phi1 的拟合项按 phi2 的正负分两组 fit1 -lambda1 * ((img - c11)**2 * H2 (img - c10)**2 * (1 - H2)) \ lambda2 * ((img - c01)**2 * H2 (img - c00)**2 * (1 - H2))多相版本对初始化更敏感两个初始圆如果完全重合r11和r00会有一个接近零均值计算不稳定。常见做法是让两个初始圆错开一个偏左上、一个偏右下或者一个圆大一个圆小。另外多相版本的迭代次数通常要翻倍因为两个水平集互相影响收敛更慢。初始化策略本身也值得单独说。默认用中心圆是最省事的但如果目标不在图像中心曲线要花很多轮才能移过去还可能卡在中间。更好的做法是用 Otsu 阈值先粗分割把前景区域的质心和大致半径算出来用这个圆做初始 φ。代码就三行from skimage.filters import threshold_otsu thresh threshold_otsu(img) binary img thresh ys, xs np.where(binary) cx, cy xs.mean(), ys.mean() r np.sqrt(binary.sum() / np.pi) * 0.8 # 略小于等效圆留扩张余地这样初始曲线已经贴近目标迭代 100 轮左右就能收敛比盲猜中心圆快一倍以上。如果目标有多个就取最大连通域的质心或者干脆用多相版本。最后一个习惯每次调参只改一个变量改完存一张结果图文件名带上参数值。Chan-Vese 的参数空间不算大但组合起来也有几十种不记录的话跑两天就忘了哪组最好。我一般用fcv_mu{mu}_nu{nu}_dt{dt}.png这种命名回头对比一目了然。希望帮到你。本文还有配套的精品资源点击获取