☰
MATLAB灰度图着色:基于稀疏线性方程组的颜色扩散方法
2026/10/13 11:32:38 网站建设 项目流程

简介:面向图像处理与计算机视觉学习者的MATLAB着色项目,以灰度图像恢复彩色为核心任务,适合需要掌握灰度转换、颜色推理与图像分析流程的开发者参考。项目按开源方式组织,压缩包大小仅891KB,便于快速获取与部署。目前已有138人学习浏览,可作为灰度着色方向的入门与进阶实践参考。内容围绕灰度处理与着色关键技术展开,涵盖rgb2gray灰度变换公式、图像聚类(如imsegkmeans)分割主要区域、高斯滤波与边缘检测等预处理手段,并涉及基于纹理、深度及像素邻域信息的多种着色策略,以及最近邻法、色彩扩散或CNN深度学习模型的实现思路。这些方法分别利用局部纹理特征、深度信息和邻域颜色传播机制完成色彩推理;读者可借此梳理从灰度读取、预处理、算法设计到结果展示的完整流程,也能在开源许可下修改和扩展算法,积累MATLAB图像处理工程经验。

1. 灰度图不只是调滤镜:这是一套能自己扩散颜色的 MATLAB 着色项目

手头如果是上世纪的老照片、扫描存档的灰度医学影像,或者传感器只输出了亮度通道的工业图,难点往往不在“把图变清晰”,而在“让灰度图拥有可信的颜色”。前一阵我在开源社区翻到一个叫 Colorization-project 的 MATLAB 着色项目,把源码拆开后发现,它并没有用现成的调色滤镜,而是靠用户在图上画几笔颜色标记,再让算法通过一组稀疏线性方程组把颜色从标记点扩散到全图。如果你有 MATLAB R2018b 以上版本,想给灰度图做可控着色,同时想理解背后 YUV 空间与邻域权重是怎么配合的,这份源码值得下载后亲手跑一遍。它能解决的核心问题是:在分离亮度与色度的前提下,用可交互的少量标记完成整张图的颜色重建。

2. 算法选型:YUV 空间、邻域权重与稀疏方程组

2.1 为什么先转 YUV,而不是直接在 RGB 上涂颜色

很多第一次接触着色的人会直接把色块画在 RGB 图像上,然后用插值去填周围区域,结果往往是颜色晕染严重,物体边缘出现紫边或绿边。原因在于 RGB 三个通道互相纠缠,亮度信息分散在三个分量里,人为改变某一通道的数值时,很难保证不破坏原有的明暗层次。

这份项目的做法是先把颜色转移到 YUV 空间,Y 通道直接保留原灰度图的亮度,U 和 V 两个色度通道分别代表蓝色差与红色差。因为人眼对亮度变化最敏感,对色度变化相对不敏感,所以把亮度固定住,只求解 U、V 两个通道,既降低了计算量,也让最终结果的明暗关系和原图保持一致。相比 HSV 空间,YUV 没有 H 通道的周期性问题,不会出现颜色从色相环另一端绕回来的情况,在线性代数求解框架里也更稳定。

我在实际使用中一般会把 YUV 转换单独封装成一个函数,输入 RGB 标记矩阵,输出三个分量,这样主流程里只需要关心 U、V 两个待求量。项目中并没有用复杂的深度学习模型,而是沿用经典的优化着色思路,这也是它适合入门和二次修改的原因。

2.2 5×5 邻域权重:颜色是靠什么扩散到周围像素的

颜色从标记点扩散到全图,并不是简单的高斯模糊,而是构建一个像素之间的相似度图。对于灰度图中的任意两个相邻像素 i 和 j,定义它们之间的权重为:

w_ij = exp(-(Y_i - Y_j)^2 / (2 * sigma^2)) / d_ij

其中 Y_i 是像素 i 的灰度值,d_ij 是两个像素的几何距离。这个公式的含义很直白:如果两个像素的亮度接近,比如同属于一块平滑的皮肤区域,那么权重趋于 1,颜色可以顺利扩散过来;如果两个像素之间隔着一条明显边缘,亮度差很大,权重就接近 0,颜色会被阻挡在边界内。sigma 控制亮度差异的缩放尺度,取值过小时,即使很轻微的纹理变化也会阻断颜色扩散;取值过大时,物体边界形同虚设。

代码里默认使用 5×5 邻域窗口,也就是每个像素与周围 24 个像素建立连接。邻域半径 r=2 时效果比较均衡,r 增大到 3 会让颜色传播范围更广,但跨边界的风险也会上升。你可以把权重矩阵打印成热力图看一下,会发现它其实是一张“颜色通道的可通行地图”。

2.3 组装稀疏线性方程组:从邻域权重到矩阵求解

得到所有像素两两之间的权重后,需要组装一个大型线性系统。设 N 为图像总像素数,D 为对角矩阵,对角线上的值是每个像素所有邻域权重之和,W 为权重矩阵,M 为标记位置构成的对角矩阵,那么求解色度通道的方程可以写成:

(D - W + lambda * M) * u = lambda * u_marker

这里 u_marker 是标记位置的 U 通道颜色值,非标记位置为 0。lambda 是标记锚定系数,它决定标记点的约束强度。之所以在左侧加上 lambda * M,是因为如果没有这个东西,整个方程只依赖像素间的相对关系,颜色可能存在任意常数偏移,结果不稳定。加入标记约束后,解被牢牢钉在用户画下的颜色附近,同时仍然允许邻域平滑传播。

核心代码骨架如下,我自己重构过的版本可以跑通,并且保留了稀疏结构避免内存爆炸:

function [A, b] = build_wsystem(I, marks, r, sigma, lam) % 构建着色线性系统 % I 灰度图,double,范围 0~1 % marks 标记矩阵,大小与 I 相同,非标记位置为 0 % r 邻域半径,默认 2 % sigma 亮度差缩放参数,默认 60 % lam 标记锚定系数,默认 0.01 [H, W] = size(I); m = H * W; L = spalloc(m, m, m * (2*r+1)^2); % 预分配稀疏矩阵 idx = reshape(1:m, H, W); % 每个像素的全局线性索引 for dx = -r:r for dy = -r:r if dx == 0 && dy == 0 continue; end xs = max(1, min(W, (1:W) + dx)); % 边界钳制 ys = max(1, min(H, (1:H) + dy)); [Xg, Yg] = meshgrid(xs, ys); j = idx(sub2ind([H, W], Yg(:), Xg(:))); i = idx(:); diffI = I(i) - I(j); wgt = exp(-diffI.^2 / (2 * sigma^2)) ./ max(sqrt(dx^2+dy^2), 1); L = L + sparse(i, j, wgt, m, m); end end D = spdiags(sum(L, 2), 0, m, m); M = spdiags(double(marks(:) ~= 0), 0, m, m); A = D - L + lam * M; b = lam * marks(:); end

这个函数的核心价值在于把双重循环得到的邻域关系写成了稀疏矩阵 L。需要注意 sub2ind 在这里的作用是把二维的行列坐标转成一维索引,因为稀疏矩阵的每一行代表一个像素,列代表该像素的邻域像素。预分配时按照 m 乘以邻域像素总数估算非零元素数量,可以显著减少动态增长带来的性能损耗。

求解时我一般用预条件共轭梯度法 pcg,而不是直接反解矩阵。因为直接求逆的时间复杂度接近 O(N^3),对 30 万像素以上的图像完全不可接受:

% 求解 U 通道示例 u0 = marks_yuv(:, :, 2); % 从 YUV 标记矩阵中取 U 通道 [A, bu] = build_wsystem(gray_img, u0, 2, 60, 0.01); % 预处理提升收敛速度 Lpre = ichol(A, struct('type', 'nofill')); % pcg 迭代求解,容差 1e-6,最多迭代 300 次 u = pcg(A, bu, 1e-6, 300, Lpre); u = reshape(u, H, W);

参数说明:pcg 的第二个参数是右端项 bu,第三个是收敛容差,越小结果越精确但迭代越多;第四个是最大迭代次数,防止极端情况下无限循环;第五个是预处理矩阵。ICHOL 预处理适合这种对角占优的稀疏矩阵,能明显加快收敛。V 通道用完全相同的流程再跑一次。

3. 跑通项目:从标记矩阵到彩色输出的完整实操

3.1 代码结构:从标记文件到 YUV 合成的主链路

拿到开源项目后,建议先看目录结构。这份着色项目的主要文件可以按数据流顺序排成一条线:

文件作用输入输出
run_demo.m演示入口灰度图与标记文件合成彩色图
get_strokes.m交互式画标记用户在图像上点击标记坐标与颜色
make_mark_matrix.m标记矩阵生成坐标与颜色数组三维标记矩阵
build_wsystem.m线性系统构建灰度图、标记、参数稀疏矩阵 A 和右端项 b
colorize.m主着色流程灰度图、标记矩阵彩色图像

这里有一个容易忽略的细节:标记矩阵不是二值图,而是大小与原图一致的三通道矩阵,在标记位置记录颜色值,非标记位置为 0。run_demo.m 的顶层逻辑是先读取灰度图,再读取或生成标记矩阵,调用 colorize.m 后输出结果。

3.2 生成标记矩阵:手工坐标与 GUI 画笔两种方式

最简单的标记方式是在脚本里直接指定像素坐标和颜色,适合快速验证。我准备了一段可以直接粘贴的代码:

% prepare_marks.m % 手工生成标记矩阵:在指定坐标写入指定颜色 gray_img = im2double(imread('data/portrait_gray.png')); [H, W] = size(gray_img); % 三个标记点:第1列为 x,第2列为 y pos = [80, 120; 200, 180; 260, 300]; % 三行对应三个标记点的 RGB 颜色(0~1) color = [0.85, 0.65, 0.55; 0.75, 0.55, 0.45; 0.25, 0.30, 0.35]; marks = zeros(H, W, 3); for k = 1:size(pos, 1) x = round(pos(k, 1)); y = round(pos(k, 2)); % 这里把颜色写入一个 3x3 邻域内,方便后续扩散启动 xr = max(1, x-1):min(W, x+1); yr = max(1, y-1):min(H, y+1); for c = 1:3 marks(yr, xr, c) = color(k, c); end end save('data/mark_map.mat', 'marks');

这段代码的关键点是:MATLAB 的矩阵索引第一维是行 y,第二维是列 x,很多人习惯先写 x 再写 y,结果标记点全部偏移到转置位置。标记点的颜色值会被写入 3×3 邻域,这是为了给扩散过程一个初始“种子区域”。如果你用的是 GUI 方式,get_strokes.m 内部会用类似 roipoly 的工具让你在图上画区域,本质上也是生成这样一个三维矩阵。

3.3 主流程:YUV 分离、双通道求解与彩色合成

当标记矩阵准备好后,主着色脚本的结构非常清晰。先读灰度图,再把标记矩阵从 RGB 转到 YUV,然后分别求解 U、V,最后合成为彩色图像。

% colorize_demo.m % 主流程:灰度图 + 标记矩阵 -> 彩色结果 clear; clc; close all; gray_img = im2double(imread('data/portrait_gray.png')); load('data/mark_map.mat', 'marks'); % RGB 标记转 YUV,只保留 U/V 分量 marks_yuv = rgb2yuv(marks); u_marker = marks_yuv(:, :, 2); v_marker = marks_yuv(:, :, 3); % 构建并求解 U 通道 [H, W] = size(gray_img); [A, bu] = build_wsystem(gray_img, u_marker, 2, 60, 0.01); Lpre = ichol(A, struct('type', 'nofill')); u = pcg(A, bu, 1e-6, 300, Lpre); u = reshape(u, H, W); % 构建并求解 V 通道 [A, bv] = build_wsystem(gray_img, v_marker, 2, 60, 0.01); v = pcg(A, bv, 1e-6, 300, Lpre); v = reshape(v, H, W); % Y 通道沿用原灰度图,合成 YUV 后转回 RGB yuv_out = cat(3, gray_img, u, v); color_out = yuv2rgb(yuv_out); % 拼接显示与保存 imshowpair(gray_img, color_out, 'montage'); imwrite(color_out, 'result_colorized.png');

这里有几个参数值得单独说明。r=2 对应 5×5 邻域窗口,sigma=60 是亮度差缩放系数,在 8 位灰度图归一化到 0~1 后,60 这个值相当于允许约 0.15 的灰度差产生明显权重。lambda=0.01 是标记锚定强度,如果觉得标记点的颜色被周围像素“稀释”得厉害,可以相应调大。pcg 的容差 1e-6 已经足够,肉眼几乎看不出继续减小容差带来的差异。

3.4 结果对比与保存:注意输出色彩空间

合成阶段最容易出错的地方是通道顺序。cat(3, gray_img, u, v) 三个参数分别对应 Y、U、V,不能顺手写成 cat(3, u, gray_img, v)。很多奇怪的整体偏色都源于这种通道错位。我在第一次跑的时候也踩过,后来习惯在 yuv2rgb 函数内部加断言,检查输入第三个维度的尺寸是否都是 H 乘 W,能提前拦截大部分低级错误。

保存结果时建议用 PNG 格式而不要用 JPEG,因为 PNG 是无损压缩,不会在颜色扩散区域附近产生压缩伪影。项目里还有一个 compare_with_ref.m 脚本,可以把结果和参考彩色图并排显示,方便目测颜色还原度。

4. 避坑与排查:四个常见翻车点与解决办法

4.1 标记明明画了,出图却还是黑白原样

现象:标记点在原图上显示正常,但着色结果和灰度图几乎一样,U、V 通道全解出来是零。

原因:最常见的情况是标记矩阵被读成了逻辑类型或 uint8 类型,参与运算时被隐式取整。比如 marks 是 uint8,标记位置写入的是 255,而非 1,在构建 b = lam * marks(:) 时数值是 255 的整数倍,但求解结束后 reshape 出来的数组里大部分值是 0 与 255 之间的浮点数,合成时被意外钳制到 0。另一种可能是你只给 Y 通道画了标记,U、V 通道全是 0,相当于只约束了亮度却没约束色度。

解决:在读取标记后统一执行 marks = im2double(marks),确保数值范围是 0~1。然后打印一下 u_marker 和 v_marker 的非零元素数量,如果少于标记点数,说明通道拆分有问题。我在调试时会在求解前加一行 assert(nnz(u_marker) > 0),没有标记就直接报错。

4.2 颜色越过物体边界,糊成一片

现象:人脸部分的肤色扩散到了背景墙上,或者衣服的颜色渗进了皮肤区域,画面像被水彩晕染过。

原因:sigma 设置过大,导致边缘两侧的亮度差无法有效阻止颜色传播。对于归一化到 0~1 的灰度图,sigma=60 其实偏大,因为大多数自然图像的边缘亮度差在 0.2~0.8 之间,60 的平方值让 exp(-diffI^2 / (2*3600)) 在 diffI=0.5 时仍然有约 0.97 的权重,边缘几乎完全不设防。

解决:把 sigma 调整到适当范围。皮肤与背景的亮度差通常在 0.2 到 0.5 之间,希望边缘处权重降到 0.5 以下,需要让 diffI^2 / (2*sigma^2) 大于 0.69。以 diffI=0.3 计算,sigma 取 30 左右比较合适。也可以先画一条标记线跨过边界,观察哪一个 sigma 值下颜色被阻挡,再反向确定参数区间。

4.3 出图偏灰偏淡,像盖了一层雾

现象:彩色结果有颜色,但整体饱和度明显偏低,看起来不像手绘标记那么鲜艳。

原因:求解过程中只优化了 U、V 通道的平滑性,但标记矩阵中的颜色如果是在 RGB 空间手工指定的,转换到 YUV 后 U、V 的绝对值很小,而 Y 通道保持原灰度不变,最终合成时色度的权重感不足。另一个常见原因是把标记画在了亮度极端的区域,比如纯黑或纯白位置,这些区域 UV 分量天然接近 0。

解决:标记颜色尽量取自真实场景中的中间调区域。如果你在代码里用 [0.85, 0.65, 0.55] 这种肤色值,归一化后 YUV 的 U、V 大约只有 0.05 到 0.1,视觉上就会偏淡。可以把 lambda 从 0.01 调到 0.05,增强标记颜色的强制力,让解更接近标记色。还有一种做法是在合成后做一次饱和度拉伸,例如把 U、V 各乘 1.3 再转回 RGB,但这是事后调整,不如从标记源头解决。

4.4 中等图像也要跑好几分钟,像卡死一样

现象:500×400 的灰度图,跑完 U、V 两个通道需要几分钟,甚至内存占用接近几个 GB。

原因:build_wsystem.m 中如果不小心用了全矩阵存放权重,内存消耗会直接到 O(N^2)。比如 20 万像素的全矩阵需要 200000^2 × 8 字节,约等于 320 GB,机器必然卡死。另一个隐藏原因是循环里反复调用 sparse 没有预分配,导致大量内存重排。

解决:严格使用稀疏矩阵,并且在一开始用 spalloc 预分配。我在重构版本里已经把双重循环中的 sparse 累加保留在 L 上,但如果你发现内存仍然很高,可以把邻域半径 r 从 2 降到 1,这将把非零元素数量从每个像素约 24 个降到约 8 个,矩阵规模直接缩小三分之二。此外,调试阶段先对图像做 0.25 倍降采样,确认参数有效后再跑全分辨率,这是最省时间的策略。

5. 参数调优:邻居半径、亮度缩放与锚定系数的工程取舍

5.1 四个核心参数各管什么

参数默认值控制对象调大后的效果
r 邻域半径2颜色可扩散的物理距离扩散更远,但边缘容易失守
sigma60亮度差异对权重的衰减速度边缘更容易被穿越
lambda0.01标记点对解法的锚定强度标记色更忠实,但可能出现硬边界
标记数量视情况提供颜色约束的信息量收敛更稳定,但交互成本上升

这里最容易误调的是 sigma 和 lambda 的组合。sigma 控制权重缩放,lambda 控制标记约束强度,两者不是重复作用。sigma 负责决定颜色传播的路径,lambda 负责决定标记点本身有多么“不可妥协”。如果 sigma 调小了,颜色很难穿过纹理区域,这时即使把 lambda 调大,也只是把标记点周围一小块区域涂满颜色,无法向远处扩散。

5.2 参数联动:一套比较稳的经验区间

我在多次跑图后形成了一套经验参数:r 固定为 2,sigma 在 20 到 40 之间调整,lambda 在 0.005 到 0.05 之间调整。噪声较大的图像,sigma 适当调小,否则噪点会成为颜色跨越边缘的桥梁。光照渐变的图像,比如黄昏天空,sigma 可以调大一些,让颜色平滑过渡。lambda 则根据标记点的数量反向调整,标记越多,每个点的约束可以越弱,否则会出现明显色块感。

画标记的位置比数量更重要。我在人像上通常只画 5 到 8 笔,分别落在额头、脸颊、嘴唇、头发和背景,而不是均匀撒点。颜色扩散需要的是“锚点”,不是“填色板”。标记点要尽量覆盖图像的主要亮度层次:暗部、中间调、亮部各至少一个,这样方程组的右端项才有足够的信息来约束整体色度解。

5.3 快速扫参流程:先降采样验证,再全分辨率出图

推荐一个高效的参数筛选流程:先把原图缩小到 25% 分辨率,固定 r=2,设置 sigma 和 lambda 各 4 个候选值,组合成 16 组参数,批量跑一遍后把所有结果拼成蒙太奇图。这个过程在 200×150 的小图上只需要几秒钟到十几秒钟。

% quick_scan.m % 极小分辨率快速验证参数组合 small_img = imresize(gray_img, 0.25); small_marks = imresize(marks, 0.25); sigma_list = [20, 30, 40, 50]; lambda_list = [0.005, 0.01, 0.02, 0.05]; montage_imgs = cell(1, 16); k = 0; for s = sigma_list for la = lambda_list k = k + 1; A = build_wsystem(small_img, small_marks(:, :, 2), 2, s, la); u = pcg(A, la * small_marks(:, :, 2)(:), 1e-6, 200); montage_imgs{k} = reshape(u, size(small_img)); end end

扫参后重点判断两个现象:一是物体边界处有没有颜色串扰,二是在大块平坦区域颜色是否过渡均匀。如果 16 组结果里边界都干净,就选饱和度最高的一组;如果某些组出现了边缘溢色,说明 sigma 偏大或 lambda 偏小。确定参数后,再把完整尺寸的图跑一遍,通常 2000×1500 左右的图像在合理参数下能在一分钟左右完成。

6. 进阶技巧:自动采样标记与批量结果验证

当同一批灰度图来自相同场景或相同设备时,手动一一点点标记太费时间,可以写一个自动采样函数:从一张相近构图的彩色参考图中,选几个固定坐标点提取颜色,再把坐标映射到灰度图上的对应位置,批量生成标记矩阵。

% auto_marks.m % 从参考彩色图采样标记色,生成标记矩阵 ref = im2double(imread('data/ref_color.png')); [H, W] = size(gray_img); pos_list = [50, 60; 120, 150; 210, 90]; % 与灰度图对齐的坐标 marks = zeros(H, W, 3); for i = 1:size(pos_list, 1) y = pos_list(i, 2); x = pos_list(i, 1); patch = ref(max(1,y-2):min(H,y+2), max(1,x-2):min(W,x+2), :); marks(y, x, :) = reshape(median(reshape(patch, [], 3), 1), 1, 1, 3); end

采样时最好在彩色参考图中取局部中值而不是单一像素点,因为单一像素容易受传感器噪点影响,局部中值更接近人眼看到的稳定颜色。批量处理时还可以为每张灰度图生成一个着色结果,同时计算与原彩色参考图的 SSIM 指标来判断整体还原质量。合成结果后把灰度图、标记覆盖图、着色结果三张图水平和方向拼接,形成直观的验收板。从那以后,我每次做批量灰度着色时都强制走一遍这个流程:先自动采样标记,再小图快速验证参数,最后全分辨率出图并统计指标。这套办法帮我省掉了大量返工时间,也希望这次的拆解能帮你在自己的灰度图上少踩几个坑。

本文还有配套的精品资源,点击获取

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询