★ FEATURED ARTICLE 科学计算【免费下载链接】cvxpyA Python-embedded modeling language for convex optimization problems.项目地址https://gitcode.com/gh_mirrors/cv/cvxpy点击查看免费下载导读本文以 CVXPY 官方示例 tv_inpainting.rst 为蓝本系统讲解如何使用 CVXPY 对灰度图与彩色图进行总变分Total Variation, TV图像修复in-painting给定一幅部分像素损坏/缺失的图像通过最小化图像的 TV 范数并约束已知像素不变重建出完整图像。读完本文你将掌握 TV 修复的数学模型、CVXPY 中tv原子atom的底层实现原理、灰度与彩色两套完整可运行的求解代码以及面对超大规模锥规划问题时如何选择求解器并读懂 SCS 的迭代日志。一、问题背景什么是图像修复Inpainting图像修复指在已知部分像素值的前提下推测并填补缺失像素从而恢复整幅图像。典型应用场景包括文字/Logo 遮挡恢复如示例中被白色文字“This is the Loki test image...”覆盖的图片随机噪声/划痕去除如示例中约 70% 像素被随机置零的彩色图老旧照片修复、视频去字幕、医学影像补全等。在数学上图像修复是一个不适定ill-posed问题缺失像素有无数种补法。为此需要引入先验正则——总变分它偏好“分段平滑”piecewise smooth的图像即在像素值变化剧烈处允许跳变保留边缘在平坦区域则抑制噪声。这正是 TV 修复在图像去噪、去遮挡中效果出色的原因。二、灰度图 TV 修复数学模型2.1 图像与已知像素的表示一幅灰度图表示为 $m \times n$ 的强度矩阵 $U^{\text{orig}}$像素值通常落在 $[0, 255]$。已知像素下标集合为 $\mathcal{K} \subset {1,\ldots,m} \times {1,\ldots,n}$已知值记为 $U^{\text{orig}}_{ij}, (i,j)\in\mathcal{K}$。我们要重建的矩阵为 $U \in \mathbf{R}^{m\times n}$且必须满足保真约束$$U_{ij} U^{\text{orig}}_{ij}, \quad \forall (i,j) \in \mathcal{K}$$2.2 $\ell_2$ 总变分定义CVXPY 示例采用$\ell_2$ 总变分对矩阵按离散梯度取 2-范数再求和$$\mathop{\bf tv}(U) \sum_{i1}^{m-1}\sum_{j1}^{n-1}\left|\begin{bmatrix} U_{i1,j}-U_{ij} \ U_{i,j1}-U_{ij}\end{bmatrix}\right|_2$$即对每个像素计算“右邻像素差、下邻像素差”组成的梯度向量取其欧几里得范数再对全图求和。注意范数不加平方——这是 TV 与 Tikhonov 正则$|\nabla U|_2^2$的本质区别不加平方的 $\ell_2$ 范数使目标函数对边缘处的梯度惩罚是线性的从而允许图像中存在锐利边缘而不会像平方惩罚那样过度模糊。2.3 优化问题重建图像 $U$ 通过求解如下凸优化问题得到$$ \begin{aligned} \underset{U}{\text{minimize}} \mathop{\bf tv}(U) \ \text{subject to} U_{ij} U^{\text{orig}}_{ij}, \quad (i,j)\in\mathcal{K} \end{aligned} $$目标函数 $\mathop{\bf tv}(U)$ 是凸的范数求和约束是线性等式约束因此整个问题是一个凸优化问题可用 CVXPY 直接建模并交给锥求解器conic solver求解。三、灰度图修复CVXPY 实现与代码逐行解析3.1 加载图像与构造 Known 矩阵import matplotlib.pyplot as plt import numpy as np # 加载原始图与损坏图。 u_orig plt.imread(data/loki512.png) u_corr plt.imread(data/loki512_corrupted.png) rows, cols u_orig.shape # known 为 1 表示像素已知0 表示像素被损坏。 known np.zeros((rows, cols)) for i in range(rows): for j in range(cols): if u_orig[i, j] u_corr[i, j]: known[i, j] 1 %matplotlib inline fig, ax plt.subplots(1, 2, figsize(10, 5)) ax[0].imshow(u_orig, cmapgray) ax[0].set_title(Original Image) ax[0].axis(off) ax[1].imshow(u_corr, cmapgray); ax[1].set_title(Corrupted Image) ax[1].axis(off);关键点u_orig与u_corr是形状为(rows, cols)的 NumPy 数组示例图片为 512×512known是 0/1 掩码矩阵像素未被损坏两图相等记为 1被损坏记为 0掩码矩阵known后续既作为约束的乘子也用于从u_corr中“抠出”已知像素值。3.2 建模并求解# 使用总变分修复重建原始图像。 import cvxpy as cp U cp.Variable(shape(rows, cols)) obj cp.Minimize(cp.tv(U)) constraints [cp.multiply(known, U) cp.multiply(known, u_corr)] prob cp.Problem(obj, constraints) # 使用 SCS 求解。 prob.solve(verboseTrue, solvercp.SCS) print(optimal objective value: {}.format(obj.value))代码中三个要素对应数学模型U cp.Variable(shape(rows, cols))重建图像变量形状与图像一致cp.tv(U)TV 原子即目标函数 $\mathop{\bf tv}(U)$cp.multiply(known, U) cp.multiply(known, u_corr)等式约束。cp.multiply是逐元素 Hadamard 乘法见 cvxpy/atoms/affine/binary_operators.py乘上 0/1 掩码后约束等价于“在已知像素处 $U$ 等于原图值在未知像素处 $00$ 无约束”巧妙地将已知像素约束与缺失像素自由度写入同一个等式。prob.solve(verboseTrue, solvercp.SCS)使用SCSSplitting Conic Solver求解并打印详细日志求解结果存入U.value最优目标值由obj.value读取。3.3 SCS 求解日志解读灰度图示例在 512×512 图上运行 SCS v2.0.2 的输出摘要如下Lin-sys: sparse-indirect, nnz in A 1554199, CG tol ~ 1/iter^(2.00) eps 1.00e-05, alpha 1.50, max_iters 5000, normalize 1, scale 1.00 Variables n 523265, constraints m 1045507 Cones: primal zero / dual free vars: 262144 soc vars: 783363, soc blks: 261121 ... Status: Solved Solve time: 2.55e02s ... cx 11044.2661, -by 11044.2813 optimal objective value: 11044.28989542425这些数字揭示出问题的真实规模与结构变量数 $n 523265$512×512 262144 个像素变量加上锥分解引入的松弛变量soc vars 783363 等总计约 52 万约束数 $m 1045507$包含 262144 个 primal zero / dual free 变量对应保真等式约束的松弛以及 783363 个二阶锥SOC分量、261121 个 SOC 块——每个像素对应一个二维 SOC 块正是 $\ell_2$ 范数 $|[U_{i1,j}-U_{ij};,U_{i,j1}-U_{ij}]|_2$ 被 SCS 转换为 SOC 约束的结果求解耗时约 255 秒该示例运行时的硬件环境下的实测值收敛到primal res ≈ 9.0e-06、dual res ≈ 8.2e-06、rel gap ≈ 6.9e-07均在默认容差eps1e-05之下判定Status: Solved最优目标值 ≈ 11044.29即重建图像的 TV 值。文档明确提示这里选择 SCS 是因为它能扩展到比 ECOS 更大的问题规模。ECOS 基于内点法对 52 万变量、105 万约束的锥规划内存开销过大。3.4 展示修复结果与差异图fig, ax plt.subplots(1, 2, figsize(10, 5)) # 展示修复后的图像。 ax[0].imshow(U.value, cmapgray); ax[0].set_title(In-Painted Image) ax[0].axis(off) img_diff 10*np.abs(u_orig - U.value) ax[1].imshow(img_diff, cmapgray); ax[1].set_title(Difference Image) ax[1].axis(off);修复结果直接来自U.valueVariable求解后的数值差异图为10 * |u_orig - U.value|将差异放大 10 倍以便肉眼观察。修复图与原始图几乎一致差异图仅在被文字遮挡的区域显示微弱残影——说明 TV 修复成功抹去了文字而保留了底层图像结构。四、彩色图 TV 修复三通道变量与随机掩码4.1 数学模型扩展彩色图表示为 $m\times n\times 3$ 的 RGB 矩阵 $U^{\text{orig}}$每个像素 $U^{\text{orig}}_{ij} \in \mathbf{R}^3$ 是一个 RGB 向量。TV 定义与灰度版形式相同但每个梯度分量变为三维向量$$\mathop{\bf tv}(U) \sum_{i1}^{m-1}\sum_{j1}^{n-1}\left|\begin{bmatrix} U_{i1,j}-U_{ij} \ U_{i,j1}-U_{ij}\end{bmatrix}\right|_2$$这里向量的每个“分量”本身是三维 RGB 差向量范数取在拼接后的六维向量上对应到 CVXPY 实现则是对 R/G/B 三个通道的梯度统一堆叠取 2-范数见下文源码解析。4.2 随机损坏掩码的构造与灰度版“文字遮挡”不同彩色版通过随机丢弃 70% 像素制造损坏import matplotlib.pyplot as plt import numpy as np np.random.seed(1) # 加载图像。 u_orig plt.imread(data/loki512color.png) rows, cols, colors u_orig.shape # known 为 1 表示像素已知0 表示像素被损坏。 # known 矩阵随机初始化。 known np.zeros((rows, cols, colors)) for i in range(rows): for j in range(cols): if np.random.random() 0.7: for k in range(colors): known[i, j, k] 1 u_corr known * u_orig %matplotlib inline fig, ax plt.subplots(1, 2, figsize(10, 5)) ax[0].imshow(u_orig, cmapgray); ax[0].set_title(Original Image) ax[0].axis(off) ax[1].imshow(u_corr); ax[1].set_title(Corrupted Image) ax[1].axis(off);np.random.seed(1)固定随机种子保证可复现每个像素以 30% 概率保留np.random.random() 0.7三个通道共享同一保留/丢弃决策u_corr known * u_orig被丢弃的像素值置 0损坏图呈黑色斑点状。4.3 三变量建模与求解# 使用总变分修复重建原始图像。 import cvxpy as cp variables [] constraints [] for i in range(colors): U cp.Variable(shape(rows, cols)) variables.append(U) constraints.append(cp.multiply(known[:, :, i], U) cp.multiply(known[:, :, i], u_corr[:, :, i])) prob cp.Problem(cp.Minimize(cp.tv(*variables)), constraints) prob.solve(verboseTrue, solvercp.SCS) print(optimal objective value: {}.format(prob.value))三个矩阵变量R、G、B 各用一个(rows, cols)的cp.Variable放入variables列表逐通道约束对每个通道i用该通道的掩码known[:, :, i]构造保真等式约束cp.tv(*variables)tv原子接受多个矩阵参数将三个通道作为“第三维”统一处理这正是彩色 TV 的实现方式详见下一节源码。求解器同样选择 SCS文档明确指出ECOS 和 CVXOPT 无法扩展到如此大规模的问题彩色版问题规模比灰度版更大。4.4 彩色版求解规模与日志要点WARN: A-p (column pointers) not strictly increasing, column 523264 empty WARN: A-p (column pointers) not strictly increasing, column 785408 empty WARN: A-p (column pointers) not strictly increasing, column 1047552 empty ... Lin-sys: sparse-indirect, nnz in A 3630814, CG tol ~ 1/iter^(2.00) Variables n 1047553, constraints m 2614279 Cones: primal zero / dual free vars: 786432 soc vars: 1827847, soc blks: 261121 ... Status: Solved Solve time: 6.99e02s ... optimal objective value: 11465.652787130613变量数升至约 105 万262144×3 个通道像素变量加上锥松弛约束数约 261 万日志开头的三条WARN: A-p ... column ... empty是 SCS 对稀疏矩阵中空列的提示属于无害告警求解约 700 秒收敛primal res ≈ 9.0e-06、dual res ≈ 9.7e-06、rel gap ≈ 2.4e-06判定Status: Solved最优目标值约 11465.65。4.5 彩色修复结果的可视化import matplotlib.pyplot as plt import matplotlib.cm as cm %matplotlib inline rec_arr np.zeros((rows, cols, colors)) for i in range(colors): rec_arr[:, :, i] variables[i].value rec_arr np.clip(rec_arr, 0, 1) fig, ax plt.subplots(1, 2, figsize(10, 5)) ax[0].imshow(rec_arr) ax[0].set_title(In-Painted Image) ax[0].axis(off) img_diff np.clip(10 * np.abs(u_orig - rec_arr), 0, 1) ax[1].imshow(img_diff) ax[1].set_title(Difference Image) ax[1].axis(off)三个变量的.value分别对应 R/G/B 通道重新拼装成(rows, cols, 3)数组rec_arr用np.clip(rec_arr, 0, 1)将浮点像素值裁剪回合法显示范围plt.imread读入的 PNG 像素归一化到 $[0,1]$差异图同样放大 10 倍并裁剪到 $[0,1]$。视觉上修复图与原始图几乎一致但 RGB 差异图显示许多像素的通道值仍有可见差异——这是随机丢失 70% 像素后重建的必然结果。五、源码纵深CVXPY 中tv原子是如何实现的TV 修复的核心依赖是 CVXPY 的tv原子其完整实现位于 cvxpy/atoms/total_variation.py并通过 cvxpy/atoms/init.py 的from cvxpy.atoms.total_variation import tv导出为cp.tv。5.1 向量与矩阵的分支处理value Expression.cast(value) if value.ndim 0: raise ValueError(tv cannot take a scalar argument.) # 向量使用 L1 范数。 elif value.ndim 1: return norm(value[1:] - value[0:value.shape[0]-1], 1) # 矩阵使用 L2 范数。 elif value.ndim 2: ... else: raise ValueError(tv cannot have input arrays with more than 2 dimensions.)一维向量相邻元素差的一范数 $\text{tv}(x) \sum_i |x_{i1}-x_i|$即 L1 总变分二维矩阵离散梯度的 $\ell_2$ 范数求和即本文 2.2 节的公式标量ndim 0与三维以上输入会抛出ValueError——彩色图正是通过“多个矩阵参数”而非三维数组来绕过该限制。5.2 彩色 TV 的核心多矩阵参数堆叠rows, cols value.shape args map(Expression.cast, args) values [value] list(args) diffs [] for mat in values: diffs [ mat[0:rows-1, 1:cols] - mat[0:rows-1, 0:cols-1], mat[1:rows, 0:cols-1] - mat[0:rows-1, 0:cols-1], ] length diffs[0].shape[0]*diffs[1].shape[1] stacked vstack([reshape(diff, (1, length), orderF) for diff in diffs]) return sum(norm(stacked, p2, axis0))实现分四步取离散梯度对每个矩阵第一个参数value加上*args中的其余通道矩阵计算两个差分切片——mat[0:rows-1, 1:cols] - mat[0:rows-1, 0:cols-1]是右邻差对应公式中 $U_{i,j1}-U_{ij}$mat[1:rows, 0:cols-1] - mat[0:rows-1, 0:cols-1]是下邻差对应 $U_{i1,j}-U_{ij}$按列优先展平每个差分矩阵用reshape(..., orderF)展平成一行length (rows-1)*(cols-1)是差分元素总数堆叠vstack将所有通道、所有方向的差分行堆成一个大矩阵每一“列”恰好包含 $U_{i1,j}-U_{ij}$ 与 $U_{i,j1}-U_{ij}$ 在各通道上的分量逐列取 2-范数并求和norm(stacked, p2, axis0)对每列求 $\ell_2$ 范数灰度图列长为 2彩色图列长为 $2\times36$正是 4.1 节公式的向量化实现最后sum求和。这也解释了求解日志中的锥结构norm(..., p2, axis0)中的每个 $\ell_2$ 范数项都会被 CVXPY 的 DCP 到锥dcp2cone约简转换成一个二阶锥约束最终由 SCS 处理相关约简代码见 cvxpy/reductions/dcp2cone/canonicalizers/pnorm_canon.py。灰度图soc blks: 261121$261121 (512-1)\times(512-1)$与像素差分位置数完全一致可作为实现正确性的交叉验证。5.3 仓库中的自动化测试印证TV 修复流程在仓库测试中有最小化复现cvxpy/tests/test_examples.py 的test_inpainting用 20×20 的随机图像、30% 保留率掩码构造了与本文完全相同的模型U cvx.Variable((rows, cols)) obj cvx.Minimize(cvx.tv(U)) constraints [cvx.multiply(Known, U) cvx.multiply(Known, Ucorr)] prob cvx.Problem(obj, constraints) prob.solve(solvercvx.SCS)该测试作为TestExamples套件的一部分随 CI 运行保证cp.tv、cp.multiply与 SCS 链路长期可用。六、求解器选择为什么用 SCS 而不是 ECOS/CVXOPT求解器算法类型在本问题上的表现SCS一阶算子分裂法ADMM 类内存占用低、可扩展到百万级变量/约束迭代多但每步代价小适合大规模锥规划ECOS内点法IPM精度高、迭代少但需要构造并分解大型稠密/稀疏正规方程52 万变量规模下内存与时间不可承受CVXOPT内点法同 ECOS同样无法扩展到本问题规模文档明确的两条结论灰度版“SCS scales to larger problems than ECOS does”彩色版“ECOS and CVXOPT dont scale to this large problem”。此外从仓库求解器注册表 cvxpy/reductions/solvers/defines.py 可以看到SCS、ECOS、CVXOPT 均为 CVXPY 内置支持的锥求解器用户只需pip install scs等安装对应求解器包即可详见 doc/source/install/index.rst。实用建议追求速度与规模选 SCS或同样一阶的 Clarabel、COSMO 等现代锥求解器仓库中均有对应接口见 cvxpy/reductions/solvers/conic_solvers追求高精度小规模选 ECOS/CVXOPT/MOSEK调试阶段可先用小图如测试中的 20×20快速验证建模正确性再上全尺寸图像。七、完整可运行代码灰度版将以下代码保存为脚本把两张示例图片放到当前目录下的data/子目录即可端到端复现import matplotlib.pyplot as plt import numpy as np import cvxpy as cp # 1. 加载图像 u_orig plt.imread(data/loki512.png) u_corr plt.imread(data/loki512_corrupted.png) rows, cols u_orig.shape # 2. 构造已知像素掩码 known np.zeros((rows, cols)) for i in range(rows): for j in range(cols): if u_orig[i, j] u_corr[i, j]: known[i, j] 1 # 3. 建模 U cp.Variable(shape(rows, cols)) obj cp.Minimize(cp.tv(U)) constraints [cp.multiply(known, U) cp.multiply(known, u_corr)] prob cp.Problem(obj, constraints) # 4. 求解SCS 可扩展至大规模 prob.solve(verboseTrue, solvercp.SCS) print(optimal objective value: {}.format(obj.value)) # 5. 可视化 fig, ax plt.subplots(1, 2, figsize(10, 5)) ax[0].imshow(U.value, cmapgray) ax[0].set_title(In-Painted Image) ax[0].axis(off) img_diff 10*np.abs(u_orig - U.value) ax[1].imshow(img_diff, cmapgray) ax[1].set_title(Difference Image) ax[1].axis(off) plt.show()彩色版仅需把第 3 步替换为 4.3 节的循环建模即可。更详细的完整示例与其余应用案例水填充、信道容量、鲁棒卡尔曼滤波等可在 doc/source/examples/applications 目录下找到。八、总结与延伸阅读本文完整复现了 CVXPY 官方文档的灰度与彩色 TV 图像修复示例从 $\ell_2$ 总变分的数学定义出发推导出以“TV 最小化 已知像素保真约束”为核心的凸优化模型随后给出两套可运行的 CVXPY 代码解读了 SCS 求解日志中的锥结构与规模信息最后深入 cvxpy/atoms/total_variation.py 源码说明tv原子如何通过离散梯度差分、按列堆叠与逐列 2-范数求和将数学公式精确翻译为 DCP 表达式并最终转换为二阶锥规划。延伸方向TV 修复的目标函数可以扩展为加权 TV、各向异性 TV 等变体CVXPY 的tv原子已支持向量 L1 与矩阵 L2 两种形态若需对更大图像更快求解可考虑仓库中同样内置的 Clarabel/COSMO 等一阶求解器cvxpy/reductions/solvers/conic_solvers原始文档位置doc/source/examples/applications/tv_inpainting.rst完整示例索引见 doc/source/examples/applications。赞分享科学计算【免费下载链接】cvxpyA Python-embedded modeling language for convex optimization problems.项目地址https://gitcode.com/gh_mirrors/cv/cvxpy点击查看免费下载相关推荐Kornia 图像去噪实战用可微全变差Total Variation与 PyTorch 优化器实现端到端去噪Kornia 图像去噪实战用可微全变差Total Variation与 PyTorch 优化器实现端到端去噪 图像去噪的目标是在去除噪声的同时尽量保留图像计算机视觉人工智能深度学习图像处理使用 Diffusers 将预训练文生图模型适配为图像修复Inpainting任务实战指南使用 Diffusers 将预训练文生图模型适配为图像修复Inpainting任务实战指南 本文面向希望复用现有 Stable Diffusion 文生图权人工智能媒体生成深度学习音频LaMa实战指南从环境搭建到图像修复全流程详解LaMa实战指南从环境搭建到图像修复全流程详解 引言解决图像修复的效率与质量难题 你是否还在为图像修复任务中遇到的以下问题而困扰修复大尺寸图像时边缘模糊、人工智能计算机视觉深度学习图像处理上一篇OpenModScan源码深度解析C与Qt如何构建跨平台Modbus客户端完整指南下一篇StemDeck DAW式混音深度指南Mute、Solo、循环、波形缩放8个必用的键盘快捷键清单创作声明：本文部分内容由AI辅助生成（AIGC），仅供参考
阅读完成 · 觉得有帮助？
