简介面向图像处理学习者、MATLAB 编程用户与课程实验需求的一份凸集投影法超分辨重建代码重点解决低分辨率图像在物理约束下逐步恢复细节的问题。压缩包共 1 个文件类型为 .m 脚本大小约 1KB轻量直接适合阅读算法主干、嵌入小规模验证或作为课堂演示素材。内容围绕初始化、投影算子、逐约束迭代与停止条件组织并涉及非负性、能量和频谱约束以及频域扩展、正则化抑制噪声等思路读者可据此拆解投影函数写法、迭代顺序与阈值判定观察参数变化对收敛速度和重建清晰度的影响。已有 3872 人学习便于后来者复现实验、核对参数设置并积累图像恢复与排错经验。1. 从一段抖动的连拍说起POCS 超分辨重建能解决什么同一台相机对着静止场景连拍十几张画面里其实藏着肉眼看不见的亚像素位移凸集投影法POCS超分辨重建就是利用这部分信息的算法。它把每一张低分辨率帧看成对高分辨率图像的一个约束也就是一个凸集重建结果必须同时落在所有这些凸集的交集里。于是算法反复在各个集合上做投影逐步把高分辨率图像修正出来。这套 MATLAB 实现的好处在于插值、滤波、FFT 都是现成函数一两百行就能跑通一版能看的。适合手上有一组同场景连拍或视频抽帧、想验证重建效果的人做成像仿真的学生、做工业检测和遥感预处理的工程师以及准备把算法接入自己 MATLAB 流水线、需要先跑通基线再改的人。前提只有一个帧与帧之间真的存在亚像素级差异否则算法无法恢复出额外细节。2. 工程结构与运行在 MATLAB 里把 POCS 跑起来动手前先明确一件事POCS 不是一个函数而是一条链子。读帧、配准、建模、投影、评估任何一环出问题最后都会表现为“结果糊”或“结果有网格”而原因往往并不相同。所以工程组织的第一原则是把每一环拆成独立文件参数集中到一个结构体里出问题时能单独替换某一环做对照实验。2.1 环境依赖与版本注意点运行环境建议 MATLAB R2018b 及以上必须装 Image Processing Toolbox因为imresize、imfilter、fspecial、psnr、ssim都在这个工具箱里。如果帧序列在 20 张以上、放大倍数又设到 4迭代会明显变慢这时 Parallel Computing Toolbox 可以用于并行计算残差但投影更新本身是串行的不能整体parfor。有一个问题需要提前处理fspecial(gaussian, ...)在新版本 MATLAB 里被标记为不推荐部分版本会给出警告未来可能被移除。稳妥做法是自己用meshgrid生成高斯核或者用imgaussfilt直接对图像做高斯滤波。我一般选前者因为超分辨里需要显式拿到 PSF 矩阵去做转置滤波直接调滤波函数反而不方便function h gauss_kernel(sigma) % GAUSS_KERNEL 生成归一化二维高斯核替代 fspecial(gaussian) % sigma 为模糊标准差单位是像素核尺寸取 6*sigma1覆盖 ±3sigma sz 2*ceil(3*sigma) 1; ax -(sz-1)/2 : (sz-1)/2; [xx, yy] meshgrid(ax, ax); h exp(-(xx.^2 yy.^2) / (2*sigma^2)); h h / sum(h(:)); % 必须归一化否则重建结果整体亮度会漂 end这段的关键在最后一行。高斯核不归一化时正向模型T不再是能量守恒的算子投影回去的残差会带着系统性偏差迭代几十轮后整幅图会逐渐变亮或变暗而且这种漂移很均匀容易被误判成“幅度约束没生效”。另一个版本相关点是interp2的插值方法。spline精度高但会在强边缘附近产生过冲cubic稳一些但锐度略差。超分辨这种毫米级较真的场景我倾向spline配幅度裁剪把过冲直接截断掉。2.2 目录结构与各文件职责一份能二次开发的 POCS 工程我一般按下面这样切。文件不多但每个都能单独拿出来测试文件职责关键参数demo_pocs_sr.m入口脚本串起读图、重建、保存无default_config.m所有可调参数的唯一来源全部load_frames.m读序列、统一尺寸、转 double 归一化灰度化开关estimate_shift.m相位相关法亚像素配准窗函数类型gauss_kernel.m生成归一化高斯 PSFpsfSigmaimshift.m亚像素平移interp2实现边界模式degrade.m正向降质平移→模糊→下采样放大倍数Lpocs_sr.mPOCS 主迭代lambda/delta/toleval_metrics.mPSNR / SSIM 对比评估参考图路径把degrade.m单独拆出来是这套结构里最值钱的一步。它既是正向成像模型又是验证工具拿一张已知的高分辨率图用degrade人工造出低分辨率序列再送进pocs_sr如果重建结果比双三次插值还差那问题一定在配准或投影而不在数据。这条“自造数据自测”的路径能省掉大量猜测后面第 6 章还会用到。2.3 最小可跑流程四行代码串起整条链路入口脚本刻意写得很短所有决策都推给default_config这样换数据集时只改一个文件% demo_pocs_sr.m —— 最小可跑入口 clear; clc; close all; addpath(genpath(fileparts(mfilename(fullpath)))); cfg default_config(); % 参数集中管理 [LR, refIdx] load_frames(cfg); % 读入低分辨率帧序列cell 数组 fprintf(读入 %d 帧尺寸 %dx%d\n, numel(LR), size(LR{1},1), size(LR{1},2)); X pocs_sr(LR, refIdx, cfg); % 核心重建 bic imresize(LR{refIdx}, cfg.scale, bicubic); % 基线对照 figure; subplot(1,3,1); imshow(LR{refIdx}); title(低分辨率参考帧); subplot(1,3,2); imshow(bic); title(双三次插值); subplot(1,3,3); imshow(X); title(POCS 重建); imwrite(im2uint8(X), cfg.outPath);LR用 cell 数组而不是三维矩阵是因为后续可能混入不同尺寸的帧比如视频抽帧时做过裁剪用 cell 能容纳这种情况。refIdx默认取第 1 帧实际使用时建议选清晰度最高、运动模糊最小的一帧做参考因为它是高分辨率网格的初值来源参考帧模糊后面所有投影都从模糊的状态开始修正。参数文件长这样每个都写了取值区间照着调就行function cfg default_config() % DEFAULT_CONFIG POCS 超分辨的全部可调参数 cfg.scale 2; % 放大倍数 L2 或 4 比较稳 cfg.psfSigma 1.2; % 高斯 PSF 标准差按实际光学模糊估 cfg.lambda 0.35; % 松弛因子0.1~0.5太大震荡太小慢 cfg.deltaScale 1.5; % 数据一致性阈值系数 cfg.noiseSigma 0.01; % 归一化到 [0,1] 后的噪声标准差 cfg.maxIter 40; % 外迭代次数20~50 足够 cfg.tol 1e-4; % 残差相对变化小于它提前退出 cfg.range [0 1]; % 幅度约束区间 cfg.refIdx 1; % 用第几帧做初值 cfg.dataDir data/frames; cfg.outPath result/sr_out.png; end这里最容易随意填的是psfSigma。它不是“随便给个 1”而是应该由成像系统的实际模糊量决定如果你的低分辨率图本身就是清晰的、只是被降采样那 PSF 应该接近一个 δ 函数psfSigma取 0.5 左右如果是相机镜头拍出来的通常会大于 1。这个值偏大重建图会过锐、边缘出现白边偏小则怎么迭代都不够锐。常用的标定办法是拿一张理想边缘图看边缘过渡占几个像素sigma ≈ 过渡宽度 / 2.5。3. 正向模型与配准成不成八成看这两步很多 POCS 实现最后效果不行问题并不在投影公式写错了而是正向模型和实际成像对不上或者配准精度只到整像素。这一章把这两块拆开讲因为它们是后面所有迭代的地基。3.1 把“我要反演什么”写成算子超分辨的出发点是这样一个假设手上的低分辨率帧y_k是高分辨率图像X经过一系列已知的物理过程得到的。完整链条是三步顺序不能乱把X反向平移s_k对齐到第 k 帧的视角卷上光学系统的点扩散函数h模拟模糊按倍数L隔点采样模拟探测器。写成算子就是y_k D · H · S_k · X其中S_k是平移、H是卷积、D是下采样。POCS 迭代中每一步都要用这个正向算子把当前的高分辨率估计“走一遍”和真实低分辨率帧比残差function y degrade(x, h, L, dy, dx) % DEGRADE 正向成像模型平移 - 模糊 - 下采样 % x 高分辨率估计double, 已归一化 % h 归一化高斯 PSF % L 下采样倍数 % dy/dx 第 k 帧相对参考帧的位移行/列方向 xs imshift(x, -dy, -dx); % 反向平移站到第 k 帧的位置 xb imfilter(xs, h, replicate, same); % 光学模糊边界复制不引入黑边 y xb(1:L:end, 1:L:end); % 隔点采样等价 L 倍下采样 endimfilter而不是conv2是因为conv2的same默认零填充图像四周会被当成“外面是黑的”模糊之后边界一圈系统性偏暗投影回去就形成一圈固定的暗环。replicate复制边缘像素能把这个偏差压到最小。xb(1:L:end, 1:L:end)这个写法对应的是“点采样”模型也就是每个低分辨率像素只取高分辨率网格上的一个点。真实传感器是面积积分更准确的写法是对L×L块取均值。两种模型在L2时差别不大L4时点采样会让重建结果偏锐、带一点混叠如果发现结果里有细密的高频纹理把这里改成块均值再试一次% 面积积分版本替换上面的隔点采样 y zeros(size(xb,1)/L, size(xb,2)/L); for i 1:L for j 1:L y y xb(i:L:end, j:L:end); end end y y / (L*L);代价是慢一些但和真实相机的成像更接近。选哪个取决于你的数据来源仿真数据用点采样真实相机用面积积分。3.2 亚像素配准相位相关加抛物线插值位移s_k估不准后面全白搭。整像素精度的配准会让重建结果退化成“把参考帧插值放大”因为亚像素信息根本没被利用上。这里用相位相关法它对亮度变化不敏感计算量靠 FFT还能直接给出亚像素修正。function [dy, dx] estimate_shift(ref, mov) % ESTIMATE_SHIFT 相位相关法估计 mov 相对 ref 的亚像素平移 % 输出 dy/dx 为行/列方向位移正号表示 mov 中的内容向右下移动 F fft2(ref) .* conj(fft2(mov)); F F ./ (abs(F) eps); % 相位归一化只保留相位信息 cc real(ifft2(F)); % 互相关面峰值位置即整像素位移 [~, k] max(cc(:)); [py, px] ind2sub(size(cc), k); dy (py-1) parabola(cc(:, px), py); % 行方向亚像素修正 dx (px-1) parabola(cc(py, :)., px); % 列方向亚像素修正 % FFT 位移是循环的落在 [N/2, N) 的真实位移是负的要折回来 if dy size(cc,1)/2, dy dy - size(cc,1); end if dx size(cc,2)/2, dx dx - size(cc,2); end end function d parabola(v, i) % 三点抛物线拟合用峰值左右各一点估出亚像素偏移量 if i 1 || i numel(v), d 0; return; end a v(i-1); b v(i); c v(i1); den a - 2*b c; if abs(den) 1e-12, d 0; else, d 0.5*(a - c)/den; end endF ./ (abs(F) eps)这一行是相位相关的核心。不做归一化就是普通互相关图像亮度整体变化会导致峰值偏移归一化之后只看相位对曝光差异免疫。eps是防零除别省。抛物线的三点拟合假设相关峰附近近似二次函数在峰值明显的情况下误差在 0.1 像素以内够用。如果序列噪声大、相关面比较平可以升级成五点拟合或者上采样相关面再做插值代价是计算量翻几倍。还有一个隐含约定容易搞混imshift(x, -dy, -dx)里的负号。estimate_shift返回的是mov相对ref的位移而degrade是要把高分辨率图从“参考帧坐标系”搬到“第 k 帧坐标系”方向相反所以传负值。这个符号如果写反正向模型就变成了反向平移残差会大得离谱且怎么迭代都降不下来很容易被误判成参数没调好。3.3 配准结果的自检三个数先看一眼配准完不要急着进迭代先打印三个数位移量级、相关峰值的旁瓣比、以及配准后两帧的残差均值。位移量级如果所有dy/dx都接近 0说明帧之间根本没有亚像素抖动或者你的序列是从同一张图复制出来的这种情况下 POCS 不可能带来增益最多和插值持平。峰值旁瓣比max(cc(:)) / mean(abs(cc(:)))大于 8 通常说明配准可信低于 4 就要怀疑帧间有旋转或尺度变化纯平移模型不适用。残差均值用估计的位移把两帧对齐后做差均值应该显著小于不对齐时的差。如果对齐前后几乎一样说明位移估错了。三个数都是几十毫秒的事情比跑完 40 轮迭代再盯着结果图发呆要划算得多。这也是我一般建议在pocs_sr开头加一段if cfg.debug, ... end的原因。4. 投影算子与主循环POCS 迭代的代码骨架到这一步正向模型有了位移有了剩下的就是“怎么把残差变成对高分辨率估计的修正”。POCS 的名字里“凸集”和“投影”两个词说的是同一件事的两面。4.1 三个凸集各管一件事标准的 POCS 超分辨至少用三个约束凸集定义管什么数据一致性集C_D{X : |y_k - T_k X| ≤ δ}让重建结果重新降质后接近真实低分辨率帧幅度约束集C_A{X : X ∈ [0, 1]}防止迭代中出现负值或过曝白点能量约束集C_M{X : sum(X) S}防止整体亮度漂移可选δ是这个算法里最需要理解的参数。它代表“允许多大的残差”理论值跟噪声标准差同量级。δ 0意味着强制每一帧的数据完全一致这在无噪声的仿真数据上没问题但真实数据里等于把噪声也当成信号去拟合结果就是高频噪点被逐帧放大最后出现明显的颗粒感甚至棋盘格。工程上取δ (1~3) × σ_noise宁可松一点。幅度约束看着平淡实际很关键。投影修正的残差e有正有负累加几十轮后边缘附近很容易冲出[0,1]如果不裁掉下一轮degrade里的插值会产生更极端的值形成正反馈。裁剪是免费的稳定性。4.2 主循环一行投影公式对应一段代码核心更新式是X ← X λ · T_k^T[ sign(r) · max(|r| - δ, 0) ]。逐项翻译成代码function X pocs_sr(LR, refIdx, cfg) % POCS_SR 凸集投影法超分辨重建主循环 ref LR{refIdx}; X imresize(ref, cfg.scale, bicubic); % 参考帧双三次上采样做初值 X min(max(X, cfg.range(1)), cfg.range(2)); h gauss_kernel(cfg.psfSigma); delta cfg.deltaScale * cfg.noiseSigma; % 数据一致性阈值 shifts estimate_all_shifts(LR, refIdx); % 一次性算好所有帧的位移 prevE inf; for it 1:cfg.maxIter energy 0; for k 1:numel(LR) ysim degrade(X, h, cfg.scale, shifts(k,1), shifts(k,2)); res LR{k} - ysim; % 当前帧的残差 energy energy sum(res(:).^2); % 软阈值只把超出 δ 的那部分残差投回去噪声不被放大 res sign(res) .* max(abs(res) - delta, 0); % 残差升采样回高分辨率网格再经 PSF 转置滤波高斯核对称同一个核 e imresize(res, cfg.scale, bicubic); e imfilter(e, h, replicate, same); X X cfg.lambda * e; X min(max(X, cfg.range(1)), cfg.range(2)); % 幅度约束 end energy energy / numel(LR); if cfg.verbose fprintf(iter %2d energy%.6f\n, it, energy); end if abs(prevE - energy) / (prevE eps) cfg.tol break; % 残差不再下降提前退出 end prevE energy; end end几个位置值得逐个说明。sign(res) .* max(abs(res) - delta, 0)是软阈值也是“投影到C_D”这个动作的具体实现。它的意思是残差在±δ以内的部分视为噪声直接丢弃超出的部分才认为是模型和真实数据的不一致需要修正。去掉这一行、直接用res去更新就是经典的“硬一致性反投影”在仿真数据上收敛飞快在真实数据上几轮之后就开始起噪点。imresize(res, cfg.scale, bicubic)是D^T即下采样的转置。严格意义上转置应该用“插零 低通”这样才是精确的伴随算子重建结果的重建误差最小。imresize是近似好处是平滑、不会引入插零带来的高频镜像工程上更稳。我试过两种imresize版本在L2时 PSNR 只低 0.1~0.2 dB但迭代过程平顺很多不会在第 5 轮出现一次明显的能量跳变。cfg.lambda松弛因子控制步长。取 1 是完全投影收敛快但容易在两个凸集之间来回跳取 0.1~0.5 是欠松弛慢但稳。默认 0.35 是个折中如果发现迭代曲线锯齿明显就往下调。帧的遍历顺序也有讲究。默认按k顺序走实际上按位移量级从小到大排更好先处理位移小的帧高分辨率网格和它最接近投影修正幅度小、更准再逐步处理位移大的。这样做相当于从易到难的松弛策略通常能少迭代 5~10 轮。4.3 收敛判据与参数标定tol判据看的是残差能量的相对变化而不是绝对残差。因为不同数据集的绝对残差量级差很多用相对量更通用。1e-4是个经验值太松会提前停在一个还没收敛的状态太紧基本跑满maxIter白等。真正的调参顺序我一般是这样先固定lambda0.35、delta1.5σ、maxIter40跑一遍看迭代曲线。曲线一路平缓下降、最后趋平说明参数合理曲线前几轮掉得很快然后开始上升说明delta太小或者lambda太大曲线几乎不降八成是位移估错了或者 PSF 差太远。标定psfSigma最省事的办法是用合成数据扫一遍。取一张高清图用不同的sigma做degrade造出低分辨率序列再建一张sigma的查找表跑重建看哪个sigma输入能得到最高 PSNR。这个扫描一次十几分钟但能让你对“我的数据对应的合理 PSF 大概是多少”有个量化的感觉比凭肉眼猜强太多。5. 避坑与排查跑不出效果时先看这几处下面这些坑我在不同的数据集上或多或少都遇到过。症状不同但都能在十分钟内定位前提是知道该看哪里。5.1 五个高频翻车现场现象一重建结果和双三次插值几乎一模一样细节毫无提升。原因配准退化。最常见的是所有位移都估成了 0或者帧序列本身就是同一张图复制出来的比如为了凑数把单帧重采样成多帧。这种情况下 POCS 的所有凸集约束几乎是重复的交集就是插值本身自然没有增益。解决打印shifts矩阵看位移的分布。如果全部小于 0.05 像素先确认数据源是否真的有多帧独立采集。另一种可能是相位相关被图像的大面积平坦区域主导峰值不明显这时可以对图像做一次高通预滤波减去imgaussfilt的结果再配准。现象二迭代十几轮后出现规则网格纹理越迭代越明显。原因delta设成了 0 或者极小硬一致性反投影把噪声当信号反复拟合。网格纹理的周期通常等于放大倍数L这是降采样算子D的周期性在残差上的体现。解决把delta提到1.5σ以上同时对残差做软阈值4.2 的代码里已经包含。如果还想更稳可以在每 5 轮投影后做一次轻度平滑或者直接降lambda。现象三图像四周一圈明显发暗或发亮。原因imshift用了零填充imfilter也没指定边界模式边界像素在正向模型里被系统性地低估残差在边界处永远是同一个符号投影几十轮就积出一圈亮/暗边。解决imshift里interp2的越界值改成复制边缘或者干脆用replicate模式的相关函数imfilter加replicate。更彻底的做法是维护一个边界掩膜投影时对离边界不足 PSF 半径的像素降权。现象四重建结果整体偏亮或偏暗细节是对的。原因PSF 没归一化或者下采样用了块求和但 PSF 是按点采样标的两条链的能量守恒关系对不上。解决先检查sum(h(:))是否等于 1再检查degrade里下采样有没有引入L²的能量变化。块均值版本会把能量压低L²倍必须除以L*L补回来第 3 章的代码里已经有这一步。现象五跑到一半内存爆掉或者一晚上没跑完。原因大图 多帧 每轮重新生成 PSF 和重复imresize。imfilter对double类型的整幅图做卷积内存占用是原图的几倍。解决PSF 和位移在循环外算好代码里已经这样做了把图像转成single内存直接减半精度损失在超分辨这种场景可以忽略帧数超过 30 时考虑分块处理或者抽稀。还有一点常被忽略imresize的bicubic在每次内层循环都调用一次如果res尺寸固定可以预先算好插值索引复用省下的时间比想象中多。5.2 一套十分钟的排查顺序出问题时不要一上来就改参数按这个顺序走第一步用合成数据验证整条链路。取一张高清图degrade造 8 帧低分辨率序列跑pocs_sr。如果合成数据上 PSNR 能超过双三次插值 1 dB 以上说明代码没问题问题在你的真实数据如果合成数据也上不去问题在代码或参数。第二步看迭代曲线。能量不下降就是正向模型和配准的问题跟投影公式无关。第三步单独拆出degrade和estimate_shift做单元测试。取高分辨率图用已知位移(0.3, 0.7)做degrade再把结果送回estimate_shift看返回的位移是不是接近(0.3, 0.7)。这一步能直接把“配准不准”和“投影不收敛”两类问题分开。第四步才动参数。按delta → lambda → psfSigma的顺序调每次只改一个改完记录 PSNR。三个参数一起调最后你永远不知道是哪个起了作用。提示把每次实验的cfg、PSNR、迭代轮数写进一个 CSV跑够二十组之后你会对“什么数据配什么参数”形成直觉。这个东西比任何调参攻略都管用。6. 二次开发把 POCS 接进自己的流程这套代码真正价值在于它的可拆解性而不只是“能跑出一张更清晰的图”。三个改造方向我按性价比排一下。第一个是自适应阈值改动最小、收益最直接。把固定的delta换成按当前残差估的噪声% 放在 pocs_sr 的内层循环里每轮更新一次 r res(:) - median(res(:)); sigmaHat 1.4826 * median(abs(r)); % MAD 稳健估计 delta cfg.deltaScale * sigmaHat; % 每轮自适应MAD 比标准差抗异常值图像边缘处的残差不会被误当成噪声。改完之后你会发现同一套参数在噪声水平差异很大的两组数据上都能用不用每次手动调delta。第二个是加正则。POCS 本身对噪声敏感除了软阈值还可以在每 N 轮之后做一次轻度去噪作为额外的凸集约束。最省事的做法是用各向异性扩散或者非局部均值跑一两次就够别做太重否则会把刚恢复的高频细节一起抹掉。重一点的方案是把整个问题改写成min ||y - Tx||² λ·TV(x)用 ADMM 解代码量翻三倍但 PSNR 通常还能再提 0.5~1 dB适合做算法对比实验。第三个是换正向模型。如果你的低分辨率帧不是简单的平移关系比如手持拍摄带轻微旋转或者序列来自不同传感器那degrade里就得把平移换成仿射或者单应矩阵estimate_shift也要换成特征点匹配加fitgeotrans。这种改动工作量大但结构上不需要动主循环——这正是把degrade单独拆成文件的好处。验证方法就一条永远用合成数据说话。我现在的习惯是任何改动先在“高清图 → 人工降质 → 重建”的闭环上跑一遍PSNR 和 SSIM 都要看因为 PSNR 高的版本有时候 SSIM 反而低那是过锐了。两个指标同向提升才算真的改好了。从那以后我每次动投影公式都强制走一遍合成数据的闭环哪怕只是改了一个符号——血泪经验符号写反的时候人眼是看不出来的只会觉得“今天参数怎么调都不对”。希望帮到你。本文还有配套的精品资源点击获取