FSIM图像质量评价：从相位一致性到Python实现

简介：这套FSIM（特征相似性）计算代码为图像质量评估和图片相似性对比提供了轻量级参考实现。与PSNR、SSIM等传统指标相比，FSIM通过相位一致性与梯度幅值刻画结构特征，能更细腻地反映视觉差异，适合图像处理、计算机视觉方向的学生、研究者或算法工程师开展算法效果验证与特征相似性实验。压缩包仅2.21MB，共5个文件，包含Python版FSIM实现、配套测试脚本、Matlab版FeatureSIM函数以及两张BMP测试图片，用户可根据实际工程环境选用对应源码，并通过自带测试样本快速理解调用方式。资源下载页已有1400人学习/下载，代码结构简洁、可直接运行输出相似度数值；两份源码既可独立学习比对，也能作为基础模块嵌入到图像质量评价、相似图片检索或模型效果评估任务中，有助于节省从零编写FSIM算法的时间。

1. 从 PSNR 到 FSIM：为什么图片质量评价要看特征

做图片质量评价的人，大概率都遇过这种尴尬：一张 JPEG 从质量 90 压到 85，PSNR 只掉了 0.8dB，SSIM 还在 0.99，肉眼却已经能看出边缘发糊。换 FSIM 计算代码跑一遍，分数从 0.983 掉到 0.927，差异一下就拉开了。FSIM（Feature 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 Congruency，PC）理论：人眼感知到的特征点，通常出现在傅里叶分量相位最一致的位置。边缘、角点这类结构化信息并不依赖特定尺度或方向，而是由相位关系决定。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 end

wavelength 从 3 像素递增到 24 像素，覆盖从细边缘到粗轮廓的信息；sigmaOnf 越大径向带宽越宽、对频率偏移越不敏感，0.75 是多数实现里手感较好的值。方向项用min(abs(angle-theta), pi-abs(angle-theta))处理角度周期，避免 0 和 π 被切开。每个滤波器最后把直流分量置零，等效于减去图像均值，这样光照整体平移不会影响 PC 结果。四个方向的响应会按相位差组合出主能量响应 PC1 和垂直向响应 PC2。

2.3 从 PC 与 GM 到 FSIM 分数

PC 给出的是特征强度，梯度幅度（Gradient Magnitude，GM）给出的是边缘锐度。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.fft，OpenCV 负责读图和 Scharr 梯度，FFT 负责频域滤波。移植时第一步不是写滤波器，而是统一输入：MATLAB 里 imread 读进来是 uint8，FeatureSIM.m 会转 double；Python 如果直接用 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(f"cannot 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, scales=4, angles=4, base_wavelength=3.0, sigma_onf=0.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, scales=4, angles=4): rows, cols = img.shape img_f = fft2(img) # 按方向组织响应，每个方向有 scales 个尺度的复数响应 responses = np.empty((angles, scales, rows, cols), dtype=np.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(axis=1) # 沿尺度累加 pc = energy_sum.max(axis=0) / (energy_sum.sum(axis=0) + 1e-10) return pc

这里做了一个工程化简化：经典算法会用正交滤波器对做相位一致性，而不是直接能量比，两者的主趋势一致，数值上会有差异。如果你要和 FeatureSIM.m 的输出逐像素对齐，需要按 Kovesi 原版公式补上噪声估计和响应归一化。代码里用 np.empty 分配复数数组，避免 list append 的额外开销；大图下这个函数是 FSIM 的主要耗时点，瓶颈在 16 次 ifft2，可以先用小图验证逻辑，再考虑用 scipy.fft 的多线程参数。

最终汇聚和 MATLAB 完全对应：S_PC、S_G 求出后相乘得到 S_L，再用 PCm 加权平均。唯一要注意的是 eps：MATLAB 的 eps 是双精度浮点最小间隔，Python 里用 1e-10 代替即可，太小会导致除零警告。

def fsim_score(img1, img2, bank, T1=0.85, T2=1.6): if img1.shape != img2.shape: raise ValueError(f"shape mismatch: {img1.shape} vs {img2.shape}") pc1 = phase_congruency(img1, bank, scales=4, angles=4) pc2 = phase_congruency(img2, bank, scales=4, angles=4) 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% 的情况出在前三行。

另一个隐蔽差异是 meshgrid 的参数顺序。MATLAB 的 meshgrid(x, y) 生成的是按行变化的 Y、按列变化的 X；numpy 的 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。注意这不是固定阈值，具体项目和图像内容会漂，但数量级不会错。

如果两张完全不同内容的图跑到 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 scores

H、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 会趋向 0，S_L 却可能因为除零放大。处理方式是在最终加权前做一次掩码，把 PCm 低于 1e-6 的位置直接从中剔除，而不是只靠 eps 硬撑。放一张纯色图在测试集里跑一遍，很快就能验证这层防抖逻辑有没有生效，建议把它写进批量脚本的断言里。

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