四偏振图像处理：用MATLAB还原偏振角与偏振度图像

简介：面向光学成像与机器视觉学习者的偏振图像分析资源包，覆盖偏振椭圆、偏振角、偏振角图像、四偏振图像及椭圆偏振率等核心概念，包含从基础理论到MATLAB实现的可运行示例，便于快速搭建偏振信息处理流程。压缩包共2个文件，包含1个.m脚本与1个.bmp图像，整体仅268KB，轻量便于查看和调试。资源已有859人学习/下载，适合研究生、算法工程师用于偏振成像基础学习或算法验证。脚本可完成不同偏振角度图像的裁剪与合成，并计算输出偏振角与偏振度图像；配套示例图直观呈现像素级偏振分布，帮助理解椭圆偏振率对线性偏振偏离程度的量化，以及反射、折射等光学效应下的偏振变化，为光学检测、遥感分析、对比度增强等应用提供基础。整体代码结构简洁，便于在此基础上二次开发与算法改进。

1. 拆开pianzhen.zip：一份四偏振图像处理脚本的起点

拿到这个压缩包时，里面只有两个文件：pianzhen.m和一张偏振度图像.bmp。表面看像课程作业，实际处理对象并不简单：四偏振图像指的是同一场景在 0°、45°、90°、135° 四个检偏偏振方向下拍到的四张灰度图，而pianzhen.m要做的是从这四张图里逐像素解算出偏振角图像和偏振度图像。这个思路和 Sony IMX250MZR 这类分焦平面偏振相机机内计算逻辑一致，只是相机固件把结果直接输出了，这里要自己写。适合刚接触偏振成像、又想弄清楚偏振角、偏振椭圆参数怎么从普通强度灰度图里推导出来的工程师。

2. 偏振椭圆与Stokes参数：四偏振图像还原偏振角的前提

2.1 四张强度图能组合出的偏振信息

偏振光电场矢量在垂直于传播方向的平面内划出轨迹，偏振椭圆描述了电场矢量末端划过的形状，椭圆长轴的方位角就是后续输出的偏振角。普通相机只记录能量，把椭圆轨迹压成灰度，所以只看单张 0° 图像根本分不清场景表面是被吸收还是偏振方向在变化。这就是为什么处理时要同时拿到多个检偏方向下的光强，利用检偏器对不同振动方向的投影关系反推椭圆的形状和朝向。

检偏器透光轴方向与水平方向夹角为 α 时，出射光强可以写成 (I(\alpha)=\frac{S_0}{2}(1 + \frac{S_1}{S_0}\cos(2\alpha) + \frac{S_2}{S_0}\sin(2\alpha)))。这里 S0、S1、S2 是 Stokes 参数前三项。当采集了 0°、90° 两张图像时，S0 和 S1 可以直接差分得到；再补上 45°、135° 两张图像，S2 也就能确定。pianzhen.m的核心计算必然落在这三条公式上。

实际写脚本时要注意图像读进去的类型，下面是一段最基础的 Stokes 参数计算代码：

I0 = double(imread('0deg.bmp')); % 0° 通道，转 double 避免后续减法溢出 I45 = double(imread('45deg.bmp')); % 45° 通道 I90 = double(imread('90deg.bmp')); % 90° 通道 I135 = double(imread('135deg.bmp')); % 135° 通道 S0 = I0 + I90; % 总强度项 S1 = I0 - I90; % 水平分量与垂直分量的差 S2 = I45 - I135; % 对角方向分量的差

imread对 8bit BMP 默认返回 uint8 矩阵，uint8 相减遇到负值会截断成 0，所以必须先转 double。S0 描述光场总强度，在数值上等于 0° 和 90° 图像之和；S1 为正说明水平偏振分量占优，为负则垂直偏振分量占优；S2 对应 45° 与 135° 两个对角方向的差异。这里有容易混淆的习惯：有些教材把 S0 写成 ((I0+I90)/2)，这在偏振度计算里不会改变百分比结果，但如果你把 S1/S0 拿去和相机 SDK 输出的 DOLP 对比，会发现数值差了一倍。我一般直接让 S0 = I0 + I90，并在脚本注释里标明白，避免后面对接其他 SDK 时产生歧义。

2.2 从 S1、S2 推导偏振椭圆与偏振角图像

偏振椭圆的倾斜角也就是偏振角 θ，用半角公式恢复：

[\theta = \frac{1}{2}\arctan_2(S2, S1)]

由于 θ 和 θ+π 描述的椭圆取向相同，偏振角天然拥有 π 周期性。MATLAB 里最常见的写法是：

AOP = 0.5 * atan2(S2, S1); % 结果范围 [-pi/2, pi/2] AOP = mod(AOP, pi); % 映射到 [0, pi]

atan2接收两个参数而不是atan(S2./S1)的原因是后者在 S1 接近 0 时会产生除零和大角度跳变，而且无法区分 45° 与 225° 象限。用atan2得到的角度在 [-90°, 90°] 区间，再mod(AOP, pi)映射到 [0°, 180°)，这样输出图像不会出现因正负角跳变导致的黑白撕裂。手算验证：若 0° 通道最亮、90° 通道最暗，即 S1 > 0、S2 = 0，AOP = 0°；若 45° 通道最亮，则 AOP 大约为 45°。

偏振椭圆中电场矢量旋转的幅度，还需要引入椭圆率角：

[\chi = \frac{1}{2}\arcsin(S3 / S0)]

S3 代表圆偏振分量。四偏振图像里只有 I0、I45、I90、I135 时，S3 没有信息来源，所以pianzhen.m如果直接给出椭圆偏振相关结果，多半是把 S3 置为 0，或者像我一样先输出线偏振部分。要完整重建偏振椭圆，需要额外采集左右旋圆偏振通道，这个边界在第 4 章继续展开。

注意：偏振角为 0° 和 180° 在物理上等价，显示偏振角图像时不要用 0 到 360 度映射，否则同一个椭圆取向会被分成两种颜色，伪彩色图看起来像发生了角度突变。

3. 用pianzhen.m把偏振角图像拼出来：四通道输入约定与合成流程

3.1 解压后的输入约定需要先确认

解开pianzhen.zip后，压缩包里并没有原始的四张 0°、45°、90°、135° 图像，只留下pianzhen.m和一张处理后 BMP。这说明脚本的输入要么是从外部读取四个文件，要么是变量里已经放好四偏振图像。常见的做法是用 cell 数组遍历文件名，再按角度顺序读入。顺序一旦写反，S1、S2 的符号会跟着反，最终偏振角图像会整体旋转 45° 或 90°，肉眼很难直接看出来。

files = {'0deg.bmp', '45deg.bmp', '90deg.bmp', '135deg.bmp'}; imgs = cell(1, 4); for k = 1:4 imgs{k} = double(imread(files{k})); if size(imgs{k}, 3) == 3 imgs{k} = rgb2gray(imgs{k}); % 个别相机输出三通道 RGB，但偏振强度本质是灰度 end end I0 = imgs{1}; % 0° I45 = imgs{2}; % 45° I90 = imgs{3}; % 90° I135 = imgs{4}; % 135°

这个循环里对三通道图像做了降维处理，因为偏振片后面接的相机如果是彩色传感器，解码后会出现 RGB 三分量，直接使用会把 Bayer 马赛克当成有效信号。rgb2gray在这里仅用于排除色彩干扰，真正高质量偏振数据采集通常已经是黑白相机直接输出灰度，不需要这一步。文件名里我刻意保留了deg角度后缀，实际工程里命名可能是cam0.tif、cam45.tif，建议在读取后立刻打印尺寸和角度顺序，确认四张图的空间分辨率完全一致。

3.2 合成四偏振图像并计算偏振角图像

偏振处理中常用到一个概念：四偏振图像不只是四张独立的灰度图，而是被组合成一个 (H \times W \times 4) 的数据立方体。这样后续的 ROI 裁剪、滤波、逐像素计算都共享同一套空间索引。pianzhen.m里大概率就有一截类似下面的代码，先把四通道三维叠起来，再裁掉传感器边缘的阴影区域。

quad = cat(3, I0, I45, I90, I135); % 四偏振图像数据立方体 ROI = quad(100:357, 200:511, :); % 裁掉边缘暗区，尺寸由实际画面决定 I0c = ROI(:, :, 1); I45c = ROI(:, :, 2); I90c = ROI(:, :, 3); I135c = ROI(:, :, 4); S0 = I0c + I90c; S1 = I0c - I90c; S2 = I45c - I135c; AOP = 0.5 * atan2(S2, S1); AOP = mod(AOP, pi); % [0, 180°) AOP_img = uint8(AOP / pi * 255); % 映射到 0~255 灰度

cat(3, ...)沿第三维拼接，第三维索引 1 到 4 分别对应检偏角度，而不是 0° 到 135° 的数值本身。ROI 参数 100:357、200:511 是拍工业样品时常用的经验裁剪范围，用来去掉因镜头边缘照度不均造成的偏振角漂移。映射到 0~255 的AOP / pi * 255把 180° 压到 255，丢失了角度刻度但方便直接存 BMP 预览。

保存时尽量保留一份弧度原始数据，因为灰度图已经丢失了角度量纲，后续想统计某个区域平均偏振角时必须重新从 mat 文件读取。伪彩色建议用 hsv 色带，它首尾相连，恰好对应偏振角的 0° 和 180° 同值关系，避免 viridis 这类两端色差给人错误印象。

3.3 四通道配准错位的快速判断

分焦平面偏振相机在一个传感器上做 2×2 像素马赛克，实际输出的四偏振图像并不是四次独立曝光，而是空间交错采样。读出时要按 2×2 单元拆成四张子图，再对齐到同一像素网格；pianzhen.m如果处理的是旋转检偏器方案，则不存在马赛克问题。判断脚本有没有配准错位，一个快速办法是看 S0 图像上是否出现规则网格条纹。S0 本身是总强度，理论上应该接近普通灰度图，不该有棋盘状高频分量。一旦看到这类条纹，就说明四张子图之间有亚像素偏移，需要重新做相位相关配准再继续算偏振角图像，否则偏振角会在物体边缘出现规律性误差。

4. 偏振度图像与椭圆偏振率：线性通道能算到哪一步

4.1 从偏振椭圆投影关系计算偏振度图像

压缩包里那张偏振度图像.bmp体现了整个处理流程的最终目标之一。偏振度衡量每个像素位置光场中偏振成分占整体强度的比例，对只有四个线偏振通道的输入，可以直接算的是线偏振度 DOLP：

[ \mathrm{DOLP} = \frac{\sqrt{S1^2 + S2^2}}{S0} ]

这个公式来自偏振椭圆长轴与短轴能量之差，长轴方向在 0° 和 90° 通道之间体现，短轴方向在 45° 和 135° 通道之间体现。写成 MATLAB 一行就是：

DOLP = sqrt(S1.^2 + S2.^2) ./ (S0 + eps); % eps 防止除零 DOLP(DOLP > 1) = 1; % 物理上限为 1 DOLP(DOLP < 0) = 0; % 数值异常截断

eps是 MATLAB 内置机器精度常量，加到 S0 上能避免黑色区域除零得到 NaN。截断操作主要是兜底，因为输入图像若有坏像素，S0 接近 0 时 DOLP 可能被放大成几十。真正的难点在噪声：暗光环境下 S0 较小，DOLP 噪声会被平方项放大，导致偏振度图像看起来满是雪花。我一般会在计算后加一个强度蒙版，只对信号足够的像素保留偏振度结果：

mask = S0 > 30; % 8bit 图像，阈值常用 20~50，动态范围不同要重新标定 DOLP = DOLP .* mask; % 低照度像素直接置 0

阈值 30 不是固定值。如果相机是 12bit 输出，动态范围从 0 到 4095，阈值要放到 200 以上；反之 8bit 输出下取 30 左右比较合理。判断标准是蒙版不要吃掉暗部物体轮廓，也不要留下纯噪声区域。偏振度图像通常被用来识别材质边界，金属表面镜面反射的 DOLP 明显高于漫反射区域，所以这张图比普通灰度图更容易突出缺陷边缘。

4.2 椭圆偏振率图像的真实起点

项目正文里提到的椭圆偏振率，量化的是偏振态偏离线偏振的程度。用 Stokes 第三项 S3 表示就是对 (\chi = \arcsin(S3/S0)/2) 的求解。但这里有个绕不开的信息边界：0°、45°、90°、135° 四个检偏方向都是线偏振片，不可能测出圆偏振分量，S3 在数学上无解。如果pianzhen.m强行输出椭圆偏振率图像，常见做法是先创建一个与 S0 同尺寸的零矩阵：

S3 = zeros(size(S0)); % 四线偏振通道无法获得 S3 chi = 0.5 * asin(min(max(S3 ./ (S0 + eps), -1), 1)); % 椭圆率角

这时代码能跑通，但输出图像全黑，因为线偏振光的椭圆率角就是 0。真正要得到椭圆偏振率图像，必须改变采集方案：在镜头前加一个可旋转四分之一波片，或者使用带有圆偏振像素的专用偏振相机。拿旋转波片方案举例，需要额外拍一张经过右旋圆偏振通道的强度图 IR：

IL = 0.5 * (I0 + I90); % 近似总强度的一半 S3_est = 2 * IR - IL; % 右旋圆偏振强度与总圆偏振贡献的差 chi = 0.5 * asin(min(max(S3_est ./ (S0 + eps), -1), 1)); % 椭圆率角

2 * IR - IL的来源是 S3 定义为右旋与左旋圆偏振强度之差，仅采集 IR 时用总强度一半做基准近似，这是工业现场常用的简化处理。得到 chi 后，偏振椭圆的短轴与长轴之比就是 (\tan(\chi))。若后续把 AOP 和 chi 组合进同一个 HSV 图像，H 通道放偏振角、V 通道放椭圆率，就能在一张图里同时表达偏振椭圆的方向和圆度，这也是很多商业偏振相机的可视化格式。

5. 进阶：用合成偏振靶标反向验证pianzhen.m的输出

偏振成像调试最大的坑是缺少真值。从真实场景拍来的四偏振图像，你不知道每个像素真实的偏振角是多少，就算pianzhen.m算出一个怪异的偏振角图像，也很难判断是算法问题还是场景本来如此。解决办法是用合成偏振靶标：先构造一个偏振角已知的图像，再反推出四偏振图像，把结果送进同一套 Stokes 解算，比较恢复出的偏振角和输入理论值是否一致。

合成靶标的数学关系就是第 2 章的光强公式 (I(\alpha)=\frac{S0}{2}(1 + \mathrm{DOLP}\cos(2(\theta-\alpha))))。下面生成一幅偏振角从 20° 渐变到 70° 的测试图：

[x, y] = meshgrid(1:256, 1:256); theta_map = deg2rad(20 + 50 * (x + y) / 512); % 对角方向 20° 到 70° 渐变 S0v = 0.8; DOLP = 0.6; I0 = S0v/2 * (1 + DOLP * cos(2 * (theta_map - 0))); I45 = S0v/2 * (1 + DOLP * cos(2 * (theta_map - deg2rad(45)))); I90 = S0v/2 * (1 + DOLP * cos(2 * (theta_map - deg2rad(90)))); I135 = S0v/2 * (1 + DOLP * cos(2 * (theta_map - deg2rad(135))));

meshgrid生成空间坐标，theta_map是每个像素的偏振椭圆长轴角度，随 x 和 y 线性变化。四张合成图之间的差异完全由 cos 项决定，模拟了理想检偏器输出。把它存成 BMP 后跑pianzhen.m的计算流程，再对比恢复出的 AOP 与theta_map的差值：

est = 0.5 * atan2(I45 - I135, I0 - I90); % 用同一组公式恢复 est = mod(est, pi); err = rad2deg(est(:) - theta_map(:)); err = abs(mod(err + pi/2, pi) - pi/2); % 处理 180° 周期折叠 fprintf('最大误差 %.4f 度\n', max(err(:)));

这段误差统计里mod再次出现，是因为偏振角 0° 与 180° 等价，直接相减会得到接近 180 的假误差。理论上合成数据无噪声时最大误差应该小于1e-3度；如果超过 0.1 度，优先检查四张图像是否在读写过程中被做了数值压缩或尺寸对齐。

还有一个更快的通道错位检查：把恢复流程里的 I0 和 I90 交换后再算一次，偏振角图像应整体增加 90°；把 I45 与 I135 交换，偏振角分布会沿 45° 方向镜像。用合成靶标跑一遍这两组交换，再和原始输出对比，就能确认pianzhen.m里四个通道的排列顺序没有接反。

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