简介面向偏振图像处理与多源图像融合的 MATLAB 代码包适合光学成像、遥感、医学影像等方向的研究人员和工程师参考。压缩包内共1个文件即 qzw3.m 脚本包体仅2KB结构紧凑便于快速阅读与复用。脚本覆盖偏振度计算、偏振相角提取以及偏振信息与原始强度图像融合的主要流程可帮助理解去反射、增强对比度、揭示隐藏结构等实际效果并观察融合后强度图像的信息增益。已有470人学习下载说明该资料在这一细分主题下具备一定参考热度。借助脚本读者既能对照算法步骤学习偏振物理量的计算方式也能将融合结果用于后续分析或作为原型嵌入自身项目。对于希望快速入门偏振融合、验证相关算法或搭建最小示例的开发者这份轻量资源提供了直观且可扩展的参考实现。1. 偏振融合不是新瓶装旧酒一份 zip 里藏着从偏振度到图像融合的完整链路拿到qzw3.zip这个包时我以为又是一堆随手整理的算法残片。解压之后翻完目录才发现它把偏振度DoP、偏振强度、斯托克斯参数计算和图像融合串成了一条能直接跑的完整链路。做工业视觉这几年普通强度图在镜面、透明件、低对比度目标面前经常翻车偏振信息恰好补上了这一块缺口。这份资源适合正在做图像融合、偏振成像或者工业检测的从业者也适合刚入门、想找一份能跑通的参考实现的人。它的价值不在于某个单点算法多高深而在于从四通道偏振原图到融合结果之间的每一环都给了可验证的代码和参数。2. 偏振数据的三件套DoP、AOP、强度图怎么算才不翻车2.1 偏振度是图像融合里的“信息增量”普通摄像头拿到的是强度图丢掉了光的偏振状态。偏振度Degree of Polarization, DoP描述的是光在某个方向上的偏振分量占总强度的比例它反映的是物体表面的反射特性、材质和粗糙度。同一块场景里金属、玻璃、水面、塑料的DoP值差异往往比强度图里的灰度差异大得多这就是融合时真正的信息增量。偏振方向也有价值。偏振角Angle of Polarization, AOP与物体表面的法线方向、入射角相关做形状恢复或者高光抑制时会用到。把强度图、DoP图、AOP图三者融合等于同时用了亮度、材质、几何三类信息比单纯加权叠加要稳得多。这份资源里默认走的就是这条路线先把四方向偏振原图转成 Stokes 参数再算出 DoP、AOP、强度图最后做多尺度融合。2.2 Stokes 参数是偏振计算的通用接口偏振相机输出的一般是 0°、45°、90°、135° 四个方向的偏振强度图也有的输出 0°、60°、120° 三方向格式。计算 DoP 和 AOP 之前先把四张图转成 Stokes 参数这是最通用的一步import numpy as np from PIL import Image def load_polarized_frames(paths): 读取四方向偏振强度图返回 float32 数组 frames [] for p in paths: img np.array(Image.open(p)).astype(np.float32) frames.append(img) return frames def stokes_from_frames(frames): 输入 [I0, I45, I90, I135]输出 S0, S1, S2 I0, I45, I90, I135 frames S0 (I0 I45 I90 I135) / 2.0 S1 I0 - I90 S2 I45 - I135 return S0, S1, S2 def dop_aop_from_stokes(S0, S1, S2): 计算偏振度 DoP 与偏振角 AOP eps 1e-6 DoP np.sqrt(S1**2 S2**2) / (S0 eps) AOP 0.5 * np.arctan2(S2, S1) return DoP, AOP这段代码的思路是先把四方向强度图转成 S0、S1、S2 三个 Stokes 分量再从这个中间表示计算 DoP 和 AOP。S0 近似等于总光强S1 体现 0° 与 90° 方向上的偏振差异S2 体现 45° 与 135° 方向上的差异。arctan2返回的是弧度制范围在[-π, π]之间后面做彩色显示或者特征提取时要注意这个范围不能直接当作 0 到 255 的灰度图处理。2.3 几个直接影响后续融合质量的参数计算偏振参数时几个小参数会直接影响融合结果。第一是eps暗区里 S0 接近 0DoP 会爆炸式放大噪声我在实践中一般把eps设成整张图 S0 最大值的 1e-3 量级而不是固定数值。第二是输入图像位深很多偏振相机输出的是 12bit 或 16bit 的 tiff直接读进来当 8bit 用会丢失暗部细节。第三是坏像素校正偏振相机的超像素结构导致相邻像素对应不同偏振方向S1、S2 对错位极其敏感如果原始图像是 Bayer 格式必须先做插值解马赛克再计算 Stokes。提示偏振度和偏振角的计算顺序应当先配准后计算。如果四张不同方向的图之间存在几个像素的偏移先配准再算Stokes否则融合阶段会出现重影后面排查起来非常耗时。3. 把 zip 里的偏振融合跑起来目录结构、环境与核心流程3.1 解包之后先读这两个文件下载到的qzw3.zip解压后不复杂核心就两部分一组偏振原图和一组可运行的融合脚本。原图通常以frame_0.tif、frame_45.tif、frame_90.tif、frame_135.tif这样的命名规则组织对应四个偏振方向的采集结果。脚本部分会看到两个关键文件stokes.py负责从四方向强度图计算偏振参数fusion.py负责将强度图、DoP图、AOP图合成最终结果。动手改代码之前先把这两个文件单独摘出来跑一遍确认当前环境能正常执行再决定往哪个方向调。3.2 环境依赖与常见版本组合这个包的核心依赖集中在 NumPy、OpenCV 和 PyTorch 三个库上。NumPy 负责数组运算OpenCV 负责图像读取、金字塔和下采样PyTorch 只在融合部分用到了它的张量操作做权重计算不涉及模型训练。我本地环境是 Python 3.9 配 NumPy 1.24、OpenCV 4.8、PyTorch 2.1直接能跑通。pip install numpy opencv-python-headless torch --index-url https://download.pytorch.org/whl/cpuopencv-python-headless比完整版更干净适合服务器上跑不需要 GUI 支持。CPU 版 PyTorch 在纯推理场景下完全够用不必非要装 CUDA 版省两百多 MB 空间。3.3 从四方向强度图到融合输入跑通整个流程之前先把原始数据转成融合阶段需要的三个输入。这个脚本做的事情是读取四张原始帧计算 Stokes 参数导出 S0 强度图、归一化后的 DoP 图、映射后的 AOP 图分别存成 npy 文件方便后续融合阶段快速加载。import numpy as np from PIL import Image from glob import glob def prepare_inputs(data_dir, output_dir): # 按角度排序读取四张原始帧 paths sorted(glob(f{data_dir}/frame_*.tif)) frames [np.array(Image.open(p)).astype(np.float32) for p in paths] # 计算 Stokes 参数 I0, I45, I90, I135 frames S0 (I0 I45 I90 I135) / 2.0 S1 I0 - I90 S2 I45 - I135 # 计算偏振度加 epsilon 防止暗区除零 eps 1e-3 * np.max(S0) DoP np.sqrt(S1**2 S2**2) / (S0 eps) DoP np.clip(DoP, 0, 1) # 归一化到 0-1 # 计算偏振角映射到 0-255 便于可视化 AOP 0.5 * np.arctan2(S2, S1) AOP_norm ((AOP np.pi/2) / np.pi * 255).astype(np.uint8) # 强度图做一次简单增强方便融合 S0_norm (S0 - S0.min()) / (S0.max() - S0.min()) np.save(f{output_dir}/intensity.npy, S0_norm) np.save(f{output_dir}/dop.npy, DoP) np.save(f{output_dir}/aop.npy, AOP_norm) print(fSaved: intensity, dop, aop to {output_dir}) prepare_inputs(./data/raw, ./data/processed)这里有几个参数值得细说。eps的基准是 S0 的最大值而不是固定常数因为不同场景下的光强差异可以相差两个数量级固定eps1e-6在暗光场景会失效。np.clip(DoP, 0, 1)是强制归一化物理上 DoP 不会超过 1但噪声和计算误差可能把它推到 1 以上融合前钳位能避免极端值拉爆权重。AOP 映射到[0, 255]只是为了可视化实际做特征提取时建议保留弧度值因为线性映射会改变角度之间的相对距离影响后续的梯度分析和边缘检测。4. 融合策略选型像素级、特征级与决策级的三条路4.1 像素级融合的权重公式偏振融合最常用的是像素级加权平均因为它简单、可解释性强、参数少。核心问题是权重怎么设。有人直接把三个输入求平均结果往往是灰蒙蒙一片因为 DoP 图的数值范围和强度图不匹配。正确的做法是先归一化再做融合权重按场景用途来定import numpy as np def fuse_by_weight(intensity, dop, aop, w_intensity0.5, w_dop0.3, w_aop0.2): 像素级加权融合输入均为 float32范围 0-1 # 三个输入先做各自的最大最小值归一化 def normalize(x): return (x - x.min()) / (x.max() - x.min() 1e-6) intensity_norm normalize(intensity) dop_norm normalize(dop) aop_norm normalize(aop.astype(np.float32)) # AOP 如果是 uint8 先转 float # 归一化权重总和为 1 w_total w_intensity w_dop w_aop fused (w_intensity * intensity_norm w_dop * dop_norm w_aop * aop_norm) / w_total return np.clip(fused, 0, 1)这套权重组合默认给强度图最大占比 0.5因为强度图在大多数场景里仍然是信息主体DoP 给 0.3 突出材质差异AOP 给 0.2 做几何轮廓补充。如果你做的是镜面检测场景建议把 DoP 权重提到 0.5因为高光区域的偏振度特征最显著。如果是远距离遥感场景AOP 因为噪声大可以降到 0.1 以下。像素级融合的短板在暗部。强度图经过归一化后暗部被压缩DoP 在暗部本来就不稳定两个弱信号叠加后容易出现色块断层。遇到这种情况我不会继续调权重而是切到金字塔融合用多尺度分解把高低频分开处理。4.2 金字塔融合把低频、边缘和噪声分开处理拉普拉斯金字塔融合是偏振融合里性价比最高的方案。思想很简单图像分解成不同尺度的带通分量高频部分保留边缘和细节低频部分保留整体亮度分布融合时对每一层单独做决策。这样做的好处是强度图的高频边缘和你想要保留的 DoP 纹理不会互相污染。import cv2 import numpy as np def laplacian_pyramid(img, levels4): 构建拉普拉斯金字塔返回从高层到底层的列表 gaussian img.astype(np.float32) pyramid [] for _ in range(levels): # 下采样缩小一半得到低频近似 down cv2.pyrDown(gaussian) # 上采样回原尺寸和原图做差得到带通分量 up cv2.pyrUp(down, dstsize(gaussian.shape[1], gaussian.shape[0])) layer gaussian - up pyramid.append(layer) gaussian down pyramid.append(gaussian) return pyramid def fuse_pyramids(pyr_a, pyr_b, levels4): 对每一层取局部能量更大的像素 fused [] for la, lb in zip(pyr_a, pyr_b): # 用拉普拉斯能量作为显著性度量 energy_a cv2.convertScaleAbs(la)**2 energy_b cv2.convertScaleAbs(lb)**2 # 对能量做一次高斯模糊避免逐像素震荡 energy_a cv2.GaussianBlur(energy_a, (5, 5), 0) energy_b cv2.GaussianBlur(energy_b, (5, 5), 0) mask energy_a energy_b fused.append(np.where(mask, la, lb)) return fused def reconstruct(fused_pyr): 从金字塔重建融合结果 img fused_pyr[-1] for layer in reversed(fused_pyr[:-1]): img cv2.pyrUp(img) img cv2.add(img, layer) return np.clip(img, 0, 1)levels这个参数直接影响融合行为。默认 4 层对应图像 16 倍尺度跨度从单个像素到整体亮度都覆盖到了。如果原图分辨率是 1024×10244 层金字塔里最底层只剩 64×64是纯低频的“底子”这时候高频层保存的是边缘、噪声和细纹理。GaussianBlur这一步容易被跳过但它在逐像素能量比较里起了稳定作用没有它的融合结果会在边缘处出现像素级闪烁视觉上像是细碎的白点。4.3 特征级、决策级什么时候才值得切换像素级融合适合两三年内的小型项目参数直观、坑少。特征级融合比如 PCA 或稀疏表示需要在融合前分别从三个输入里提取特征向量计算量上了一个台阶但带来一个实际好处可以显式滤掉 DoP 里的暗部噪声。决策级融合一般配合语义分割或显著性检测先判断每个区域“用哪张图最可信”再按区域合成。如果你的场景里目标区域有明确的物理解释比如水面、镜面、草地决策级融合效果最稳。如果分不清场景里是什么决策级反而会因为误分割引入硬边界。从这份资源的代码结构看它主要支持像素级和金字塔融合特征级部分只提供了一组 PCA 权重的算例没有现成的稀疏表示模块。我一般建议的顺序是先跑加权平均得到基线结果再升级到金字塔最后根据实际场景判断值不值得上特征级。5. 偏振融合避坑实录四个高发问题与排查顺序5.1 现象一融合结果整体发灰、对比度起不来现象三路输入各自都正常加权融合后的图像像蒙了一层雾直方图集中在中间 30% 区间。原因最常见的是没有做归一化或者归一化顺序错了。强度图是 12bit 原始数据直接转 float32DoP 图已经是 0-1AOP 图还是弧度值三个不同量纲的数组直接乘权相加数值大的通道把其他通道的信息压没了。第二个常见原因是权重分配不合理强度图权重过高导致 DoP 的局部对比度贡献被淹没。解决先统一量纲再做融合三个输入各自独立做 min-max 归一化做完之后再乘权重。这是最省事且有效的一步八成以上的发灰问题都能靠它解决。如果归一化之后还是灰把w_dop提高 0.1、w_intensity降 0.1 试试DoP 图的纹理对比度对整体清透感贡献最大。5.2 现象二融合结果边缘出现双影、轮廓不干净原因双影几乎都是偏振相机采集阶段出了问题。四方向偏振片采集不是同时完成的如果场景里目标在动或者手持相机有微抖动四个方向的强度图之间存在亚像素级别的错位。而 Stokes 参数对错位极敏感S1、S2 是从两张不同方向的图直接相减得到的一个像素的错位就会在边缘处产生正负交替的错误值DoP 计算出来边缘处会多出一条亮线。解决计算 Stokes 参数之前先做配准。常见的做法是对 I0、I45、I90、I135 四张图做相位相关配准得到亚像素平移量再对齐。配准精度要达到 0.1 像素量级才能满足偏振计算的要求普通特征点配准的精度不够。import cv2 import numpy as np def align_frames(frames, reference_idx0): 使用相位相关法将其他通道对齐到参考帧 ref frames[reference_idx] aligned [ref] for i, img in enumerate(frames): if i reference_idx: continue # 计算亚像素平移量 shift, response cv2.phaseCorrelate( np.float32(ref), np.float32(img) ) # 根据平移量作仿射平移 rows, cols img.shape M np.float32([[1, 0, shift[0]], [0, 1, shift[1]]]) aligned_img cv2.warpAffine( img, M, (cols, rows), flagscv2.INTER_LINEAR cv2.WARP_FILL_OUTLIERS ) aligned.append(aligned_img) return alignedcv2.phaseCorrelate返回的是亚像素级别的位移向量就算原图只偏了 0.3 个像素也能测出来。配准的参考帧我一般选 I0因为它是重建 S0 的主体。如果配准之后还是有重影排查顺序是先确认目标有没有运动再确认相机有没有机械振动实在不行把各方向曝光时间拉短从源头上减少运动模糊。5.3 现象三DoP 图大量细碎噪点融合后呈椒盐状原因eps参数设置不合理会直接毁掉 DoP 图。暗光环境下S0 在某些暗区可能只有几十的灰度值此时如果按 S0 最大值乘以固定系数去设eps暗区里的除零问题会导致 DoP 变成随机值。还有一种情况是输入图像有坏点和死像素相邻像素间的异常值在 S1、S2 减法运算中被放大。解决把eps设置成基于全局亮度的自适应值对 S0 极小的区域单独处理。我固定用eps np.max(S0) * 1e-3作为下限同时将 S0 小于eps*3的像素直接标记为无效区域融合时跳过这些像素。这几条如果都做了还是噪点多就要检查是不是原始数据本身有偏暗或者偏亮的坏像素特别是 0°、90° 方向通道如果坏点位置不一致S1 会整片异常。注意偏振相机的超像素结构意味着每个 2×2 像素块共享同一时刻的光信号坏点会带着同方向的一整组偏振通道一起坏掉。做完坏点校正之后再算Stokes否则融合结果再怎么调参数都救不回来。5.4 现象四AOP 图显示为乱七八糟的花纹不像平滑的连续图原因AOP 的范围是[-π/2, π/2]直接当作灰度图显示时角度跨越 -π/2 到 π/2 的那一行像素会从黑跳到白形成肉眼可见的“伪边缘”。这其实是显示问题不是数据错误。另外arctan2返回的角度在接近 ±π/2 的边界处会自然出现跳变如果场景里理想偏振面正好在这个位置就会看到一圈圈的条纹。解决不要把 AOP 当作灰度图显示改用颜色来表示方向。物理上 AOP 是周期量应该用色相环表示然后在融合时把 AOP 先转成单位向量分量cos 和 sin而不是直接用角度值参与融合。def aop_to_embedding(aop_rad): 把偏振角转成两个连续分量避免角度边界跳变 cos_2a np.cos(2 * aop_rad) # 二倍角因为偏振角周期性为 π sin_2a np.sin(2 * aop_rad) return cos_2a, sin_2a用二倍角的原因是偏振角是 π 周期量0 和 π 表示同一个偏振方向直接用原始角度会让 0 和 π 在数值上相差 3.14但物理上它们等价。转成 cos、sin 两个分量后这个周期性问题就消失了融合时用这两个分量代替原来的 AOP 灰阶图后续不管做特征提取还是聚类都不会遇到跳变伪影。5.5 一条排查顺序建议偏振融合出问题时的排查顺序很重要。先查数据对齐再查归一化最后查融合权重。这个顺序不是随便拍的数据对齐决定了上游物理量是否准确归一化决定了数值尺度是否合理权重决定了信息是否被有效利用。按这个顺序排查基本能在五步之内定位到问题不至于在一个无关的权重参数上调半天。6. 偏振融合效果验证三个指标与一个固定动作融合效果“看起来清楚了”是不够的主观视觉容易掩盖算法层面的问题。我用三个客观指标来验证每个指标都有明确的计算方法和物理含义。第一个是信息熵。融合结果的信息熵比任何单通道输入熵高说明融合确实带来了信息增量。计算公式是-Σ p(i) * log2(p(i))其中p(i)是像素灰度级的出现概率。理想状态下熵值高两个点左右就算合格。第二个是平均梯度也叫清晰度指标对融合结果做 Sobel 算子提取梯度幅值后取均值数值越高说明边缘越锐利纹理越清晰。第三个是边缘保持度 QAB/F在 0 到 1 之间衡量融合结果里保留了多少输入图的边缘信息。对于偏振融合场景QAB/F 在 0.6 以上就是不错的成绩。import numpy as np from scipy.stats import entropy import cv2 def evaluate_fusion(fused, source_channels): 验证融合质量的三个指标 # 指标 1: 信息熵 hist np.histogram(fused.ravel(), bins256, range(0, 1))[0] hist hist / hist.sum() entropy_val entropy(hist) # 指标 2: 平均梯度 gx cv2.Sobel(fused, cv2.CV_64F, 1, 0, ksize3) gy cv2.Sobel(fused, cv2.CV_64F, 0, 1, ksize3) grad_mag np.sqrt(gx**2 gy**2) avg_grad grad_mag.mean() # 指标 3: 与各输入的结构相似度均值 from skimage.metrics import structural_similarity as ssim ssim_scores [ ssim(fused, source, data_range1.0) for source in source_channels ] avg_ssim np.mean(ssim_scores) return { entropy: round(entropy_val, 4), avg_gradient: round(float(avg_grad), 4), avg_ssim: round(avg_ssim, 4) }这套验证脚本跑完如果熵值提升了但平均梯度下降了说明融合把细节磨平了需要增大高频层的能量选择权重。如果平均梯度高但熵值低说明结果布满颗粒噪声需要回到 DoP 计算的eps设置上找原因。三个指标之间互相校验基本能判断出问题出在融合策略还是上游数据。从那以后我每次做偏振融合都在固定动作里加一步先把四通道原图的配准误差数值和趋势打出来再跑融合。配准误差大于 0.5 像素的帧直接丢弃不参与计算与其让错位的数据带崩融合结果不如在源头就过滤掉。偏振信息说到底是个精密的物理量任何一个环节马虎后面都会用各种奇怪的方式还回来。这套流程现在成了我项目里的标配希望帮到你。本文还有配套的精品资源点击获取