简介这套FSIM特征相似性计算代码为图像质量评估和图片相似性对比提供了轻量级参考实现。与PSNR、SSIM等传统指标相比FSIM通过相位一致性与梯度幅值刻画结构特征能更细腻地反映视觉差异适合图像处理、计算机视觉方向的学生、研究者或算法工程师开展算法效果验证与特征相似性实验。压缩包仅2.21MB共5个文件包含Python版FSIM实现、配套测试脚本、Matlab版FeatureSIM函数以及两张BMP测试图片用户可根据实际工程环境选用对应源码并通过自带测试样本快速理解调用方式。资源下载页已有1400人学习/下载代码结构简洁、可直接运行输出相似度数值两份源码既可独立学习比对也能作为基础模块嵌入到图像质量评价、相似图片检索或模型效果评估任务中有助于节省从零编写FSIM算法的时间。1. 从 PSNR 到 FSIM为什么图片质量评价要看特征做图片质量评价的人大概率都遇过这种尴尬一张 JPEG 从质量 90 压到 85PSNR 只掉了 0.8dBSSIM 还在 0.99肉眼却已经能看出边缘发糊。换 FSIM 计算代码跑一遍分数从 0.983 掉到 0.927差异一下就拉开了。FSIMFeature Similarity特征相似性是图片质量评价领域公认比 PSNR、SSIM 更贴近人眼判断的参数也是图片相似性评价时经常被拉出来对比的指标。检索时经常能看到“结构相似性 FSIM 算法”这种写法这里顺带说清楚FSIM 与 SSIM 一样属于全参考评价但一个看结构函数分解一个看特征响应一致性后者在纹理重排、局部相位变化这类场景下要敏感得多。这份资源包含 MATLAB 版 FeatureSIM.m、Python 版 FSIM.py、测试脚本 test.py以及两张配套测试图 0000.bmp 和 0001.bmp适合做超分、去噪、压缩编码效果对比的同学也适合需要把“差异不明显”这类主观描述变成可量化分数的场景。2. FSIM 算法拆解相位一致性与梯度幅度如何构成相似性映射2.1 PSNR 与 SSIM 为什么看不清纹理差异先明确一个前提PSNR 是像素级差的 dB 化表达整张图所有位置的误差一视同仁SSIM 把亮度、对比度、结构三项相乘本质上做的是局部块统计。两个指标对加性噪声、轻微模糊都敏感但一旦失真表现为纹理重排、边缘相位偏移它们的反应会远小于人眼感知。这一点在做压缩算法对比时尤其明显块效应被平滑掉之后 SSIM 可能反而升高主观质量却在下降。FSIM 的出发点是相位一致性Phase CongruencyPC理论人眼感知到的特征点通常出现在傅里叶分量相位最一致的位置。边缘、角点这类结构化信息并不依赖特定尺度或方向而是由相位关系决定。FSIM 用 PC 作为底层特征图再用梯度幅度做补充两张图的相似性就建立在“特征响应是否一致”上而不是“某个区域像素是否接近”。所以拆这份 FSIM 源码时重点看三块相位一致性怎么算、梯度幅度怎么算、两个相似性映射怎么合成。下面先从 FeatureSIM.m 里最费解的相位一致性部分开始讲这部分也是 MATLAB 和 Python 版本移植时最容易出偏差的地方。2.2 相位一致性 PC 的计算路径与关键参数相位一致性的标准计算分三步构造多尺度多方向的 log-Gabor 滤波器组对图像做频域滤波得到复数响应把相邻方向的响应能量组合成 PC 值。Kovesi 的经典实现会在响应幅度上做噪声补偿FeatureSIM.m 里同样保留了这套逻辑。常见参数是 4 个尺度、4 个方向尺度用 2 的指数递增方向取 0、π/4、π/2、3π/4覆盖从细边缘到粗轮廓的信息。% FeatureSIM.m 相位一致性主循环的常见组织方式 T1 0.85; % 相位相似性小常量防止除零也限定灵敏度 baseWavelength 3; % 最细尺度对应的波长像素 sigmaOnf 0.75; % log-Gabor 径向带宽控制 numScale 4; % 尺度数量 numAngle 4; % 方向数量 [rows, cols] size(img1); [X, Y] meshgrid(1:cols, 1:rows); radius sqrt((X - cols/2).^2 (Y - rows/2).^2); radius(rows/2 1, cols/2 1) 1; % 频域中心置 1避免 log(0) angle atan2(Y - rows/2, X - cols/2); for s 1:numScale wavelength baseWavelength * 2^(s - 1); % 尺度按 2 的幂递增 logGabor exp(-(log(radius ./ wavelength)).^2 / (2 * sigmaOnf^2)); for n 1:numAngle theta (n - 1) * pi / numAngle; % 当前方向角度 spread exp(-(min(abs(angle - theta), pi - abs(angle - theta))).^2 / (2 * 0.45^2)); filter logGabor .* spread; filter(rows/2 1, cols/2 1) 0; % 去掉直流分量 % 频域乘完做逆变换得到该尺度方向的复数响应再参与能量组合 end endwavelength 从 3 像素递增到 24 像素覆盖从细边缘到粗轮廓的信息sigmaOnf 越大径向带宽越宽、对频率偏移越不敏感0.75 是多数实现里手感较好的值。方向项用min(abs(angle-theta), pi-abs(angle-theta))处理角度周期避免 0 和 π 被切开。每个滤波器最后把直流分量置零等效于减去图像均值这样光照整体平移不会影响 PC 结果。四个方向的响应会按相位差组合出主能量响应 PC1 和垂直向响应 PC2。参数常见取值作用调整方向baseWavelength3最细尺度波长图像纹理更细时调小sigmaOnf0.75log-Gabor 径向带宽越大对频率偏差越不敏感numScale4尺度数量特征尺度跨度大时增加到 5numAngle4方向数量关注斜向纹理时可加到 62.3 从 PC 与 GM 到 FSIM 分数PC 给出的是特征强度梯度幅度Gradient MagnitudeGM给出的是边缘锐度。FeatureSIM.m 里 GM 用 Scharr 算子计算这一步计算量小但作用很大PC 对平坦区域响应低GM 能补充对比度变化的信息。两张图各自得到 PC1、PC2、GM1、GM2 之后先分别算相似性映射再做一次加权平均得到最终分数。% 相位一致性映射 S_PC (2 * PC1 .* PC2 T1) ./ (PC1.^2 PC2.^2 T1); % 梯度幅度映射 S_G (2 * GM1 .* GM2 T2) ./ (GM1.^2 GM2.^2 T2); % 特征相似性映射默认 alpha beta 1这里直接相乘 S_L S_PC .* S_G; % 用 PC 最大值加权权重更高处更可能是稳定特征 PCm max(PC1, PC2) eps; FSIM sum(sum(S_L .* PCm)) / sum(sum(PCm));两个相似性映射都是比值形式分子是交叉项分母是能量项T1、T2 的作用是让分母为 0 时结果可控。T1 取 0.85、T2 取 160 是相对稳定的共识值但 T2 和图像灰度范围强相关如果把图归一化到 [0,1]T2 要缩到 1.6 左右否则梯度项几乎不起作用。FSIM 最终结果落在 (0,1] 区间同一张图得分恒为 1差异越大分数越低。FeatureSIM.m 里还有一步对 PC 响应的归一化Python 移植时最容易漏掉的就是这里漏掉后分数区分度会变差但单张图上不一定看得出来。3. Python 实现对照从 MATLAB 移植到 numpy 的关键差异3.1 依赖与输入归一化FSIM.py 的常见实现是 numpy opencv-python numpy.fftOpenCV 负责读图和 Scharr 梯度FFT 负责频域滤波。移植时第一步不是写滤波器而是统一输入MATLAB 里 imread 读进来是 uint8FeatureSIM.m 会转 doublePython 如果直接用 cv2.imread 拿 uint8 做运算log-Gabor 响应和 T2 的匹配关系会全部错位。我一般会把图像直接归一化到 [0,1] 的 float64这样调试时两个语言版本输出更容易对齐。import numpy as np import cv2 from numpy.fft import fft2, ifft2, fftfreq def read_gray(path): img cv2.imread(path, cv2.IMREAD_GRAYSCALE) if img is None: raise FileNotFoundError(fcannot read image: {path}) return img.astype(np.float64) / 255.0 # 转到 [0,1]保证 T2 语义一致 def gradient_magnitude(img): gx cv2.Scharr(img, cv2.CV_64F, 1, 0) gy cv2.Scharr(img, cv2.CV_64F, 0, 1) return np.sqrt(gx * gx gy * gy)read_gray 里直接用 IMREAD_GRAYSCALE 读灰度彩色图会在读盘阶段做加权转换比在模块内转更省事也不容易出错。归一化到 [0,1] 之后T2 建议改成 1.6 而不是 160这一点在源码注释里通常会写清楚如果你看到 Python 版输出普遍偏低优先查 T2 和输入范围是否匹配。gradient_magnitude 用两级 Scharr 核OpenCV 的 CV_64F 保证梯度不截断这一行是 uint8 溢出高发区不要偷懒用默认取值。3.2 在 numpy 里构造 log-Gabor 滤波器组numpy 里构造滤波器组和 MATLAB 最大的不同是 FFT 的象限组织numpy.fft.fftfreq 返回的频率从 0 到 0.5 再到负半轴而 MATLAB 的 meshgrid 坐标是从 1 到 N、中心点在中间。如果照抄 MATLAB 的坐标公式方向响应会整体偏移 90 度。常见做法是按 fftfreq 生成半径和角度再逐尺度逐方向构造最后同样置零直流分量。def log_gabor_bank(shape, scales4, angles4, base_wavelength3.0, sigma_onf0.75): rows, cols shape fx fftfreq(cols) # 列方向频率单位是周期/像素 fy fftfreq(rows) FX, FY np.meshgrid(fx, fy) radius np.sqrt(FX * FX FY * FY) radius[0, 0] 1.0 # 直流位置置 1避免 log(0) angle np.arctan2(FY, FX) bank [] for s in range(scales): wavelength base_wavelength * (2 ** s) radial np.exp(-(np.log(radius / (1.0 / wavelength))) ** 2 / (2 * sigma_onf ** 2)) for n in range(angles): theta n * np.pi / angles diff np.abs(np.mod(angle - theta, np.pi)) diff np.minimum(diff, np.pi - diff) angular np.exp(-(diff ** 2) / (2 * 0.45 ** 2)) g radial * angular g[0, 0] 0.0 # 去掉均值分量等价于对图像去中心化 bank.append(g) return bank注意1.0 / wavelength是频率中心log-Gabor 在频域用对数频率偏移定义径向响应所以看起来和 MATLAB 版本不同本质相同。wavelength 从 3 递增到 24对应频率中心从 0.333 降到 0.042。角度项做了 π 周期折叠因为方向滤波器不需要区分正反方向。bank 里 16 个滤波器按方向分组每组 4 个尺度。特征更碎的图像可以把 base_wavelength 调到 2但要注意最粗尺度也同步缩小覆盖频率带整体向高频移动对噪声会更敏感。3.3 相位一致性响应与 FSIM 汇聚有了滤波器组相位一致性的计算就是把图像 FFT 后和每个滤波器相乘做逆变换得到复数响应再按方向把尺度的能量组合起来。实际项目里我不会在这个函数里做太多花活基础相位一致性、噪声补偿、响应归一化三件事分离方便不同语言版本对照。下面是简化但功能完整的核心逻辑。def phase_congruency(img, bank, scales4, angles4): rows, cols img.shape img_f fft2(img) # 按方向组织响应每个方向有 scales 个尺度的复数响应 responses np.empty((angles, scales, rows, cols), dtypenp.complex128) for n in range(angles): for s in range(scales): f bank[n * scales s] responses[n, s] ifft2(img_f * f) # 相位一致性近似取最大方向能量与总能量之比 energy np.abs(responses) ** 2 energy_sum energy.sum(axis1) # 沿尺度累加 pc energy_sum.max(axis0) / (energy_sum.sum(axis0) 1e-10) return pc这里做了一个工程化简化经典算法会用正交滤波器对做相位一致性而不是直接能量比两者的主趋势一致数值上会有差异。如果你要和 FeatureSIM.m 的输出逐像素对齐需要按 Kovesi 原版公式补上噪声估计和响应归一化。代码里用 np.empty 分配复数数组避免 list append 的额外开销大图下这个函数是 FSIM 的主要耗时点瓶颈在 16 次 ifft2可以先用小图验证逻辑再考虑用 scipy.fft 的多线程参数。最终汇聚和 MATLAB 完全对应S_PC、S_G 求出后相乘得到 S_L再用 PCm 加权平均。唯一要注意的是 epsMATLAB 的 eps 是双精度浮点最小间隔Python 里用 1e-10 代替即可太小会导致除零警告。def fsim_score(img1, img2, bank, T10.85, T21.6): if img1.shape ! img2.shape: raise ValueError(fshape mismatch: {img1.shape} vs {img2.shape}) pc1 phase_congruency(img1, bank, scales4, angles4) pc2 phase_congruency(img2, bank, scales4, angles4) gm1 gradient_magnitude(img1) gm2 gradient_magnitude(img2) pc_m np.maximum(pc1, pc2) 1e-10 # 加权权重 s_pc (2 * pc1 * pc2 T1) / (pc1 ** 2 pc2 ** 2 T1) s_g (2 * gm1 * gm2 T2) / (gm1 ** 2 gm2 ** 2 T2) s_l s_pc * s_g return np.sum(s_l * pc_m) / np.sum(pc_m)函数开头先校验 shape这是批量调用时最容易翻车的地方早抛出比后面输出全是 NaN 好查。T2 默认 1.6 是配合 [0,1] 归一化的取值如果读图脚本去掉了归一化记得改回 160。整个函数保持纯函数风格bank 从外部传入批量打分时只需要构造一次滤波器组每张图复用即可。3.4 MATLAB 与 Python 实现的边界差异两版代码在同一对测试图上跑分数应在小数点后两位一致。如果差得多按下面的表逐一排查90% 的情况出在前三行。差异点MATLAB 行为Python 行为影响索引起点下标从 1 开始ndarray 从 0 开始滤波器中心错一位方向响应会翻转图像类型uint8 转 double 不归一化常归一化到 [0,1]T2 需要同步缩放卷积边界imfilter 默认补零cv2.Scharr 默认边界反射边缘 10px 内差异放大直流分量手动置零需要同样手动处理不置零则亮度偏移影响 PCeps约为 2.2e-16常用 1e-10平坦区域取值影响明显另一个隐蔽差异是 meshgrid 的参数顺序。MATLAB 的 meshgrid(x, y) 生成的是按行变化的 Y、按列变化的 Xnumpy 的 meshgrid 默认同样是笛卡尔序但如果换用 np.mgrid 就会顺手写反。方向滤波器对角度敏感行列反了之后斜向纹理的 PC 响应会互相错位整体分数仍然有区分度但和 MATLAB 版对不上。提示做移植对照时不要只对比最终分数把中间变量 S_PC、S_G 各存一份做逐元素比较。分数一样不代表特征图一致特征图一致才说明滤波器相位方向没搞反。4. 用 test.py 跑通一次 FSIM 评价流程4.1 test.py 在测什么包里的 test.py 是给第一次接触 FSIM 的人准备的入口通常做三件事读入 0000.bmp、0001.bmp调用 FSIM 计算函数打印分数。这两张图一张是参考图一张是经过处理的对比图分数能直接反映两张图的相似程度。把 test.py 里的文件路径换掉就可以把测试脚本变成你自己的对比工具也就是说你不需要重新组织代码只改路径就能验证自己的图片。实现上通常就是读图、转灰度、调 fsim_score 三步有些版本会把两张图的预览图也 show 出来方便你确认输入没有读错。4.2 MATLAB 端运行与输出如果你主要用 MATLAB把 FeatureSIM.m 和两张 bmp 放进同一目录在命令窗口直接调函数即可。FeatureSIM.m 的入参是两张图像矩阵不是文件路径所以要先 imread。score FeatureSIM(imread(0000.bmp), imread(0001.bmp)); fprintf(FSIM %.4f\n, score);imread 读入的是 uint8 三通道数据FeatureSIM.m 内部会先转灰度并转 double所以外部不需要预处理。输出是一个 0 到 1 之间的标量建议在 fprintf 里保留四位小数两位小数在 0.99 这个区间看不出差异。如果报错提示矩阵维度不一致先检查两张图的宽高是否相同FSIM 不做缩放尺寸不同必须预处理。4.3 Python 端运行与报错排查Python 版入口有两种调用方式直接跑 test.py或把 FSIM.py 当模块导入。前者适合验证环境后者适合写批量脚本。无论哪种文件路径和 read_gray 的路径处理保持一致避免中文路径下 OpenCV 的读取问题。python test.py python FSIM.py --ref 0000.bmp --dist 0001.bmp常见报错就三种。第一是 cv2.imread 返回 None多半是路径写错或权限问题read_gray 里抛 FileNotFoundError 比静默崩溃好查。第二是 shape mismatch参考图和失真图分辨率不同多数数据集里成对图片不会刚好一致建议在脚本里强制 resize 成同一尺寸。第三是输出 NaN几乎都发生在全黑或全白区域PC 分母为 0检查 phase_congruency 里有没有加 1e-10 这个常数。提示如果你看到 FSIM 大于 1.0说明归一化或 eps 的取值有问题。FSIM 的上界是 1同图输出大于 1.0001 就不要再往下分析了先查输入图像矩阵的范围。4.4 分数边界参考拿到一个分数之后怎么判断它合理我一般用下面这个经验范围做 sanity check。注意这不是固定阈值具体项目和图像内容会漂但数量级不会错。对比情况FSIM 参考范围说明同一张图1.0恒等用于自检算法实现同内容不同压缩率0.97 ~ 0.99肉眼可感知但差异不大明显失真模糊/噪声0.90 ~ 0.96主观评价多数人会说“有差距”完全不同内容0.60 ~ 0.85内容无关仅作下界参考如果两张完全不同内容的图跑到 0.9 以上基本可以判定实现出错了。最常见的原因是两张图被误当成同一张图比较或者相位一致性计算时把 bank 里的 16 个滤波器索引循环写成了同一个。用 test.py 先做一次同图自检再做一次明显差异图对比能快速确认链路完整。5. 把 FSIM 用到批量数据集上目录遍历、评分归一化与异常检测5.1 批量打分脚本FSIM 真正发挥价值的地方是批量对比把算法的输出和参考图逐对打分。常见做法是一个目录放参考集另一个目录放失真集文件名字相同。脚本里只构造一次滤波器组避免每张图都重建 16 个频域滤波器这是最容易忽略的性能瓶颈图一多能差出几倍时间。import os, glob import numpy as np def batch_fsim(ref_dir, dist_dir, suffix.bmp): ref_files sorted(glob.glob(os.path.join(ref_dir, * suffix))) if not ref_files: return {} H, W read_gray(ref_files[0]).shape # 用第一张图确定滤波器尺寸 bank log_gabor_bank((H, W)) scores {} for ref_path in ref_files: name os.path.basename(ref_path) dist_path os.path.join(dist_dir, name) if not os.path.exists(dist_path): continue # 缺图的样本先跳过计数留到后面统计 img1 read_gray(ref_path) img2 read_gray(dist_path) scores[name] fsim_score(img1, img2, bank) return scoresH、W 先由第一张参考图决定批量前用一次 read_gray 获取 shape避免每次循环里重复取。缺图的样本先跳过而不是抛异常跑完统计缺失数即可如果脚本因为一张图中断前面几个小时的耗时全白费。sorted 保证输出顺序稳定后面画曲线或者写 CSV 时不容易对错行。5.2 归一化与阈值选取批量分数出来后不同内容的参考图本身 FSIM 基数就不同直接拿 0.95 固定阈值去卡会误伤。常见做法是对每组同内容的分数做 min-max 归一化再定阈值或者干脆只看排序不看绝对值。比如 100 张测试图第 95 张之后的分数明显掉一截这个拐点就是参考阈值。def normalize_scores(scores): arr np.array(list(scores.values())) lo, hi arr.min(), arr.max() return {k: (v - lo) / (hi - lo) for k, v in scores.items()}min-max 归一化要按失真类型分开做不要把所有压缩、噪声、模糊混在一起归一化。混在一起会让阈值被最差的一组拉走局部差异反而看不出来。5.3 三个容易踩的坑批量场景下三个坑最常出现。第一参考图和失真图通道数不一致一边是 PNG 灰度、一边是 JPEG 彩色read_gray 读出来都是灰度但彩色转灰度和直接读灰度的权重不同分数会系统性偏低。统一从 RGB 读再转灰度最保险。第二尺寸不一致。FSIM 不做多尺度融合输入分辨率不同 PC 响应本来就不对齐建议统一 resize 到同一个尺寸再比。不要用 cv2.IMREAD_UNCHANGED 跳过通道处理彩色 alpha 通道进到 float 计算里会把梯度带偏。第三平坦区域 NaN。全黑、全白、大面积纯色背景在 phase_congruency 里能量极低加权时 PCm 会趋向 0S_L 却可能因为除零放大。处理方式是在最终加权前做一次掩码把 PCm 低于 1e-6 的位置直接从中剔除而不是只靠 eps 硬撑。放一张纯色图在测试集里跑一遍很快就能验证这层防抖逻辑有没有生效建议把它写进批量脚本的断言里。本文还有配套的精品资源点击获取