简介偏振图像分析工具用Python实现了完整的偏振图像处理流程面向图像处理、光学检测、遥感、生物医学成像和材料科学等方向的开发者涵盖去马赛克、斯托克斯向量计算与穆勒矩阵分析三大任务便于理解偏振光变化规律并完成图像复原、偏振参数提取和目标材质差异分析。压缩包共25个文件大小为1.49MB其中包含11个源代码文件、9张示例照片、3张示意图、1份说明文档和许可证文件。源代码包括偏振图像去马赛克、斯托克斯参数求解、穆勒矩阵计算等功能模块示例图片展示了不同偏振角度下拍摄的强度图以及经偏振角图、偏振度图处理后的效果图说明文档则对代码结构、依赖环境和使用方法进行了梳理。目前已有1516人学习/下载。借助这些代码读者可以快速搭建偏振图像分析实验环境参考示例数据理解偏振参数的计算过程并在此基础上修改或扩展算法适合希望在偏振成像方向快速上手的科研人员和工程师。1. 偏振图像分析工具的第一步从马赛克原图到物理量打开dragon_IMX250MZR_intensity.jpg这类偏振相机原始图时第一眼往往不是彩色照片而是一张灰扑扑的马赛克。普通相机拜耳阵列分的是 R/G/B 颜色偏振相机比如 Sony IMX250MZR则把 0°、45°、90°、135° 四个偏振方向的微偏振片铺在相邻 2×2 像素上原始图需要先拆成四个方向的强度通道才能继续算斯托克斯向量和穆勒矩阵。Polanalyser 就是这套流程的 Python 实现demosaicing.py负责去马赛克stokes.py计算 S0/S1/S2 与 DoLP/AoLPmueller.py处理多偏振态测量下的 16 元素矩阵反解。它适合用 Python 做偏振成像算法验证、光学检测标定或材料偏振响应研究的人不需要自己从头造轮子。2. 去马赛克把偏振拜耳阵列拆成四个独立强度2.1 偏振马赛克与普通拜耳马赛克的本质区别普通相机的拜耳阵列在 2×2 单元里放 R、G、G、B 三种滤色片去马赛克的本质是插值用周围像素估计每个位置的完整 RGB。偏振相机的 2×2 单元则放 0°、45°、90°、135° 四种微偏振片每个像素只能测到某一个方向的线偏振强度。它没有“颜色”需要恢复要的是把同一场景在四个方向的强度分别取出来。 如果你把原图直接当灰度图用2×2 内的偏振调制会一直混在像素里后续算 S1 I0 − I90 时每个像素拿到的都可能是相邻方向的串扰结果完全失真。这里还有一个很多人没注意的代价空间分辨率减半。四个方向的像素共享同一块靶面拆开后每个方向的有效图片是原图宽度和高度各一半。也就是说偏振相机的“原生分辨率”和普通相机同像素时并不等价。Polanalyser 的demosaicing.py不打算做超分辨它的职责就是把像素按微偏振片位置挑选出来排成四通道强度图。理解了这点下面读代码就不会被形状变化吓到。2.2 用 demosaicing 的核心逻辑拆分强度图下载的polanalyser-master里demosaicing.py对外暴露的函数在源码里通常就叫demosaicing。它的核心计算并不复杂如果你只想快速跑通下面这个切片版本足够import numpy as np import cv2 def demosaic_polar(raw, patternimx250mzr): 拆分偏振马赛克图返回 (H/2, W/2, 4) 的 float32 强度图。 通道顺序0度、45度、90度、135度。 if pattern imx250mzr: I0 raw[0::2, 0::2].astype(np.float32) I45 raw[0::2, 1::2].astype(np.float32) I90 raw[1::2, 0::2].astype(np.float32) I135 raw[1::2, 1::2].astype(np.float32) else: raise ValueError(unknown polarization mosaic pattern) return np.stack([I0, I45, I90, I135], axis-1) raw cv2.imread(dragon_IMX250MZR_intensity.jpg, cv2.IMREAD_UNCHANGED) if raw.ndim 3: raw raw[..., 0] # 某些原图会以三通道形式保存但只用了同一通道 img demosaic_polar(raw) print(img.shape, img.dtype) # (H/2, W/2, 4) float32逻辑说明raw[0::2, 0::2]取的是偶数行偶数列对应 0° 微偏振片的位置raw[0::2, 1::2]对应 45°raw[1::2, 0::2]对应 90°raw[1::2, 1::2]对应 135°。切片本身是视图astype(np.float32)会复制数据这样后续做减法时不会因为 uint16 溢出。pattern参数留给不同相机的微偏振片排布差异IMX250MZR 是常见的 0/45/90/135 按行排列。需要强调读取方式用cv2.imread(path, cv2.IMREAD_COLOR)会把单通道马赛克复制成三通道等值图再去切片就会让四个方向完全相同。所以IMREAD_UNCHANGED不是可选项是必须项。工具包里的examples脚本有的还会用tifffile或imageio读 12bit 原图效果类似。2.3 参数、dtype 与常见坑项目推荐做法原因读取原图cv2.imread(path, cv2.IMREAD_UNCHANGED)默认 flag 会把 12bit 压成 8bit破坏偏振强度线性度数据 dtypeuint16 读入分割后立即转 float32S1/S2 是差值会出现负数uint16 无法表示输入形状必须是 (H, W) 或 (H, W, 1)三通道等值马赛克会让四个方向完全相同通道顺序先统一成 0/45/90/135stokes.py的 I0−I90 和 I45−I135 依赖这个顺序输出形状(H/2, W/2, 4)每个方向共用一帧空间分辨率减半一个常见坑是相机厂商的 SDK 可能已经把数据排成了四通道或把 2×2 块改成“交错行”。这种情况下继续用上面的切片会得到棋盘状伪影。我一般会先对均匀光照下的平面拍一张图拆完后看四个通道均值是否接近如果某个通道明显偏亮或出现周期性条纹通常是行列偏移了一位把1::2和0::2对调即可。2.4 验证拆分是否正确means img.reshape(-1, 4).mean(axis0) print(means) # 均匀散射区域四个通道均值应接近逻辑说明对一片没有明显偏振特性的区域DoLP 接近 0四个方向的强度应该近似相等。任何一个通道显著偏离说明马赛克相位没对齐。这个检查只需要几行却是后面所有计算正确性的地基。如果直接跳到斯托克斯计算S1/S2 的错误会被可视化图掩盖成“颜色有点怪”很难定位是算法问题还是像素错位。3. 斯托克斯向量计算从四方向强度到 S0、S1、S23.1 斯托克斯参数的线性定义斯托克斯向量把光的偏振态写成四个实数对线偏振成像场景通常只关心前三个S0 是总光强S1 表示水平与垂直线偏振的差S2 表示 45° 与 135° 线偏振的差。只要拿到了去马赛克后的 I0、I45、I90、I135S0/S1/S2 就是简单的加减法def calc_stokes(img): 输入 (H/2, W/2, 4)返回 (H/2, W/2, 3)。 I0, I45, I90, I135 img[..., 0], img[..., 1], img[..., 2], img[..., 3] S0 I0 I90 S1 I0 - I90 S2 I45 - I135 return np.stack([S0, S1, S2], axis-1)逻辑说明S0 用一对正交方向相加是因为 0° 和 90° 强度之和等于总强度S1、S2 则是两组正交方向的差分。源码包里的stokes.py基本就是这几行只是额外处理了输入通道顺序和 NaN。要注意坐标系定义有些教材把 S1 写成 I90 − I0这样可视化时 S1 的符号会左右翻转AoLP 也会偏移 90°。同一个包里demosaicing.py和stokes.py必须保持同一套约定我一般会在一开始写注释固定通道顺序。3.2 DoLP 与 AoLP归一化和弧度处理只拿到 S0/S1/S2 还不够实际分析经常要偏振度 DoLP 和偏振角 AoLP。DoLP 表示偏振光占总强度的比例范围 0 到 1AoLP 表示偏振主轴方向物理上只有 0 到 π 的周期。实现如下def calc_DoLP(stokes, eps1e-6): S0, S1, S2 stokes[..., 0], stokes[..., 1], stokes[..., 2] return np.sqrt(S1 ** 2 S2 ** 2) / np.maximum(S0, eps) def calc_AoLP(stokes): S1, S2 stokes[..., 1], stokes[..., 2] aolp 0.5 * np.arctan2(S2, S1) return np.mod(aolp, np.pi) # 主轴在 [0, π) 内参数说明eps1e-6是为了避免 S0 接近 0 时除零np.mod(aolp, np.pi)必须加因为0.5 * arctan2返回范围是 −π/2 到 π/2而偏振角是 90° 周期直接转成 0 到 π 的循环量才方便做伪彩色和统计。如果想要角度显示乘180 / np.pi注意是 0 到 180 度不是 0 到 360 度。3.3 视觉化examples/visualize_AoLP.py 的思路示例目录里的visualize_AoLP.py不是简单地显示灰度图而是把 AoLP 映射到色环颜色代表偏振方向亮度和饱和度分别表示光强和偏振度。这样一张图能同时看出三个维度。核心代码如下def visualize_stokes(stokes): S0 stokes[..., 0] DoLP calc_DoLP(stokes) AoLP calc_AoLP(stokes) hsv np.zeros((S0.shape[0], S0.shape[1], 3), dtypenp.uint8) hsv[..., 0] (AoLP / np.pi * 179).astype(np.uint8) hsv[..., 1] (np.clip(DoLP, 0, 1) * 255).astype(np.uint8) hsv[..., 2] (S0 / np.max(S0) * 255).astype(np.uint8) return cv2.cvtColor(hsv, cv2.COLOR_HSV2BGR)逻辑说明OpenCV 的 8bit HSV 模式中 Hue 范围是 0 到 179所以把 AoLP 的 0 到 π 映射到 0 到 179。Saturation 放 DoLPValue 放归一化光强 S0。dragon_IMX250MZR_AoLP_color.jpg这类示例图就是这么生成的。jpge保存时质量选高一些否则色环边缘会出现条带。3.4 用 S0 做自检斯托克斯计算完成后我会先做一个最小验证把 S0 单独显示出来它应当和原图经过 2×2 池化后的结果很像只是尺寸减半。如果 S0 出现棋盘格说明去马赛克的行列错位如果 S0 基本正确但 S1/S2 一片杂乱则要检查暗场和坏点。下面的代码可以替代人眼判断S0 stokes[..., 0] S0_expected img[..., 0] img[..., 2] # I0 I90 print(np.max(np.abs(S0 - S0_expected)))逻辑说明S0_expected直接从去马赛克结果里相加理论上必须与calc_stokes输出的 S0 完全一致。这个断言能在早期发现通道顺序写反、数组被原地修改等问题。注意这里不是和img.mean(axis-1)对比因为mean除以了 4。4. 穆勒矩阵计算从一组偏振态测量反解物体的偏振响应4.1 为什么一张图算不出穆勒矩阵斯托克斯向量描述光的偏振状态穆勒矩阵描述物体把入射偏振态变成出射偏振态的 4×4 线性变换。每个像素在不同入射偏振态下会给出不同出射强度单张偏振图只对应一次入射最多得到一个方程而一个未知物体有 16 个矩阵元素要解。因此至少需要四个线性无关的入射 Stokes 向量配合四个线性无关的分析态测量。Polanalyser 的mueller.py做的事情就是把一组测量强度图按最小二乘反解成每个像素的 4×4 矩阵。这里要区分两个概念去马赛克和斯托克斯计算都是“单帧处理”穆勒矩阵需要多帧。最简实验装置是入射侧放一个可旋转偏振片产生偏振态 PSG出射侧放另一个偏振片分析 PSA中间放被测样品。你下载的源码包里examples/mueller.py就是这个流程的演示。4.2 测量方程与线性反解理想线偏振片在角度 θ 处的 Stokes 向量是(1, cos2θ, sin2θ, 0)。如果入射 Stokes 是s_in样品矩阵是 M出射经过分析器后的强度满足I a_out^T · M · s_in把 M 按行优先展平成 16 维向量m则每个测量就是一个行向量kron(a_out, s_in)与m的内积。16 组不同偏振态组成 16×16 矩阵 A强度图展平成(16, N)用最小二乘一次解出所有像素的矩阵元素。以下是通用的方程构建函数import numpy as np def build_mueller_equations(S_in_list, A_out_list): 根据入射/分析斯托克斯向量构造测量矩阵 A。 A [] for s_in in S_in_list: for a_out in A_out_list: A.append(np.kron(a_out, s_in)) return np.array(A) states [ [1, 1, 0, 0], # 0° 线偏振 [1, 0, 1, 0], # 45° 线偏振 [1, 0, 0, 1], # 右旋圆偏振 [1, -1, 0, 0], # 90° 线偏振 ] A build_mueller_equations(states, states) print(condition number:, np.linalg.cond(A))逻辑说明states同时作为入射和分析器状态覆盖了线偏振和圆偏振保证 A 可逆。np.kron(a_out, s_in)的排列顺序与M.flatten()的展开方式一致kron中参数顺序不能反。如果只用 0/45/90/135 线偏振A 会奇异因为所有输入输出状态都落在 S30 的子空间穆勒矩阵的第三行和第三列永远解不出来。真实测量至少需要一组四分之一波片产生圆偏振态。4.3 像素级反解与内存控制拿到 A 矩阵后按测量顺序把强度图堆成(16, H, W)的数组再展平求解# I_meas 形状 (16, H, W)顺序必须和 build_mueller_equations 的循环顺序一致 I_flat I_meas.reshape(16, -1) M_flat np.linalg.lstsq(A, I_flat, rcondNone)[0] M M_flat.reshape(4, 4, H, W) # 查看某个区域的平均穆勒矩阵 print(M[..., :3, :3].mean(axis(2, 3)))参数说明I_flat是 16 行、N 列的矩阵每列对应一个像素的 16 次测量。lstsq在这种情况下比inv(A) I_flat更稳因为偏振片不理想时 A 的条件数会变大。M_flat.reshape(4, 4, H, W)把每个像素恢复成 4×4 矩阵第 0 维是矩阵行第 1 维是矩阵列后面两维是空间坐标。内存方面16 张 1200 万像素的 float32 图就有约 768MB直接把全部像素丢进lstsq很容易爆内存。我会改成按行或按块施解比如每次取 4096 个像素block 4096 M_block np.empty((4, 4, block, I_meas.shape[2])) for start in range(0, I_meas.shape[1], block): end min(start block, I_meas.shape[1]) M_block[..., :end - start, :] np.linalg.lstsq(A, I_flat[:, start:end], rcondNone)[0].reshape(4, 4, end - start, -1)逻辑说明把高度方向切块每块内所有列共享同一个 A 矩阵lstsq一次处理一个块耗时几乎线性增加但内存占用可控。M_block的最后一维仍然是宽度宽度太大时可以再继续切。4.4 偏振态选取与误差来源入射/分析 Stokes说明[1, 1, 0, 0]0° 线偏振能量集中在 S1[1, 0, 1, 0]45° 线偏振能量集中在 S2[1, 0, 0, 1]右旋圆偏振提供 S3 信息[1, -1, 0, 0]90° 线偏振与第一个状态互补选状态时要看 A 的条件数。条件数接近 1 表示方程组对噪声不敏感大于 100 时即使测量强度有 1% 噪声矩阵元素可能被放大到 100% 误差。用np.linalg.cond(A)在测量前验证是值得的。实际采集还要注意三件事一是暗场偏振相机长时间曝光后暗电流会叠加拍摄前盖住镜头拍一帧减去二是光强饱和穆勒矩阵反解是线性运算饱和像素会直接破坏方程组三是偏振片旋转误差哪怕 1° 的角度偏差也会在圆偏振项上放大最好用相机标定出的偏振效率而非理想模型。5. 进阶批量出图与 AoLP 相位缠绕的统计技巧5.1 一个可以直接跑的批量脚本把前面的函数串起来并加上合理的输出命名就能对一组偏振图批量生成 S0、DoLP、AoLP 图from pathlib import Path import cv2 import numpy as np import polanalyser as pa raw_root Path(/data/polar_raw) out_root Path(/data/polar_out) out_root.mkdir(exist_okTrue) for raw_path in sorted(raw_root.glob(*.jpg)): raw cv2.imread(str(raw_path), cv2.IMREAD_UNCHANGED) if raw.ndim 3: raw raw[..., 0] img pa.demosaicing(raw) stokes pa.calc_Stokes(img) dolp pa.calc_DoLP(stokes) aolp pa.calc_AoLP(stokes) prefix out_root / raw_path.stem cv2.imwrite(str(prefix) _S0.png, stokes[..., 0].astype(np.uint16)) cv2.imwrite(str(prefix) _DoLP.png, (np.clip(dolp, 0, 1) * 65535).astype(np.uint16)) cv2.imwrite(str(prefix) _AoLP.png, (aolp / np.pi * 180).astype(np.uint8))参数说明S0用 uint16 保存因为它是强度累加值8bit 会丢失暗部细节DoLP本身是 0 到 1 的小数乘 65535 转成 16bitAoLP只有 0 到 180 度用 8bit 足够但要注意 0 和 180 在显示上都是黑色不适合直接看细节需要配合伪彩色映射。5.2 统计 ROI 角度时先平均 S1/S2AoLP 是循环量0° 和 180° 物理上是同一个方向但直接对角度做平均会得出完全错误的结果。比如两个像素分别是 1° 和 179°直接平均是 90°而真实主轴方向接近 0°。正确做法是先对 S1、S2 取平均再反算角度roi (slice(200, 400), slice(300, 500)) S1_roi stokes[..., 1][roi] S2_roi stokes[..., 2][roi] mean_aolp 0.5 * np.arctan2(np.mean(S2_roi), np.mean(S1_roi)) mean_aolp_deg np.degrees(np.mod(mean_aolp, np.pi))逻辑说明S1/S2 是斯托克斯空间的线性分量平均后仍然是有物理意义的合成偏振方向。这个方法比调用scipy.stats.circmean更直接因为它没有损失偏振度信息。当区域内 DoLP 很低时S1/S2 的绝对值都很小此时mean_aolp会被噪声主导应该在 ROI 统计时把DoLP低于阈值的像素先排除。保存图片时的最后一个细节是AoLP只适合做定性展示做定量分析一定要保留原始 S0/S1/S2不要从伪彩色图反推角度。本文还有配套的精品资源点击获取
偏振图像分析实战：去马赛克、斯托克斯与穆勒矩阵解析
简介偏振图像分析工具用Python实现了完整的偏振图像处理流程面向图像处理、光学检测、遥感、生物医学成像和材料科学等方向的开发者涵盖去马赛克、斯托克斯向量计算与穆勒矩阵分析三大任务便于理解偏振光变化规律并完成图像复原、偏振参数提取和目标材质差异分析。压缩包共25个文件大小为1.49MB其中包含11个源代码文件、9张示例照片、3张示意图、1份说明文档和许可证文件。源代码包括偏振图像去马赛克、斯托克斯参数求解、穆勒矩阵计算等功能模块示例图片展示了不同偏振角度下拍摄的强度图以及经偏振角图、偏振度图处理后的效果图说明文档则对代码结构、依赖环境和使用方法进行了梳理。目前已有1516人学习/下载。借助这些代码读者可以快速搭建偏振图像分析实验环境参考示例数据理解偏振参数的计算过程并在此基础上修改或扩展算法适合希望在偏振成像方向快速上手的科研人员和工程师。1. 偏振图像分析工具的第一步从马赛克原图到物理量打开dragon_IMX250MZR_intensity.jpg这类偏振相机原始图时第一眼往往不是彩色照片而是一张灰扑扑的马赛克。普通相机拜耳阵列分的是 R/G/B 颜色偏振相机比如 Sony IMX250MZR则把 0°、45°、90°、135° 四个偏振方向的微偏振片铺在相邻 2×2 像素上原始图需要先拆成四个方向的强度通道才能继续算斯托克斯向量和穆勒矩阵。Polanalyser 就是这套流程的 Python 实现demosaicing.py负责去马赛克stokes.py计算 S0/S1/S2 与 DoLP/AoLPmueller.py处理多偏振态测量下的 16 元素矩阵反解。它适合用 Python 做偏振成像算法验证、光学检测标定或材料偏振响应研究的人不需要自己从头造轮子。2. 去马赛克把偏振拜耳阵列拆成四个独立强度2.1 偏振马赛克与普通拜耳马赛克的本质区别普通相机的拜耳阵列在 2×2 单元里放 R、G、G、B 三种滤色片去马赛克的本质是插值用周围像素估计每个位置的完整 RGB。偏振相机的 2×2 单元则放 0°、45°、90°、135° 四种微偏振片每个像素只能测到某一个方向的线偏振强度。它没有“颜色”需要恢复要的是把同一场景在四个方向的强度分别取出来。 如果你把原图直接当灰度图用2×2 内的偏振调制会一直混在像素里后续算 S1 I0 − I90 时每个像素拿到的都可能是相邻方向的串扰结果完全失真。这里还有一个很多人没注意的代价空间分辨率减半。四个方向的像素共享同一块靶面拆开后每个方向的有效图片是原图宽度和高度各一半。也就是说偏振相机的“原生分辨率”和普通相机同像素时并不等价。Polanalyser 的demosaicing.py不打算做超分辨它的职责就是把像素按微偏振片位置挑选出来排成四通道强度图。理解了这点下面读代码就不会被形状变化吓到。2.2 用 demosaicing 的核心逻辑拆分强度图下载的polanalyser-master里demosaicing.py对外暴露的函数在源码里通常就叫demosaicing。它的核心计算并不复杂如果你只想快速跑通下面这个切片版本足够import numpy as np import cv2 def demosaic_polar(raw, patternimx250mzr): 拆分偏振马赛克图返回 (H/2, W/2, 4) 的 float32 强度图。 通道顺序0度、45度、90度、135度。 if pattern imx250mzr: I0 raw[0::2, 0::2].astype(np.float32) I45 raw[0::2, 1::2].astype(np.float32) I90 raw[1::2, 0::2].astype(np.float32) I135 raw[1::2, 1::2].astype(np.float32) else: raise ValueError(unknown polarization mosaic pattern) return np.stack([I0, I45, I90, I135], axis-1) raw cv2.imread(dragon_IMX250MZR_intensity.jpg, cv2.IMREAD_UNCHANGED) if raw.ndim 3: raw raw[..., 0] # 某些原图会以三通道形式保存但只用了同一通道 img demosaic_polar(raw) print(img.shape, img.dtype) # (H/2, W/2, 4) float32逻辑说明raw[0::2, 0::2]取的是偶数行偶数列对应 0° 微偏振片的位置raw[0::2, 1::2]对应 45°raw[1::2, 0::2]对应 90°raw[1::2, 1::2]对应 135°。切片本身是视图astype(np.float32)会复制数据这样后续做减法时不会因为 uint16 溢出。pattern参数留给不同相机的微偏振片排布差异IMX250MZR 是常见的 0/45/90/135 按行排列。需要强调读取方式用cv2.imread(path, cv2.IMREAD_COLOR)会把单通道马赛克复制成三通道等值图再去切片就会让四个方向完全相同。所以IMREAD_UNCHANGED不是可选项是必须项。工具包里的examples脚本有的还会用tifffile或imageio读 12bit 原图效果类似。2.3 参数、dtype 与常见坑项目推荐做法原因读取原图cv2.imread(path, cv2.IMREAD_UNCHANGED)默认 flag 会把 12bit 压成 8bit破坏偏振强度线性度数据 dtypeuint16 读入分割后立即转 float32S1/S2 是差值会出现负数uint16 无法表示输入形状必须是 (H, W) 或 (H, W, 1)三通道等值马赛克会让四个方向完全相同通道顺序先统一成 0/45/90/135stokes.py的 I0−I90 和 I45−I135 依赖这个顺序输出形状(H/2, W/2, 4)每个方向共用一帧空间分辨率减半一个常见坑是相机厂商的 SDK 可能已经把数据排成了四通道或把 2×2 块改成“交错行”。这种情况下继续用上面的切片会得到棋盘状伪影。我一般会先对均匀光照下的平面拍一张图拆完后看四个通道均值是否接近如果某个通道明显偏亮或出现周期性条纹通常是行列偏移了一位把1::2和0::2对调即可。2.4 验证拆分是否正确means img.reshape(-1, 4).mean(axis0) print(means) # 均匀散射区域四个通道均值应接近逻辑说明对一片没有明显偏振特性的区域DoLP 接近 0四个方向的强度应该近似相等。任何一个通道显著偏离说明马赛克相位没对齐。这个检查只需要几行却是后面所有计算正确性的地基。如果直接跳到斯托克斯计算S1/S2 的错误会被可视化图掩盖成“颜色有点怪”很难定位是算法问题还是像素错位。3. 斯托克斯向量计算从四方向强度到 S0、S1、S23.1 斯托克斯参数的线性定义斯托克斯向量把光的偏振态写成四个实数对线偏振成像场景通常只关心前三个S0 是总光强S1 表示水平与垂直线偏振的差S2 表示 45° 与 135° 线偏振的差。只要拿到了去马赛克后的 I0、I45、I90、I135S0/S1/S2 就是简单的加减法def calc_stokes(img): 输入 (H/2, W/2, 4)返回 (H/2, W/2, 3)。 I0, I45, I90, I135 img[..., 0], img[..., 1], img[..., 2], img[..., 3] S0 I0 I90 S1 I0 - I90 S2 I45 - I135 return np.stack([S0, S1, S2], axis-1)逻辑说明S0 用一对正交方向相加是因为 0° 和 90° 强度之和等于总强度S1、S2 则是两组正交方向的差分。源码包里的stokes.py基本就是这几行只是额外处理了输入通道顺序和 NaN。要注意坐标系定义有些教材把 S1 写成 I90 − I0这样可视化时 S1 的符号会左右翻转AoLP 也会偏移 90°。同一个包里demosaicing.py和stokes.py必须保持同一套约定我一般会在一开始写注释固定通道顺序。3.2 DoLP 与 AoLP归一化和弧度处理只拿到 S0/S1/S2 还不够实际分析经常要偏振度 DoLP 和偏振角 AoLP。DoLP 表示偏振光占总强度的比例范围 0 到 1AoLP 表示偏振主轴方向物理上只有 0 到 π 的周期。实现如下def calc_DoLP(stokes, eps1e-6): S0, S1, S2 stokes[..., 0], stokes[..., 1], stokes[..., 2] return np.sqrt(S1 ** 2 S2 ** 2) / np.maximum(S0, eps) def calc_AoLP(stokes): S1, S2 stokes[..., 1], stokes[..., 2] aolp 0.5 * np.arctan2(S2, S1) return np.mod(aolp, np.pi) # 主轴在 [0, π) 内参数说明eps1e-6是为了避免 S0 接近 0 时除零np.mod(aolp, np.pi)必须加因为0.5 * arctan2返回范围是 −π/2 到 π/2而偏振角是 90° 周期直接转成 0 到 π 的循环量才方便做伪彩色和统计。如果想要角度显示乘180 / np.pi注意是 0 到 180 度不是 0 到 360 度。3.3 视觉化examples/visualize_AoLP.py 的思路示例目录里的visualize_AoLP.py不是简单地显示灰度图而是把 AoLP 映射到色环颜色代表偏振方向亮度和饱和度分别表示光强和偏振度。这样一张图能同时看出三个维度。核心代码如下def visualize_stokes(stokes): S0 stokes[..., 0] DoLP calc_DoLP(stokes) AoLP calc_AoLP(stokes) hsv np.zeros((S0.shape[0], S0.shape[1], 3), dtypenp.uint8) hsv[..., 0] (AoLP / np.pi * 179).astype(np.uint8) hsv[..., 1] (np.clip(DoLP, 0, 1) * 255).astype(np.uint8) hsv[..., 2] (S0 / np.max(S0) * 255).astype(np.uint8) return cv2.cvtColor(hsv, cv2.COLOR_HSV2BGR)逻辑说明OpenCV 的 8bit HSV 模式中 Hue 范围是 0 到 179所以把 AoLP 的 0 到 π 映射到 0 到 179。Saturation 放 DoLPValue 放归一化光强 S0。dragon_IMX250MZR_AoLP_color.jpg这类示例图就是这么生成的。jpge保存时质量选高一些否则色环边缘会出现条带。3.4 用 S0 做自检斯托克斯计算完成后我会先做一个最小验证把 S0 单独显示出来它应当和原图经过 2×2 池化后的结果很像只是尺寸减半。如果 S0 出现棋盘格说明去马赛克的行列错位如果 S0 基本正确但 S1/S2 一片杂乱则要检查暗场和坏点。下面的代码可以替代人眼判断S0 stokes[..., 0] S0_expected img[..., 0] img[..., 2] # I0 I90 print(np.max(np.abs(S0 - S0_expected)))逻辑说明S0_expected直接从去马赛克结果里相加理论上必须与calc_stokes输出的 S0 完全一致。这个断言能在早期发现通道顺序写反、数组被原地修改等问题。注意这里不是和img.mean(axis-1)对比因为mean除以了 4。4. 穆勒矩阵计算从一组偏振态测量反解物体的偏振响应4.1 为什么一张图算不出穆勒矩阵斯托克斯向量描述光的偏振状态穆勒矩阵描述物体把入射偏振态变成出射偏振态的 4×4 线性变换。每个像素在不同入射偏振态下会给出不同出射强度单张偏振图只对应一次入射最多得到一个方程而一个未知物体有 16 个矩阵元素要解。因此至少需要四个线性无关的入射 Stokes 向量配合四个线性无关的分析态测量。Polanalyser 的mueller.py做的事情就是把一组测量强度图按最小二乘反解成每个像素的 4×4 矩阵。这里要区分两个概念去马赛克和斯托克斯计算都是“单帧处理”穆勒矩阵需要多帧。最简实验装置是入射侧放一个可旋转偏振片产生偏振态 PSG出射侧放另一个偏振片分析 PSA中间放被测样品。你下载的源码包里examples/mueller.py就是这个流程的演示。4.2 测量方程与线性反解理想线偏振片在角度 θ 处的 Stokes 向量是(1, cos2θ, sin2θ, 0)。如果入射 Stokes 是s_in样品矩阵是 M出射经过分析器后的强度满足I a_out^T · M · s_in把 M 按行优先展平成 16 维向量m则每个测量就是一个行向量kron(a_out, s_in)与m的内积。16 组不同偏振态组成 16×16 矩阵 A强度图展平成(16, N)用最小二乘一次解出所有像素的矩阵元素。以下是通用的方程构建函数import numpy as np def build_mueller_equations(S_in_list, A_out_list): 根据入射/分析斯托克斯向量构造测量矩阵 A。 A [] for s_in in S_in_list: for a_out in A_out_list: A.append(np.kron(a_out, s_in)) return np.array(A) states [ [1, 1, 0, 0], # 0° 线偏振 [1, 0, 1, 0], # 45° 线偏振 [1, 0, 0, 1], # 右旋圆偏振 [1, -1, 0, 0], # 90° 线偏振 ] A build_mueller_equations(states, states) print(condition number:, np.linalg.cond(A))逻辑说明states同时作为入射和分析器状态覆盖了线偏振和圆偏振保证 A 可逆。np.kron(a_out, s_in)的排列顺序与M.flatten()的展开方式一致kron中参数顺序不能反。如果只用 0/45/90/135 线偏振A 会奇异因为所有输入输出状态都落在 S30 的子空间穆勒矩阵的第三行和第三列永远解不出来。真实测量至少需要一组四分之一波片产生圆偏振态。4.3 像素级反解与内存控制拿到 A 矩阵后按测量顺序把强度图堆成(16, H, W)的数组再展平求解# I_meas 形状 (16, H, W)顺序必须和 build_mueller_equations 的循环顺序一致 I_flat I_meas.reshape(16, -1) M_flat np.linalg.lstsq(A, I_flat, rcondNone)[0] M M_flat.reshape(4, 4, H, W) # 查看某个区域的平均穆勒矩阵 print(M[..., :3, :3].mean(axis(2, 3)))参数说明I_flat是 16 行、N 列的矩阵每列对应一个像素的 16 次测量。lstsq在这种情况下比inv(A) I_flat更稳因为偏振片不理想时 A 的条件数会变大。M_flat.reshape(4, 4, H, W)把每个像素恢复成 4×4 矩阵第 0 维是矩阵行第 1 维是矩阵列后面两维是空间坐标。内存方面16 张 1200 万像素的 float32 图就有约 768MB直接把全部像素丢进lstsq很容易爆内存。我会改成按行或按块施解比如每次取 4096 个像素block 4096 M_block np.empty((4, 4, block, I_meas.shape[2])) for start in range(0, I_meas.shape[1], block): end min(start block, I_meas.shape[1]) M_block[..., :end - start, :] np.linalg.lstsq(A, I_flat[:, start:end], rcondNone)[0].reshape(4, 4, end - start, -1)逻辑说明把高度方向切块每块内所有列共享同一个 A 矩阵lstsq一次处理一个块耗时几乎线性增加但内存占用可控。M_block的最后一维仍然是宽度宽度太大时可以再继续切。4.4 偏振态选取与误差来源入射/分析 Stokes说明[1, 1, 0, 0]0° 线偏振能量集中在 S1[1, 0, 1, 0]45° 线偏振能量集中在 S2[1, 0, 0, 1]右旋圆偏振提供 S3 信息[1, -1, 0, 0]90° 线偏振与第一个状态互补选状态时要看 A 的条件数。条件数接近 1 表示方程组对噪声不敏感大于 100 时即使测量强度有 1% 噪声矩阵元素可能被放大到 100% 误差。用np.linalg.cond(A)在测量前验证是值得的。实际采集还要注意三件事一是暗场偏振相机长时间曝光后暗电流会叠加拍摄前盖住镜头拍一帧减去二是光强饱和穆勒矩阵反解是线性运算饱和像素会直接破坏方程组三是偏振片旋转误差哪怕 1° 的角度偏差也会在圆偏振项上放大最好用相机标定出的偏振效率而非理想模型。5. 进阶批量出图与 AoLP 相位缠绕的统计技巧5.1 一个可以直接跑的批量脚本把前面的函数串起来并加上合理的输出命名就能对一组偏振图批量生成 S0、DoLP、AoLP 图from pathlib import Path import cv2 import numpy as np import polanalyser as pa raw_root Path(/data/polar_raw) out_root Path(/data/polar_out) out_root.mkdir(exist_okTrue) for raw_path in sorted(raw_root.glob(*.jpg)): raw cv2.imread(str(raw_path), cv2.IMREAD_UNCHANGED) if raw.ndim 3: raw raw[..., 0] img pa.demosaicing(raw) stokes pa.calc_Stokes(img) dolp pa.calc_DoLP(stokes) aolp pa.calc_AoLP(stokes) prefix out_root / raw_path.stem cv2.imwrite(str(prefix) _S0.png, stokes[..., 0].astype(np.uint16)) cv2.imwrite(str(prefix) _DoLP.png, (np.clip(dolp, 0, 1) * 65535).astype(np.uint16)) cv2.imwrite(str(prefix) _AoLP.png, (aolp / np.pi * 180).astype(np.uint8))参数说明S0用 uint16 保存因为它是强度累加值8bit 会丢失暗部细节DoLP本身是 0 到 1 的小数乘 65535 转成 16bitAoLP只有 0 到 180 度用 8bit 足够但要注意 0 和 180 在显示上都是黑色不适合直接看细节需要配合伪彩色映射。5.2 统计 ROI 角度时先平均 S1/S2AoLP 是循环量0° 和 180° 物理上是同一个方向但直接对角度做平均会得出完全错误的结果。比如两个像素分别是 1° 和 179°直接平均是 90°而真实主轴方向接近 0°。正确做法是先对 S1、S2 取平均再反算角度roi (slice(200, 400), slice(300, 500)) S1_roi stokes[..., 1][roi] S2_roi stokes[..., 2][roi] mean_aolp 0.5 * np.arctan2(np.mean(S2_roi), np.mean(S1_roi)) mean_aolp_deg np.degrees(np.mod(mean_aolp, np.pi))逻辑说明S1/S2 是斯托克斯空间的线性分量平均后仍然是有物理意义的合成偏振方向。这个方法比调用scipy.stats.circmean更直接因为它没有损失偏振度信息。当区域内 DoLP 很低时S1/S2 的绝对值都很小此时mean_aolp会被噪声主导应该在 ROI 统计时把DoLP低于阈值的像素先排除。保存图片时的最后一个细节是AoLP只适合做定性展示做定量分析一定要保留原始 S0/S1/S2不要从伪彩色图反推角度。本文还有配套的精品资源点击获取
