简介：这份代码与数据包聚焦多帧图像超分辨率重构，面向具备图像处理基础、希望从理论走向MATLAB实现的学习者和开发者，着重解决低分辨率多帧图像融合、恢复与质量评估的落地问题。压缩包共399个文件、约15.38MB，以大量低分辨率与参考位图构成实验数据，配合矩阵数据文件保存中间结果，核心脚本与函数文件对应算法实现，另有图形、文档和说明素材，便于复现实验并理解处理流程。目前已有1381人学习下载，具备一定参考价值。读者可围绕多帧融合、去噪去模糊、超分辨率模型构建、训练预测以及峰值信噪比和结构相似度评估等环节，对照代码梳理实现逻辑，并借助现成图像数据快速跑通流程、观察结果，适合课程设计、项目实践与算法入门研究。
1. 从一组抖动连拍说起：多帧超分不是「多叠几张」
暗光下手机连拍十来张，每张都糊。把十张图直接平均，得到的是更糊的一张。这是很多人对多帧图像超分辨率重构的第一层误解——这类方法的关键在帧与帧之间那零点几个像素的错位，而非「多」这个字。这些错位里藏着单帧已经丢掉的采样信息，只有把它们对齐到一张更密的高分辨率网格上，亚像素的多样性才会变成真实的细节。
这份代码走的是传统路线，不依赖任何预训练权重：亚像素配准 → shift-and-add 融合到高分辨率网格 → 迭代反投影修正，另外配了合成数据生成器（有真值，能算 PSNR/SSIM）和真实连拍序列入口。适合两种人：想把多帧超分每一环都摸清楚的开发者，以及需要一条轻量基线去对比深度学习方案的实验者。关键词里那句「绝对能运行」，靠的就是合成数据那条链路——真值已知，结果可验。
2. 技术选型与目录骨架：为什么是「配准＋融合＋反投影」
2.1 三条超分路线的取舍
单帧超分大致三条路。第一条是插值类（双三次、Lanczos），一秒出图，但放大四倍后边缘是软的，纹理是糊的，它只是把像素摊开，没有增加任何信息。第二条是深度学习类，效果最好，代价是要权重、要显存、要配对训练集，调不动的时候你连它为什么糊都看不出来，是个黑匣子。第三条就是这份代码走的多帧传统路线：靠同一场景的多次采样补齐高频，数学上干净，每一步都能单独打印出来看。
选第三条的原因很实际。多帧超分的核心矛盾是「配准精度」——位移估计差 0.2 个像素，融合出来就是双重边缘，再迭代也救不回来。如果整套流程包在一个网络里，你没法判断到底是配准错了还是融合错了。传统管线把这两件事拆成两个文件，中间结果能存盘、能可视化、能单独喂已知位移做单元验证，这对二次开发的价值远大于那点峰值指标。
所以在工程上我一般这样切分职责：register.py只负责吐位移，fuse.py只负责把帧搬到 HR 网格上做加权累加，backproject.py只负责迭代修正，degrade.py是正向退化模型，既用来造合成数据，也用来做反投影里的模拟降质。这四块互相之间只传数组，不传状态，换掉任何一块都不影响其他三块。
2.2 依赖清单与目录骨架
依赖刻意压到最少：numpy做数值，opencv-python做配准和重采样，scikit-image只用来算 PSNR/SSIM，pyyaml读配置，matplotlib出对比图。没有 CUDA，没有编译环节，纯 CPU 跑 512×512 的合成序列在十几秒量级。
目录结构如下，按「入口 / 库 / 数据 / 输出」四层分开：
mfsr/ ├── run_demo.py # 合成数据一键跑通：造帧 -> 配准 -> 融合 -> IBP -> 出图 ├── run_real.py # 真实连拍序列入口，需要自己提供位移或用配准估计 ├── mfsr/ │ ├── __init__.py │ ├── degrade.py # 正向退化：亚像素平移 + 高斯模糊 + 下采样 + 噪声 │ ├── register.py # 相位相关亚像素配准 │ ├── fuse.py # shift-and-add 加权融合到 HR 网格 │ ├── backproject.py # 迭代反投影 IBP │ ├── metrics.py # PSNR / SSIM │ └── io_utils.py # 读写、YCrCb 拆通道、结果拼图 ├── configs/ │ └── demo.yaml # scale / iters / lam / sigma 全在这 ├── data/ │ ├── synthetic/ # 合成低分辨率帧 + gt.png │ └── real/ # 真实手持连拍序列 └── outputs/ # 每次运行一个时间戳子目录把配置抽到 yaml 而不是 argparse 一把梭，是踩过坑之后的习惯。多帧超分要调的东西至少有六七个（倍数、帧数、迭代数、步长、高斯核、边界模式），命令行参数一长串贴来贴去容易漏，而实验记录必须能原样复现。现在outputs/下每次运行都会把当时用的demo.yaml复制一份进去，回头看结果不对，先看配置副本。
2.3 输入数据的组织方式
真实连拍序列和合成序列的输入结构要统一，否则两套入口会分裂出两套逻辑。约定是这样的：一个目录下若干张同尺寸灰度图，文件名排序即帧序，frame_000.png为参考帧；如果目录里存在gt.png，就跑指标，没有就只出图。位移有两种来源，一是从shifts.txt读（合成数据用，格式是每行dx dy，单位是高分辨率像素），二是现场用相位相关估计（真实数据用）。
单位这件事必须一开始就钉死：所有对外暴露的位移都以高分辨率像素为单位，映射到 HR 网格时直接用，只有配准模块内部在低分辨率图上估出来的值需要乘scale。早期版本没统一，一半代码乘了 scale 一半没乘，融合出来的图整体偏移小半个像素，PSNR 反而比双三次还低，查了半天才定位到符号和单位上。
彩色图的处理方式是常见做法：转 YCrCb，只对 Y 通道做多帧重构，Cb/Cr 用双三次放大后合并。人眼对亮度细节最敏感，色度通道放大四倍也没人看得出来，硬要三通道全跑一遍纯属浪费三倍时间。
3. 核心链路拆解：一帧低分辨率图怎么落到高分辨率网格
3.1 亚像素配准：相位相关怎么挖出 0.3 个像素
配准是整个管线里最不能出错的一环。常见做法有光流（稠密但慢，对亮度变化敏感）、特征点匹配（稀疏，亚像素精度靠拟合）、相位相关（快，对纯平移严格成立）。多帧超分的场景通常是手持连拍，帧间主要就是平移加轻微旋转，相位相关足够用，而且它能直接给出亚像素结果，不用二次拟合。
import cv2 import numpy as np def estimate_shift(ref, mov, window=None): """估计 mov 相对 ref 的亚像素平移，返回 (dx, dy, response)。 ref / mov: 同尺寸单通道 float32 图，建议先减均值做亮度归一化 window: 汉宁窗，抑制 FFT 的边界泄漏；不传则自动生成 response: 相关峰的值，越接近 1 越可信，低于 0.2 基本可以判定配准失败 """ if window is None: # 尺寸参数是 (width, height)，别传反了 window = cv2.createHanningWindow(ref.shape[::-1], cv2.CV_32F) ref32 = ref.astype(np.float32) mov32 = mov.astype(np.float32) # 亮度归一化：连拍时自动曝光会让帧间整体亮度漂移，不归一化会污染相关峰 ref32 = ref32 - ref32.mean() mov32 = mov32 - mov32.mean() (dx, dy), response = cv2.phaseCorrelate(ref32, mov32, window) return dx, dy, response逻辑上分三步：先减均值，去掉帧间整体亮度差这个干扰项；再乘汉宁窗，把图像边界渐变到零，否则 FFT 会认为图像在边界处发生了剧烈跳变，相关面上会冒出假的次峰；最后做相位相关，取峰值位置作为位移。参数上唯一需要注意的是createHanningWindow的尺寸顺序是(宽, 高)，而shape给的是(高, 宽)，所以这里要[::-1]反转一次——这个细节错了不会报错，只会让窗的形状和图像不匹配，配准精度莫名其妙地掉。
response这个返回值容易被忽略，但它是最省事的质量门禁。我在真实序列上跑的时候，会先打印所有帧的 response，正常帧一般在 0.5 以上，个别运动模糊严重的帧会掉到 0.1 附近，这种帧直接剔除比强行融合更好。
符号约定必须实测。cv2.phaseCorrelate返回的是「把mov对齐到ref需要移动的量」还是反过来的，不同版本、不同调用顺序下容易搞混，且搞反了不会报错，只会让融合图越来越糊。我一般的做法是造一组已知位移的合成帧，跑一遍配准，看返回值和设定值是否同号，两分钟的事，能省掉后面一整天的自我怀疑。
3.2 高分辨率网格映射：加权累加为什么要带 mask
融合这一步的核心动作是：把每张低分辨率帧按它的位移放大scale倍后搬到 HR 网格上，然后叠加。听起来是简单的warpAffine加平均，但直接平均会出问题——warpAffine的边界默认填 0，边缘一圈会被拉黑，帧数越多黑边越明显，而且因为每帧位移不同，黑边位置还会漂移，平均之后就变成一圈脏兮兮的暗带。
def shift_and_add(frames, shifts, scale=4, interp=cv2.INTER_LINEAR): """把多帧低分辨率图按亚像素位移累加到高分辨率网格上。 frames: [H, W] float32 列表 shifts: [(dx, dy), ...] 单位是高分辨率像素（就是 register 结果 * scale） scale: 放大倍数 """ h, w = frames[0].shape H, W = h * scale, w * scale acc = np.zeros((H, W), np.float32) # 加权累加器 wsum = np.zeros((H, W), np.float32) # 权重累加器 ones = np.ones((h, w), np.float32) # 用一个全 1 图同步 warp，得到的即权重 for img, (dx, dy) in zip(frames, shifts): M = np.array([[1, 0, dx], [0, 1, dy]], np.float32) warped = cv2.warpAffine(img, M, (W, H), flags=interp, borderMode=cv2.BORDER_REFLECT) mask = cv2.warpAffine(ones, M, (W, H), flags=interp, borderMode=cv2.BORDER_REFLECT) acc += warped wsum += mask # 除的是每个像素实际被覆盖的次数，不是帧数 return acc / np.maximum(wsum, 1e-6)关键在mask这一路。全 1 图经过同样的仿射变换和同样的插值核之后，每个像素上的值恰好等于该位置被这一帧覆盖的权重，累加起来就是总覆盖度。最后用acc / wsum而不是acc / len(frames)，边缘和中间一视同仁，都是各自的实际证据量。这也顺带解决了另一个问题：如果某一帧因为配准失败被剔除，帧数变了但公式不用改。
borderMode这里选BORDER_REFLECT而不是默认的BORDER_CONSTANT，配合 mask 一起用，边缘既不会拉黑，也不会因为补零导致权重虚高。代价是极边缘处会有一点镜像重复的纹理，但在 4 倍放大下基本看不出来，比黑边好看太多。
融合结果还有一层隐含假设：帧间亮度一致。如果连拍时自动曝光在动，某个区域在 A 帧偏亮、B 帧偏暗，累加平均后会出现块状的明暗断层。解决办法是在累加前对每帧做一个简单的均值-方差归一化，或者把每帧直方图匹配到参考帧。这一步不是必需的，但真实数据上建议开着。
3.3 迭代反投影：用正向退化模型把模糊「推」回去
shift-and-add 出来的结果已经比双三次好，但它本质是加权插值，高频仍然偏软。迭代反投影（IBP）的思路很直白：假定当前 HR 估计是对的，用正向退化模型把它降质回每一帧 LR，跟真实帧比，差多少就往回投影多少，迭代收敛。
def _forward(hr, dx, dy, scale, sigma=1.2): """正向退化：HR -> 平移 -> 模糊 -> 下采样，模拟第 i 帧是怎么来的""" M = np.array([[1, 0, -dx], [0, 1, -dy]], np.float32) shifted = cv2.warpAffine(hr, M, (hr.shape[1], hr.shape[0]), flags=cv2.INTER_LINEAR, borderMode=cv2.BORDER_REFLECT) # 高斯核在 HR 尺度上要按倍数放大，否则 LR 端的模糊量对不上 blur = cv2.GaussianBlur(shifted, (0, 0), sigma * scale / 2.0) h, w = hr.shape return cv2.resize(blur, (w // scale, h // scale), interpolation=cv2.INTER_AREA) def _back(err, dx, dy, scale): """把 LR 上的残差放大并反向平移，投回 HR 网格""" up = cv2.resize(err, None, fx=scale, fy=scale, interpolation=cv2.INTER_LINEAR) M = np.array([[1, 0, dx], [0, 1, dy]], np.float32) return cv2.warpAffine(up, M, (up.shape[1], up.shape[0]), flags=cv2.INTER_LINEAR, borderMode=cv2.BORDER_REFLECT) def ibp(hr_init, frames, shifts, scale=4, iters=20, lam=0.7, tol=1e-3): """lam 是步长，0.5~0.8 之间比较稳；再大边缘会振铃""" hr = hr_init.astype(np.float32).copy() prev = float("inf") for _ in range(iters): err_map = np.zeros_like(hr) for img, (dx, dy) in zip(frames, shifts): sim = _forward(hr, dx, dy, scale) err_map += _back(img - sim, dx, dy, scale) hr = np.clip(hr + lam * err_map / len(frames), 0, 255) cur = float(np.mean(np.abs(err_map))) if prev - cur < tol: # 残差不再下降就提前停，省掉一半时间 break prev = cur return hr三个地方值得单独说。第一，_forward里高斯核的 sigma 要乘scale / 2。因为退化模型是在 HR 尺度上做模糊再下采样，而实际成像的模糊发生在 LR 尺度，不换算的话正演出来的 LR 跟真实 LR 差着一圈光晕，残差全是系统偏差，迭代越跑越偏。这个换算关系是血泪经验，很多人第一次跑 IBP 觉得没效果，问题就出在这。
第二，lam别贪大。0.7 是个稳妥值，1.0 理论上是「一步到位」，但实际因为正演模型不精确，全步长会导致边缘过冲，出现一圈白边。想更快收敛，先把sigma调准，而不是把lam调大。
第三，提前停止很有必要。默认 20 轮，实际上多数序列 8～12 轮残差就压不动了，后面全是震荡。tol=1e-3这个阈值配合mean(abs(err_map))判断，简单够用。
3.4 评价与可视化：PSNR 算不对等于白跑
from skimage.metrics import peak_signal_noise_ratio as psnr from skimage.metrics import structural_similarity as ssim def evaluate(sr, gt): """sr / gt 都是 [H, W] float32，灰度，范围 0~255""" sr = sr.astype(np.float32) gt = gt.astype(np.float32) # 先做一次整数像素对齐：重构结果常带整体亚像素偏移， # 不校正直接算 PSNR 会低 2~3 dB，这个锅不该由算法背 (dx, dy), _ = cv2.phaseCorrelate(gt, sr) M = np.array([[1, 0, round(dx)], [0, 1, round(dy)]], np.float32) sr_aligned = cv2.warpAffine(sr, M, (sr.shape[1], sr.shape[0]), flags=cv2.INTER_LINEAR, borderMode=cv2.BORDER_REFLECT) return { "psnr": float(psnr(gt, sr_aligned, data_range=255)), "ssim": float(ssim(gt, sr_aligned, data_range=255)), }指标这一步最容易自欺欺人。重构结果经常整体偏移零点几个像素，直接拿原图和它算 PSNR，会莫名其妙低两三个 dB，然后你以为算法不行，其实是评价方式不对。所以先估一次位移、取整、对齐，再算。同理，SSIM 对亮度偏移不敏感但对结构错位极其敏感，两个指标一起看更靠谱。
可视化我一般拼成三行四列：第一行双三次放大、shift-and-add、IBP 结果；第二行对应的误差热力图（跟真值绝对差放大 4 倍）；第三行放局部放大裁剪。热力图比 PSNR 更能说明问题——IBP 到底是在恢复边缘，还是只是在平坦区压噪声，一眼就看出来。
4. 跑通与复现：从合成数据到真实连拍的两条入口
4.1 环境准备
先把依赖装干净。不要用系统 Python，也不要混着装 conda 和 pip，版本冲突导致cv2导入失败是新手最常见的翻车点。
# 建议 Python 3.9+，纯 CPU 即可 python -m venv .venv source .venv/bin/activate # Windows 用 .venv\Scripts\activate pip install numpy opencv-python scikit-image pyyaml matplotlib # 版本不求新，opencv-python 4.x 任意一个近期版本都行 python -c "import cv2, numpy, skimage; print(cv2.__version__)"最后那行验证命令别省。它同时检查 cv2、numpy、skimage 三个包能不能正常导入，装完就报错比跑到一半报错好定位得多。见过装的是opencv-python-headless和opencv-python混用，导致imshow找不到的情况，那属于依赖污染，重装虚拟环境最快。
4.2 一键入口：run_demo.py 怎么跑
合成数据入口是这条管线的「绝对能运行」保证——它不依赖任何外部素材，自己造真值、自己退化、自己重构、自己算指标。
# run_demo.py 核心逻辑（省略参数解析） import cv2, numpy as np, yaml from mfsr.degrade import make_lr_frames from mfsr.register import estimate_shift from mfsr.fuse import shift_and_add from mfsr.backproject import ibp from mfsr.metrics import evaluate cfg = yaml.safe_load(open("configs/demo.yaml")) scale = cfg["scale"] # 1) 载入真值（没有就跑 make_gt.py 生成一张测试图） gt = cv2.imread("data/synthetic/gt.png", cv2.IMREAD_GRAYSCALE).astype(np.float32) # 2) 造 8 帧，位移是高分辨率像素单位，故意做成非整数 shifts_hr = [(-1.5, 0.7), (0.8, -1.2), (1.9, 1.4), (-0.6, -1.8), (2.4, -0.3), (-2.1, 0.9), (0.3, 2.2), (-1.1, -0.5)] frames = make_lr_frames(gt, shifts_hr, scale=scale, sigma=cfg["sigma"], noise=cfg["noise"], seed=7) # 3) 配准：LR 位移 * scale 才是 HR 网格上的位移 shifts_lr = [(img, *estimate_shift(frames[0], img)[:2]) for img in frames[1:]] shifts_hr_est = [(dx * scale, dy * scale) for _, dx, dy in shifts_lr] # 4) 融合 + 反投影 sr0 = shift_and_add(frames, shifts_hr_est, scale=scale) sr1 = ibp(sr0, frames, shifts_hr_est, scale=scale, iters=cfg["iters"], lam=cfg["lam"]) # 5) 三种结果对比 bicubic = cv2.resize(frames[0], None, fx=scale, fy=scale, interpolation=cv2.INTER_CUBIC) for name, img in [("bicubic", bicubic), ("saa", sr0), ("ibp", sr1)]: print(name, evaluate(img, gt))make_lr_frames内部的退化顺序是：亚像素平移（INTER_CUBIC插值）→ 高斯模糊 →INTER_AREA下采样 → 加高斯噪声。这个顺序和真实成像接近，也是做超分实验的通用做法。参数sigma控制模糊强度，noise控制噪声标准差，两个值决定了任务难度——sigma=0时 IBP 提升会非常明显，sigma调到 2.0 以上，传统方法基本就只能跟双三次打平。
跑完你会看到三行指标，正常区间是：双三次 PSNR 24～26 dB，shift-and-add 26～28 dB，IBP 27～30 dB。如果你的 IBP 比双三次还低，先查位移符号，再看sigma有没有在下采样前生效。
4.3 换成自己的数据
真实序列的流程一样，只是位移靠估。落地时按这个顺序走：
第一步，把连拍序列重命名成frame_000.png、frame_001.png……放data/real/。选帧有讲究，优先挑时间上分散的，相邻两帧位移太小（都小于 0.2 个 LR 像素）等于没有新增采样。
第二步，看每帧的配准 response，把低于 0.2 的帧剔掉。这一步比调任何参数都管用。
第三步，跑run_real.py，它会用estimate_shift现场估位移，做 shift-and-add 和 IBP，输出拼图和中间结果。没有gt.png就不算 PSNR，改看局部裁剪的视觉效果和热力图。
第四步，调整iters。真实数据一般不需要 20 轮，噪声大的序列跑多了反而放大噪点，我通常在 8～12 轮之间停手。
提示：真实序列不要指望一次到位。先只跑前 4 帧看融合是否清晰，清晰了再加帧，这是最省时间的排查顺序。
5. 避坑与排查：配准、融合、评价三段路上最容易翻车的五个点
5.1 配准段：符号、单位、亮度
现象：融合结果比双三次还糊，而且帧数越多越糊。原因：位移符号反了，每帧都是往反方向平移后再平均，等于主动做了运动模糊。 解决：造一组已知位移为(1.5, 0.5)的合成帧，调用estimate_shift，看返回的dx是正还是负。如果符号相反，在shift_and_add的M里把dx改成-dx，或者直接对所有估计值取负。这个验证必须做，不能靠猜。
现象：配准看起来对，但融合图整体往一个方向偏了小半个像素，指标一直上不去。原因：单位搞混了。estimate_shift在 LR 图上估出的位移是 LR 像素单位，直接喂给 HR 网格的warpAffine就少了scale倍的放大。 解决：在配准模块的出口统一乘scale，对外只暴露 HR 单位。这条约定写进函数 docstring，比事后调试便宜。
现象：某几帧融进去之后，画面出现块状明暗断层。原因：连拍时自动曝光在动，帧间存在增益差异，加权平均后局部亮度被拉花。 解决：在shift_and_add累加前对每帧做均值-方差归一化，或者把每帧直方图匹配到参考帧。归一化只改亮度和对比度，不改结构，对配准和融合都无害。
5.2 融合与评价段：空网格、过冲、指标错位
现象：输出图边缘有一圈黑边或者放射状的脏条。原因：warpAffine用了BORDER_CONSTANT补零，边缘像素被当成真实数据累加，除以帧数后仍然偏暗。 解决：改成BORDER_REFLECT，并且必须用wsum而不是帧数做分母。如果只改边界模式不改除法方式，边缘会被过度加权，一样出问题。
现象：IBP 迭代几轮之后，边缘出现一圈白边，PSNR 反而下降。原因：lam过大或iters过多，正演模型不精确时残差会围绕边缘震荡，形成过冲。 解决：把lam压到 0.5～0.7，加上残差提前停止条件；仍然过冲的话，在每轮之后加一次轻微的中值滤波或双边滤波，代价是细节损失一点，但稳定性提升明显。
另外补一条容易被忽略的：指标错位。重构结果和真值之间存在整体亚像素偏移时，PSNR 会平白掉 2～3 dB，让人误判算法失效。算指标前先估一次位移、取整、对齐，再算。这一步在evaluate里已经内建，别自己另写一套。
6. 二次开发：把固定管线改成可插拔实验台
这套代码最值得动手改的地方，是让配准和融合变成接口，而不是写死的函数调用。我在自己的实验里把register.py收敛成一个签名(ref, mov) -> (dx, dy, score)，然后并排放了三个实现：相位相关、cv2.findTransformECC的MOTION_EUCLIDEAN（能估旋转）、以及一个把帧上采样后做块匹配的土办法。换配准只要改配置里的一行字符串，跑同一批合成数据，谁高谁低一目了然。这个改动花不到一小时，换来的对比实验能力是长期的。
验证方法上，有三个对照实验必须跑。第一组固定sigma=0，看管线在无模糊条件下的上限，这一组的 PSNR 应该接近理想值，如果不行说明融合或反投影有 bug。第二组扫sigma从 0.5 到 2.2，画出 PSNR 随模糊强度下降的曲线，这条曲线就是这套方法的适用边界。第三组固定其他参数扫帧数 2/4/8/12，你会看到一个很典型的收益递减——通常到 6～8 帧后，配准误差累积的速度超过了新增采样的收益，帧数再多纯属浪费。这个拐点值不值得跑，直接决定了你在真实采集端要不要连拍二十张。
如果你想接深度学习方案，最省事的接法是保留这份代码的数据层和评价层，只把fuse.py换成一个轻量网络：以shift_and_add的结果加上各帧作为输入，真值用合成数据自动生成，训练目标就是 L1。好处是配准仍然可解释、可调试，网络只负责最后那点高频恢复，训崩了你随时能退回传统基线。
| 改造点 | 改动量 | 收益 |
|---|---|---|
| 配准接口化 | 小 | 三类配准横向对比，定位瓶颈 |
| 融合替换为学习模块 | 中 | 高频恢复提升，基线仍可回退 |
| 加 TV 正则进 IBP | 小 | 抑制噪声放大，代价是细节略软 |
| 多通道并行 | 小 | 彩色图端到端，但耗时翻三倍 |
从那以后我每次接一个新的重构管线，都强制先用合成数据跑一遍符号验证和sigma=0的上限测试，再碰真实序列。这半小时看起来是浪费，实际上能挡掉后面一整天的自我怀疑。希望帮到你。
本文还有配套的精品资源，点击获取
