
简介面向光学成像与机器视觉学习者的偏振图像分析资源包覆盖偏振椭圆、偏振角、偏振角图像、四偏振图像及椭圆偏振率等核心概念包含从基础理论到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的核心计算必然落在这三条公式上。检偏角度图像变量参与组合物理含义0°I0S0 S1水平方向线偏振的投影强度90°I90S0 - S1垂直方向线偏振的投影强度45°I45S0 S245° 方向线偏振的投影强度135°I135S0 - S2−45° 方向线偏振的投影强度实际写脚本时要注意图像读进去的类型下面是一段最基础的 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 写成 ((I0I90)/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 0AOP 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 预览。输出方式指令适用场景灰度图imwrite(uint8(AOP/pi*255), AOP.bmp)快速保存原始角度映射伪彩色imagesc(AOP); colormap(hsv); colorbar观察偏振角空间分布弧度数据save(AOP.mat, AOP)后续定量分析继续使用角度制AOP_deg rad2deg(AOP);与其他传感器标定数据对比保存时尽量保留一份弧度原始数据因为灰度图已经丢失了角度量纲后续想统计某个区域平均偏振角时必须重新从 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。真正要得到椭圆偏振率图像必须改变采集方案在镜头前加一个可旋转四分之一波片或者使用带有圆偏振像素的专用偏振相机。拿旋转波片方案举例需要额外拍一张经过右旋圆偏振通道的强度图 IRIL 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里四个通道的排列顺序没有接反。本文还有配套的精品资源点击获取