简介：偏振图像分析工具包面向光学成像、遥感、生物医学成像等领域的开发者和研究人员，整合了去马赛克、Stokes向量计算与Mueller矩阵分析等核心功能，可帮助用户从拜耳格式原始数据出发完成偏振信息提取与可视化。包内共32个文件，以Python脚本为主（11个py），配合示例图片、图表和说明文档，覆盖从demosaicing预处理到stokes/mueller计算及结果绘图的完整流程，压缩包仅2.29MB，轻量易部署。已有411人学习使用。通过该工具，读者不仅能直接运行现成脚本处理自己的偏振图像，还可参考其代码结构理解偏振光学模型，并借助SVG示意图与示例输出快速验证算法效果，适合希望深入掌握偏振成像技术的入门及进阶者。
1. 偏振图像分析工具：把传感器棋盘格还原成物理量
偏振相机拍到的不是一张普通照片。它的每个像素前面都贴着一片不同角度的微偏振片，在一块 2×2 的棋盘格里，四个像素分别只让 0°、45°、90°、135° 方向的光通过。于是 raw 图天生就是欠采样的：某个像素的强度只代表一个偏振方向，另外三个方向必须靠邻域插值补出来。这一步做不好，后面算出来的 Stokes 向量和 Mueller 矩阵全都会偏。
偏振图像分析工具要做的，就是把这种四分之一分辨率的强度图还原成像素级对齐的物理量。先做去马赛克，把 2×2 超级像素展开成四个全分辨率通道；再按 Stokes 定义合成 S0、S1、S2；如果还要测双折射、退偏或散射特性，就用多组偏振态主动照明，反演 Mueller 矩阵。三步各有误差源，上一步的误差还会继续向后传播。这篇文章按这套流水线讲清楚每级做什么、参数怎么选、哪些坑是偏振场景特有的。
适合读这篇文章的，是做工业缺陷检测、生物组织成像、遥感目标识别或者光学材料评测的工程师。如果你手里已经有一台分焦平面偏振相机，或者刚拿到一个类似“去马赛克 + Stokes + Mueller”的处理模块，这里讲的就是你实际会重建的那条链路。
2. 偏振去马赛克如何处理 2×2 超级像素：三种插值法与串扰抑制
2.1 分焦平面传感器的 2×2 偏振方向分布与像素命名
分焦平面（Division of Focal Plane，DoFP）是目前工业偏振相机的主流结构。常见排布是：偶数行偶数列放 0°，偶数行奇数列放 45°，奇数行偶数列放 135°，奇数行奇数列放 90°。四个方向各占原始分辨率的四分之一，所以 raw 图本质上是四张半分辨率图像的交错拼接。
不同厂商的排布顺序可能不同，有些传感器把 135° 放在偶数行奇数列，处理前必须先查数据手册。我一般会在读取 raw 后先做一个通道自检：把镜头对准均匀照明的散射面，分别统计四个方向通道的均值。若某个通道均值明显偏离其他三个，先怀疑排布定义错了，而不是传感器坏了。
排布定义直接影响插值公式。在写代码之前，把四个通道的物理方向固定成一个元组polar_order=(0, 45, 135, 90)，后续所有索引都从这里推导，能避免一半的坐标错误。
2.2 双线性、双三次和边缘导向插值：取舍表与适用场景
偏振去马赛克和彩色拜耳去马赛克有一个本质差异：彩色通道之间是光谱冗余，而偏振通道之间是方向冗余，目标物体上的某个亮点在 0° 通道可能很亮、在 90° 通道却很暗。因此，插值算法的目标不是“看起来平滑”，而是尽可能不引入虚假的偏振差分信号。
| 算法 | 计算成本 | 边缘表现 | 对偏振参数的影响 | 适用场景 |
|---|---|---|---|---|
| 双线性 | 低 | 边缘模糊 | S1、S2 被低估，DOLP 偏小 | 快速预览、作为基线 |
| 双三次 | 中 | 略好 | 改善有限，仍会平滑跳变 | 默认选择，多数情况可用 |
| 边缘导向 | 高 | 保留强度边缘 | 在偏振差分图上容易过冲 | 纹理清晰、边缘锐利的工业件 |
| 频域滤波 | 高 | 依赖载波频率 | 对图案频率敏感，易产生振铃 | 特定光学图案，不推荐通用使用 |
双线性的低估可以从频域理解：0° 通道每两个像素才有一个采样点，插值相当于低通滤波，把高频的偏振差分信号平滑掉了。对一个 DOLP=0.8 的均匀偏振区域，双线性插值后 DOLP 可能降到 0.6 左右。这不是随机噪声，是系统偏差。
我一般把双三次作为默认，边缘导向留给特定样品。如果你只需要测 Mueller 矩阵，插值引起的高频损失可以用后续的多帧平均补偿，且双线性的线性特性让误差更容易建模。
2.3 用 Python 实现可复现的偏振去马赛克最小脚本
下面这段代码处理单张 uint16 raw，输出四个与输入同尺寸的通道。它用 2×2 复制加均值核实现双线性插值，边缘用最近邻延拓，避免 scipy.ndimage.zoom 带来的半像素偏移。
import numpy as np from scipy.ndimage import convolve def demosaic_linear(raw, polar_order=(0, 45, 135, 90)): """DoFP 偏振去马赛克，双线性实现。 raw: 偶数行偶数列为 polar_order[0]，以此类推。 返回四个与 raw 同尺寸的浮点强度图。 """ h, w = raw.shape # 按物理方向索引拆出半分辨率通道 I0 = raw[0::2, 0::2].astype(np.float32) # 0° I45 = raw[0::2, 1::2].astype(np.float32) # 45° I135 = raw[1::2, 0::2].astype(np.float32) # 135° I90 = raw[1::2, 1::2].astype(np.float32) # 90° def upsample(img): """双线性上采样：复制像素 + 2x2 均值滤波。""" H, W = img.shape # 右、下各补一行，使复制后尺寸正好翻倍 img = np.pad(img, ((0, 1), (0, 1)), mode='edge') bigger = np.repeat(np.repeat(img, 2, axis=0), 2, axis=1) kernel = np.ones((2, 2), dtype=np.float32) / 4.0 return convolve(bigger, kernel, mode='nearest') return upsample(I0), upsample(I45), upsample(I135), upsample(I90)逻辑说明：先按 2×2 步长切片，把交错排布的四通道拆开；随后每个半分辨率图先补边再像素复制，复制结果经过 2×2 均值核，等效于一次双线性插值。mean 核的权重各为 0.25，没有引入方向偏置。
参数说明：raw 建议先转 float32，uint16 的除法和卷积会丢精度；polar_order 只是可读性约束，真正决定通道位置的是切片坐标。若你的传感器排布不同，改切片下标即可。
边缘处理很关键。mode='edge' 会让图像最外两行参与插值时不产生颜色失真，而默认的 zero padding 会在暗场边缘制造偏振条纹，肉眼可能看不到，但 Stokes 图上会多出一圈假偏振。
2.4 从 Stokes 角度检查插值伪影：哪些坑值得注意
插值结果好不好，不看单通道平滑度，要看合成后的 S0 是否平坦。对一个均匀照明的偏振区域，S0 应该只反映光强，与偏振方向无关。插值有问题时，S0 上会出现 2×2 棋盘纹理，也就是原图每两个像素亮暗交替。
一个简单的定量指标是棋盘差：把去马赛克后的 S0 按 2×2 窗口计算方差，再与全局方差比较。如果窗口内方差显著高于全局，说明插值引入了通道间增益不一致，需要检查平场校正或排布。另一个常见坑是插值后出现负强度。负值会让 S1 的绝对值大于 S0，最终 DOLP 大于 1。发现负值时不要直接截断成 0，先检查 raw 是否做过暗帧扣除、插值核是否越界。
去马赛克完成后，四个通道的像素网格要严格一一对应。如果四个通道之间有亚像素偏移，会在 Stokes 差分时产生边缘重影。后续先用 S0 或高对比靶标检查对齐，再做合成。
3. Stokes 向量计算：从四强度图到 DOLP 与 AOLP
3.1 Stokes 参数的定义：为什么 S0、S1、S2 是线偏振相机能拿到的全部
Stokes 向量是描述光偏振态的四个实数，通常写成 [S0, S1, S2, S3]^T。S0 是总光强；S1 表示 0° 与 90° 方向上的强度差，反映水平/垂直偏振；S2 表示 45° 与 135° 方向上的强度差，反映对角偏振；S3 表示左旋与右旋圆偏振的差。
分焦平面线偏振相机的每个像素只测了三个方向上的四档强度，且所有微偏振片都是线偏振片，无法直接区分圆偏振分量。所以对 DoFP 相机，S3 信息缺失，通常当作 0 处理。这不是理论近似，而是传感器物理结构决定的。
从四方向强度到 Stokes 的换算有两种常见约定。严格定义是 S0 = I0 + I90 或 S0 = I45 + I135，但在真实传感器上四个通道的量子效率和透过率总有差异，取四通道平均更稳：S0 = (I0 + I45 + I90 + I135) / 2。S1 = I0 - I90，S2 = I45 - I135。这个约定在多数软件库里都能对上。
3.2 逐像素计算 S0、S1、S2 的 NumPy 实现
去马赛克输出四个同尺寸通道后，Stokes 合成只是几行 NumPy 运算。代码里加了三个小处理：S0 下限保护、AOLP 范围归一化、以及低于光强阈值的像素掩膜。
import numpy as np def stokes_from_demosaic(I0, I45, I90, I135, eps=1e-8): """从四个偏振方向强度计算线偏振部分 Stokes 参数。 返回 S0, S1, S2, DOLP, AOLP，AOLP 单位弧度，范围 [-pi/2, pi/2]。 """ S0 = 0.5 * (I0 + I45 + I90 + I135) S1 = I0 - I90 S2 = I45 - I135 # DOLP 需要防止零光强处除零 s0_safe = np.maximum(S0, eps) dolp = np.sqrt(S1 * S1 + S2 * S2) / s0_safe aolp = 0.5 * np.arctan2(S2, S1) # 低强度像素的偏振方向没有意义，直接置为无效 low_light_mask = S0 < (0.02 * np.nanmax(S0)) aolp[low_light_mask] = np.nan return S0, S1, S2, dolp, aolp逻辑说明：S0 取四个通道平均，能抵消一部分通道间增益不匹配；S1、S2 是差分，DoFP 传感器上相邻像素的暗电流如果不一致，会直接污染这两个分量。arctan2 自动处理 S1 为 0 时的象限问题，0.5 系数让 AOLP 角度的物理范围落在 ±90° 而不是 ±180°。
参数说明：eps 只用来避免除零，不要设置太大，否则低光区域 DOLP 会被人为压到 0。low_light_mask 用全局最大值的 2% 作为阈值，这个值在暗室和室外高动态场景需要单独调整。把无效 AOLP 设成 NaN 比设成 0 更安全，后续可视化时 NaN 像素可以直接挖掉，而 0 会变成一片假偏振方向。
3.3 DOLP 超过 1、负强度出现：误差的处置原则与噪声传播
理想情况下 DOLP 永远在 [0, 1] 区间，但实测数据经常冒出大于 1 的像素。原因主要有三个：暗帧扣除不彻底导致差分通道有偏置；插值产生负强度使 S1、S2 被放大；以及随机噪声在低光区域主导了差分结果。
对 DOLP 超过 1 的像素，常见的做法是 clip 到 1。我一般不建议直接 clip，因为这会掩盖上游的标定问题。正确做法是回到插值前检查 raw 的暗电平，再检查四个通道的直方图。若只有零星像素超 1，多为噪声，clip 无害；若成片分布且形状规则，多半是偏振片方向标定错了，clip 只会把系统误差藏起来。
Stokes 参数的噪声传播还有一个容易被忽略的点：DOLP 是 S1、S2 绝对值的非线性函数，在低信噪比区域，噪声会被绝对值运算整流成正值，产生“即便是完全非偏振光也能测出虚假 DOLP”的现象。要验证这一点，把镜头盖上拍一张暗帧，计算暗帧的 DOLP。如果暗帧 DOLP 不是接近 0 而是有一大片高值，说明你的 DOLP 计算没有做光强权重处理，低光区域结果不可信。
4. Mueller 矩阵的标定与最小二乘反演：从实验布局到矩阵约束
4.1 偏振态生成器和分析器的基本布局：只靠线偏振相机能恢复多少
Mueller 矩阵是一个 4×4 实数矩阵，把入射光的 Stokes 向量变换为出射光的 Stokes 向量。要测它，必须在样品前放一个偏振态生成器（PSG），在样品后放一个偏振态分析器（PSA）。PSG 负责产生多组已知的入射 Stokes 态，PSA 负责测量出射的 Stokes 态。
实际布局里，PSG 通常用“线偏振器 + 可旋转四分之一波片”或液晶可变延迟器实现，能产生包含圆偏振分量的入射态。PSA 这一侧如果只是把 DoFP 线偏振相机直接放在样品后面，那只能测出 Stokes 向量的前三个分量，S3 完全没有。这种情况下，Mueller 矩阵中与 S3 相关的第四行、第四列无法反演，你能得到的只是一个“降维”的线偏振 Mueller 子矩阵。
要恢复完整的 4×4 Mueller 矩阵，PSA 一侧必须有能力测量 S3，常见方案是在相机前再加一旋转波片，把圆偏振信息转化为线偏振信息。这意味着要多次采集，单帧 DoFP 相机做不到。如果你的应用只需要二向色性、线性双折射这类线偏振相关参数，降维矩阵通常够用；需要测退偏或圆二色性时，必须上完整 PSA。
4.2 输入态组合、最小二乘反演与条件数控制
反演 Mueller 矩阵的最小二乘形式很简洁。把 N 组入射 Stokes 并成 4×N 矩阵 S_in，对应出射 Stokes 并成 4×N 矩阵 S_out，关系是 S_out = M @ S_in。M 的最小二乘解是 M = S_out @ S_in^T @ (S_in @ S_in^T)^(-1)，在 Python 里用 np.linalg.pinv 更稳。
import numpy as np def invert_mueller(stokes_in, stokes_out): """从多组入射/出射 Stokes 对反演 4x4 Mueller 矩阵。 stokes_in: 4xN，每一列是一组入射 Stokes。 stokes_out: 4xN，每一列是对应出射 Stokes。 """ # 使用伪逆以避免 S_in 行相关时直接求逆失败 M_est = stokes_out @ np.linalg.pinv(stokes_in) return M_est # 四组常见入射态：0° 线偏、90° 线偏、45° 线偏、右旋圆偏 stokes_in = np.array([ [1, 1, 0, 0], [1, -1, 0, 0], [1, 0, 1, 0], [1, 0, 0, 1], ], dtype=np.float64).T # 假设一个真实 Mueller 矩阵（这里虚构，仅用于演示反演流程） M_true = np.array([ [1.0, 0.5, 0.2, 0.0], [0.4, -0.3, 0.1, 0.0], [0.1, 0.0, 0.6, 0.0], [0.0, 0.0, 0.0, 0.8], ]) stokes_out = M_true @ stokes_in M_est = invert_mueller(stokes_in, stokes_out) print(np.round(M_est, 6))逻辑说明：stokes_in 的每一列是一个已知的入射偏振态，四组态在 4×4 空间里构成一组基。stokes_out 是对应的出射测量值。伪逆在 N 大于 4 时自动做最小二乘平均，能抑制一部分随机测量噪声。
参数说明：入射态组合选择直接影响反演精度。对 stokes_in 做奇异值分解，看奇异值是否均衡；条件数是最大奇异值与最小奇异值之比，越小越好，一般要求低于 5。上面这组态的条件数约为 2.0，是很稳的配置。如果入射态里只有线偏振没有圆偏振，第四行接近全 0，条件数会急剧升高，这时 M 的第四列完全不可辨识，反演结果会出现巨大的伪值。
实测时不要只测 4 组态，建议测 8 到 12 组并按伪逆计算。多出来的约束会让单次测量误差被平均掉，代价只是采集时间增加。PSG 旋转角度定位误差是最大的系统误差来源，旋转台每步误差 0.5° 会让 S1、S2 出现百分比级别的偏差。
4.3 Mueller 矩阵的物理可允许性检查：退偏度、二向色性与相位延迟
反演出来的矩阵不一定物理可实现。Mueller 矩阵必须满足“任意物理上可实现的入射 Stokes 向量，映射后仍是物理可实现的 Stokes 向量”这一约束。最简单的检查方法是随机生成几千组实测范围内的入射态，计算输出 DOLP，若出现显著大于 1 的输出，说明反演矩阵有问题。
若想进一步分析样品特性，可以按 Lu-Chipman 分解把 Mueller 矩阵拆成三个级联部分：退偏器、二向色衰减器、相位延迟器。退偏器对应矩阵的“平均化”程度；二向色性表示样品对不同偏振方向的吸收差异；相位延迟对应双折射。对一块透明玻璃，延迟和退偏应该接近 0；对生物组织切片，退偏项通常明显非零。
检查时不要只看 M00 是否接近整体透过率。M00 只代表强度响应，而 S1、S2 所在的行列才是偏振敏感信息。把 M_est 打印出来后，先看右下 3×3 块的对称性。很多实际样品在无偏振放大时应该近似各向同性，若右下 3×3 出现明显的非对称结构，优先怀疑入射态标定错误，而不是样品本身有手性。
5. 工程化验证：偏振通道对齐、平场修正与已知样品复盘
5.1 通道配对、坏点掩膜和均匀光照标定
在跑任何 Stokes 或 Mueller 处理前，先把四通道平场做好。把相机对准积分球或干净白纸，调整光源到约 80% 满阱，采集多帧平均，得到每个通道的响应斜率。平场修正不是简单的四通道统一乘系数，因为每个 2×2 单元里的四个偏振片透过率不同，需要逐像素求增益图。
通道配对检查我常用一张高对比十字靶标：用普通白光照明，分别提取四通道图像，用相位相关计算相互偏移。若任意两通道之间有超过 0.3 像素的偏移，就要考虑做亚像素配准。配准应在去马赛克前对半分辨率通道做，不能在插值后再做，否则会引入重采样误差。
坏点处理也要放在去马赛克之前。DoFP 相机的坏点通常是单个偏振像素失效，若不处理，插值会把坏点扩散到四个通道，在 Stokes 图上形成一个 2×2 的假偏振斑。我一般用邻域同偏振方向的像素中值替换，而不是用跨方向像素插值，后者会直接污染 S1、S2。
5.2 用一片线性偏振片和一叠玻璃验证 Stokes 与 Mueller 结果的符号和幅度
验证计算链路上的符号错误，最经济的样品就是一片线性偏振片。把偏振片放在光源和相机之间，旋转到某角度 θ，理想情况下 AOLP 应该等于 θ 加上一个固定偏置。如果 AOLP 随 θ 线性变化但方向反号，说明 S1 或 S2 的定义与传感器排布对不上；如果 AOLP 不随 θ 变化，多半是四通道索引拆错了。
对 Mueller 测量，可以用一叠普通载玻片做粗略参考：几片玻璃叠在一起，相当于一个具有明显相位延迟的样品。测得的相位延迟项应随玻璃片数近似线性增加，并且旋转样品时延迟方向跟着旋转。这个测试能暴露 PSA/PSG 角度回零误差和坐标系不一致问题。
最后一次验证建议写成自动化脚本：加载平场修正后的 raw，依次跑去马赛克、Stokes 合成、Mueller 反演，然后断言 DOLP 在 [0, 1.05] 区间内、AOLP 与偏振片机械角度差小于 2°。断言失败时打印四个通道的均值，快速判断是暗帧问题还是通道定义问题。
本文还有配套的精品资源，点击获取
