MATLAB偏振成像实战：斯托克斯参数计算与图像处理方法

简介：本资源是一套基于偏振光学原理的MATLAB图像处理实践方案，面向光学工程、机器视觉及计算成像方向的本科生与研究生，解决从原始偏振图像获取斯托克斯参量、偏振度（DOLP）与偏振角（AOP）图像的核心计算问题。压缩包共9个文件，含8幅BMP格式原始偏振图像（对应偏振片在0°、45°、90°等三组角度采集的hv1/hh1/hp1/I等分量图）及1个核心MATLAB脚本pianzhenxinxiyuanshi.m，完整实现斯托克斯四维分量S0–S3的逐像素重建与可视化。资源大小仅540KB，轻量易部署，代码结构清晰、注释充分，可直接运行复现偏振度图像（DOLP.bmp）、偏振角图像（AOP.bmp）及各斯托克斯分量图（U.bmp等）。目前已有1356人学习下载，适用于课程实验、毕业设计中偏振成像算法验证与教学演示，提供从数据采集逻辑到数学推导再到图像输出的全流程支撑。 先把结论摆在前面：把一片线偏振片装在相机镜头前，依次旋转到 0°、60°、120°，固定相机和场景，分别记录三幅偏振图像；回到 MATLAB 里用矩阵运算把三幅图换算成斯托克斯分量图像，再由斯托克斯分量算出偏振度图像和偏振角图像。这套流程看着简单，但里面真正决定结果好坏的不是最后那几行代码，而是前面的角度校准、数据处理和公式选型。这篇文章按我实际调通的流程来写，把原理、MATLAB 代码、可视化方法、常见坑一次讲清楚。适合正在做图像处理课程设计、偏振成像入门实验，或者想用 MATLAB 搭建斯托克斯成像分析工具的朋友。你不需要太深的光学基础，只要会拍三幅图、能跑 MATLAB 脚本，就能完整复现。

1. 偏振成像到底在算什么：从三幅光强图到三个重要结果

1.1 为什么用三幅图，普通照片不行吗

普通相机只能记录光强，按像素输出一个灰度值，光的偏振特性在这个过程中被完全丢掉了。但你如果盯着天空、水面、玻璃、塑料表面看，就会发现不同材质的反射光在偏振特性上差别非常明显。偏振成像就是把这些被普通相机扔掉的信息重新捡回来。

要测量偏振态，最简单的办法就是在相机前放一片线偏振片，旋转到不同角度，记录不同透光方向下的光强。你可能想问：转两个角度够不够？不够。只拍 0° 和 90°，能得到两个方向上的光强差，但没法区分 45° 方向上的偏振分量，也就是说少了描述偏振方向的关键信息。至少三个角度才能解出线偏振的三个基本参数，这就是标题里“旋转三个不同角度得到三幅偏振图像”的根本原因。

从数学上看，任意透光方向 θ 下的光强可以写成：

I(θ) = 0.5 × [S0 + S1 × cos(2θ) + S2 × sin(2θ)]

这是一个只有三个未知量 S0、S1、S2 的方程。三个不同的 θ 给出三个方程，刚好能唯一解出这三个未知量。偏振片旋转角度越多当然越好，但三幅图是满足“可解”的最低要求，也是最经典的实验配置。

1.2 斯托克斯参数 I、Q、U 的物理含义

斯托克斯矢量在完全偏振光描述里通常有四个分量，但只靠偏振片旋转，拍不到圆偏振分量，所以这里只关注前三个线偏振相关分量：

• S0（I）：总光强，和普通相机拍到的灰度图含义一致。
• S1（Q）：0° 方向线偏振分量与 90° 方向线偏振分量的差值，描述偏振方向沿着水平还是垂直。
• S2（U）：45° 方向线偏振分量与 135° 方向线偏振分量的差值，描述偏振方向沿着 45° 还是 135°。

可以这样理解：S0 是光的“总量”，S1 和 S2 是光振动方向在平面两个坐标轴上的投影。有了这两个分量，就能知道线偏振光偏得有多厉害、朝哪个方向偏。

很多材料表面反射后，光的偏振态会发生明显变化。比如光滑非金属表面反射的光，偏振度往往较高；粗糙表面散射后偏振度很低。斯托克斯分量图像的意义在于：它把光强和偏振信息拆开存放，后续做去雾、材质分类、目标检测时可以直接利用这些物理量。

1.3 偏振度、偏振角与斯托克斯分量的关系

偏振度（Degree of Linear Polarization，DoLP）描述光波中偏振部分所占的比例，定义为线偏振分量总强度与总光强之比：

DoLP = sqrt(S1² + S2²) / S0

偏振度是一个 0 到 1 之间的无量纲量。DoLP = 0 表示完全非偏振光，比如很多漫反射光；DoLP = 1 表示完全线偏振光，比如某些镜面反射或经过偏振片后的出射光。在实际图像中，大部分像素的偏振度都在 0.1 到 0.6 之间，极少接近 1。

偏振角（Angle of Linear Polarization，AoLP）描述偏振方向的朝向角度，定义为：

AoLP = 0.5 × atan2(S2, S1)

这里用 atan2 而不是 atan，是因为 atan2 能根据 S1、S2 的符号自动判断角度所在的象限。AoLP 单位是弧度，转成角度要乘以 180/π。由于线偏振方向和它的反向方向无法区分，所以 AoLP 的范围被压缩到 -90° 到 90°，或者用 mod 转到 0° 到 180°。

到这里，核心目标就很清晰了：从三幅不同偏振角度下的光强图像，估算 S0、S1、S2，再算出 DoLP 和 AoLP。后面所有 MATLAB 代码都是围绕这个公式展开的。

2. 方案选型：三个旋转角度怎么选，MATLAB 整体流程怎么搭

2.1 0°/60°/120° 和 0°/45°/90° 两种经典方案的取舍

旋转偏振片时，最常用的两组角度是 0°、60°、120° 和 0°、45°、90°。两套方案都能算出前三个斯托克斯分量，但在误差特性和公式复杂度上有点区别。

先看 0°、60°、120°。这三个角度在 0° 到 180° 范围内等间隔分布，公式解出来也比较对称：

S0 = (2/3) × (I0 + I60 + I120)
S1 = (2/3) × (2I0 - I60 - I120)
S2 = (2/√3) × (I60 - I120)

这套公式对三幅图的噪声放大相对均匀，不会出现某个方向特别敏感的情况。如果你有一台带角度刻度盘的精密偏振片旋转架，我推荐优先用这组角度。

再看 0°、45°、90°。这三组角度更容易手动摆放，公式也简洁：

S0 = I0 + I90
S1 = I0 - I90
S2 = 2I45 - I0 - I90

但这里有个容易忽略的问题：S2 的计算直接依赖 I45 的绝对值稳定性，如果 I45 有一点点系统误差，S2 的偏差会放大，偏振角图像容易出现颜色跳变。另外，0° 和 90° 的光强差异比较大时，S1 的动态范围会很宽，可视化时要注意范围压缩。

两套方案都能用，关键是别混用公式。如果你用 0°、60°、120° 拍摄，却套用 0°、45°、90° 的公式，结果会错得很离谱。代码里我会用一个变量 angle_mode 来控制，避免这种低级错误。

为什么不推荐只拍 0°、90° 两张图？因为两张图只能求出 S0 和 S1，无法得到 S2。这意味着你隐含假设 S2 = 0，也就是场景里不存在 45° 方向的偏振分量。自然场景里这种假设基本不成立，所以两张图算出的偏振度会被低估，偏振方向也不正确。

2.2 完整处理流程与代码骨架

整个流程可以分成五步：拍摄准备、图像导入与预处理、偏振参数计算、结果可视化、结果保存。

拍摄准备阶段要注意三件事：相机和场景固定不动；偏振片旋转时不要碰到相机；保证光源在拍摄期间不闪烁。如果使用手持偏振片，轻微旋转可能造成微小位移，最后需要通过图像配准来对齐三张图。

图像导入后，建议先做两个预处理：一是把图像转成 double 类型，避免 uint8 整型计算截断；二是做暗场校正，把镜头完全盖住拍一张或多张暗场图，从原图中减掉，消除传感器暗电流和热噪声。

核心计算阶段就是套公式。用矩阵运算比 for 循环快得多，MATLAB 里直接对整幅图像做乘加和开方就行。计算完 S0、S1、S2，再按公式算 DoLP 和 AoLP。最后用 imagesc 配合伪彩色显示。

下面给一个可直接运行的 MATLAB 脚本骨架，假设你已经把三幅图读进来，分别叫 I0、I60、I120：

% 核心计算模块：0°/60°/120° 方案 I0 = im2double(imread('I0.tif')); I60 = im2double(imread('I60.tif')); I120 = im2double(imread('I120.tif')); % 可选的暗场校正 dark = im2double(imread('dark.tif')); I0 = max(I0 - dark, 0); I60 = max(I60 - dark, 0); I120 = max(I120 - dark, 0); % 计算斯托克斯分量 S0 = (2/3) * (I0 + I60 + I120); S1 = (2/3) * (2*I0 - I60 - I120); S2 = (2/sqrt(3)) * (I60 - I120); % 计算偏振度和偏振角 eps_val = 1e-6; DoLP = sqrt(S1.^2 + S2.^2) ./ (S0 + eps_val); AoLP = 0.5 * atan2d(S2, S1); % 单位：度，范围：[-90, 90] % 显示结果 figure; subplot(2,2,1); imagesc(S0); axis image; colormap gray; title('S0 总光强'); subplot(2,2,2); imagesc(S1,[-0.5 0.5]); axis image; colormap gray; title('S1'); subplot(2,2,3); imagesc(S2,[-0.5 0.5]); axis image; colormap gray; title('S2'); subplot(2,2,4); imagesc(DoLP,[0 1]); axis image; colormap jet; colorbar; title('DoLP 偏振度');

这里给分母 S0 加了一个很小的 eps_val，是为了避免全黑像素处 S0 = 0 导致 DoLP 变成 NaN。后面章节我会专门讲这个坑。

2.3 核心计算模块代码（含公式对照）

如果你用的是 0°、45°、90° 方案，核心代码改成这样：

% 核心计算模块：0°/45°/90° 方案 I0 = im2double(imread('I0.tif')); I45 = im2double(imread('I45.tif')); I90 = im2double(imread('I90.tif')); S0 = I0 + I90; S1 = I0 - I90; S2 = 2*I45 - I0 - I90; DoLP = sqrt(S1.^2 + S2.^2) ./ (S0 + 1e-6); AoLP = 0.5 * atan2d(S2, S1);

代码很短，但真正难的是前面那三张图拍得够不够“一致”。我见过不少同学把时间都花在 MATLAB 调代码上，结果拍出来的三幅图因为光源闪烁、角度偏差大，怎么算都是噪声。与其反复调算法，不如回拍摄现场把偏振片角度校准做好。

3. 实操环节：从图片读取到结果可视化，每一步都别踩坑

3.1 图像读取与数据预处理：为什么必须用 im2double

MATLAB 里 imread 读出来的图像默认是 uint8，范围 0 到 255。如果直接用 uint8 类型做 S1 = I0 - I90，小于 0 的像素会被自动截断到 0，S2 的结果也会完全错误。所以第一步必须用 im2double 把图像转成 double 类型，同时把像素值归一化到 0 到 1 区间。

这里还有一个容易被忽略的点：相机输出的灰度值不一定是光强的线性响应。很多消费级相机会自动做 Gamma 校正和降噪，这会让灰度值与真实光强之间呈非线性关系。偏振计算公式里的 I(θ) 是线性光强，如果直接用非线性灰度值代入，偏振度会偏低。严谨做法是先对相机做辐射定标，把灰度值反变换成线性辐亮度；如果只是课程设计，在相机设置里关闭美颜、关闭自动对比度、关闭 Gamma 增强，并选择 RAW 或者 TIFF 格式，也能减少明显问题。

暗场扣除也很重要。把镜头完全盖住，拍一张和实验时相同曝光时间、相同 ISO 的暗场图，然后从每幅偏振图像中减掉。这一步能去除传感器暗电流和固定模式噪声。注意减完可能产生负值，记得用 max(..., 0) 把负值截断。

3.2 三幅图像配准和偏振片角度校准

如果偏振片是装在手动旋转架上，旋转过程中几乎不会改变图像几何；但如果是手持偏振片贴在镜头前，或者旋转时稍微压到了镜头，三幅图像之间可能出现亚像素级平移。偏振计算是逐像素进行的，哪怕 0.5 像素的偏移，在物体边缘处也会产生很大的偏振度伪影。

解决办法有两种：一是拍摄前就把相机和偏振片固定在一个刚性平台上，保证只有偏振片内部旋转，图像几何完全不变；二是事后用 MATLAB 的 imregister 或者相位相关方法做配准。配准顺序一般以 I0 为参考，把 I60 和 I120 对齐到 I0 上。对大多数实验场景，刚性配准就够了，不需要做复杂的非线性变形。

偏振片角度校准是另一个容易被忽视的环节。0° 方向是怎么定义的？如果只是凭刻度盘随手转到 0°，那么整套结果的绝对偏振角会整体偏移一个常数，但 DoLP 不受影响。如果你只关心相对偏振分布，问题不大；如果需要绝对角度，建议先拍一个已知偏振方向的目标做标定。我实际测过，角度偏差在 1° 以内时，DoLP 误差很小，但 AoLP 会有一个明显偏移；偏差超过 5° 后，DoLP 本身也会明显失真。

3.3 可视化技巧：伪彩色、直方图、显示范围

算完 DoLP 和 AoLP，用 imagesc 直接显示，会发现图像灰蒙蒙或者颜色发噪。原因是你没有手动设定显示范围。MATLAB 的 imagesc 默认会把数据最小值和最大值映射到整个色标，但 DoLP 的值域理论上是 0 到 1，绝大多数像素可能集中在 0.05 到 0.5，直接用默认范围会把小差异压缩成一片颜色。

建议这样设置：

figure; imagesc(DoLP, [0 1]); axis image; colorbar; colormap(jet); title('DoLP 偏振度');

想看低偏振目标时，可以把上限调低，比如 [0 0.3]，这样会更清楚。AoLP 用 hsv 色标比较合适，因为偏振角是周期性的，0° 和 180° 在物理上是同一个方向，hsv 色标首尾相接，看起来更自然。但注意 AoLP 默认范围是 [-90, 90]，如果你想显示成 [0, 180]，用 mod(AoLP, 180) 转换一下，避免颜色跳变一半。

S1 和 S2 的取值范围是对称的，中心是 0，通常显示为 [-0.5, 0.5] 或 [-max(abs(S1(:))), max(abs(S1(:)))]。显示时用双色色标，比如 blue-white-red，能直观反映正负。

3.4 保存结果并叠加信息

保存 DoLP 图像时，imwrite 对 double 类型默认要求数据范围在 0 到 1，如果直接保存带 NaN 或大于 1 的值，输出可能会变成全黑或全白。可以先把数据裁剪到 [0, 1]，再转成 uint8 保存：

DoLP_save = uint8(round(DoLP * 255)); imwrite(DoLP_save, 'DoLP.png');

如果想把多幅结果拼成一张图，用 imwrite 配合 RGB 通道。例如把 S0 放到 R 通道，DoLP 放到 G 通道，AoLP 归一化后放到 B 通道，可以生成一张信息叠加的伪彩色图。这样后续检查数据时，不用反复打开多个文件。

4. 常见问题与排查思路

4.1 问题速查表

4.2 偏振度图像出现大片 NaN 或 Inf 怎么办

这是最常见的报错现象。S0 在总光强为 0 的区域等于 0，比如纯黑背景、镜头遮挡区域。此时直接算 sqrt(S1.^2 + S2.^2) ./ S0 会得到 NaN，MATLAB 显示时这些像素会变成空白或异常值。

推荐做法是提前构造一个掩码：

mask = S0 > 1e-3; % 只统计有足够光强的像素 DoLP = zeros(size(S0)); DoLP(mask) = sqrt(S1(mask).^2 + S2(mask).^2) ./ S0(mask); DoLP(~mask) = 0; % 或者设为 NaN，但不建议直接显示

如果只是简单加 1e-6，在极暗区域会出现 DoLP 虚高，因为分子本身也是噪声主导。所以更严谨的做法是用阈值掩码：只有 S0 超过某个阈值才计算偏振度，否则直接置 0。

4.3 角度误差对结果的影响有多大

偏振片的角度误差是很多人忽略的系统误差。理论公式假设三个角度精确等于 0°、60°、120°，如果实际变成 0°、58°、123°，S1、S2 的计算就会引入偏差。我做过一个简单仿真：5° 的角度误差会使 DoLP 的相对误差达到 10%，AoLP 的误差可能达到几度到十几度，具体取决于场景的偏振态分布。

所以在拍摄前，一定要确认偏振片的刻度盘和实际透光方向一致。偏振片通常标有透光轴方向，不要把刻度线的方向当成透光方向。更科学的做法是找到偏振片透光轴方向后，再对准 0° 刻度。如果条件允许，用步进电机控制的偏振轮拍摄，角度重复精度可以做到 0.1° 以内，结果稳定性会好很多。

4.4 算法提速与批处理优化建议

对单张 1080p 图像，上面的矩阵运算毫秒级完成，完全不需要优化。但如果你要处理一批图像序列，或者做实时视频偏振成像，就有必要优化。

首先把核心计算封装成函数：

function [S0, S1, S2, DoLP, AoLP] = calc_stokes(I0, I60, I120) % 输入：三幅同尺寸灰度图像(double 类型) % 输出：斯托克斯分量、偏振度、偏振角 S0 = (2/3) * (I0 + I60 + I120); S1 = (2/3) * (2*I0 - I60 - I120); S2 = (2/sqrt(3)) * (I60 - I120); S = sqrt(S1.^2 + S2.^2); DoLP = S ./ (S0 + 1e-6); AoLP = 0.5 * atan2d(S2, S1); end

然后用 parfor 对文件序列做批处理。注意 parfor 里不要频繁调用 imshow 和 imwrite，先把结果存到内存或临时变量，最后统一写盘。图像数据如果很大，可以用 single 类型计算，内存占用减半，速度略快，精度对大部分应用足够。

5. 结果怎么看：从图像到结论的实战经验

5.1 偏振度图像亮暗代表什么

DoLP 图像中，亮的区域表示偏振度高，暗的区域表示偏振度低。一般来说，镜面反射、光滑塑料表面、水面反射、天空散射光某些角度，偏振度较高；粗糙表面、漫反射布料、纸张等，偏振度很低。

在目标检测场景里，偏振度图像可以帮忙把高亮反光目标和背景分离开。比如室外拍摄时，普通灰度图里阳光下的水面和大楼玻璃都很亮，但在 DoLP 图像里它们的偏振特性和其他物体明显不同。你可以用 DoLP 做阈值分割，提取强反射区域。

不过要记住，DoLP 是“相对”量，它会受到光源角度、观察角度影响。同一个物体，换个角度拍，DoLP 可能差异很大。所以不要直接拿 DoLP 绝对值当物体的固有属性，更适合的做法是作为特征送入分类器，或者结合 S1、S2 一起分析。

5.2 偏振角图像的颜色跳变怎么解释

AoLP 图像用 hsv 色标显示时，有时候会看到颜色在某个区域突然从红色跳到蓝色，看起来像缺陷。这通常有两种原因：一种是真实的偏振方向接近 90° 边界，微小的物理变化导致显示颜色循环跳变；另一种是低偏振区域的噪声主导，AoLP 没有物理意义，颜色随机跳动。

第二种情况很常见。DoLP 接近 0 的像素，S1、S2 都接近 0，atan2(0, 0) 结果是 0 或 NaN，即使有微小噪声角度也会乱跳。解决方法是显示 AoLP 时，把 DoLP 低于阈值的像素置为灰色或黑色：

AoLP_display = AoLP; AoLP_display(DoLP < 0.05) = -90; % 或者设为 NaN imagesc(AoLP_display, [-90 90]); axis image; colorbar; colormap(hsv);

这样就能避免颜色噪声干扰，只展示可信区域的偏振方向。真实场景中 AoLP 图像能反映物体表面朝向和光源相对位置，比如光滑曲面反射区域的偏振角会沿着曲面轮廓变化，这在三维重建和材质分析里很有用。

5.3 斯托克斯分量图适合哪些应用

S0 就是总光强图，跟普通灰度图几乎一样，但它与普通相机图存在一个细微差别：已经去除了偏振片带来的角度相关调制。S1 和 S2 则分别代表水平和 45° 方向的偏振差异，它们对边缘、纹理、材质边界很敏感。

在去雾应用里，大气散射光的偏振度通常较高，可以用 S1、S2 估计大气光偏振方向，然后从 S0 中减去散射成分，提高清晰度。在工业检测里，S1、S2 图像可以突出划痕、应力双折射区域，因为这些区域会局部改变偏振态。在遥感领域，偏振度图像常用于地表分类，比如区分土壤、植被和水体。

斯托克斯分量图像的一个好处是：它们都是线性组合结果，可以直接做后续代数运算。比如目标增强时，可以用 S1、S2 构造偏振差分特征，而不需要重新拍图。理解了 S0、S1、S2 的物理意义，你可以根据自己的需求设计更多特征组合，而不仅仅局限于 DoLP 和 AoLP。

我做这类实验最大的体会是：第一步把拍摄环节的光源稳定性和角度校准做好，比后面调 MATLAB 代码重要得多。偏振成像的输入数据如果不可靠，再漂亮的伪彩色图也只是噪声的另一种形式。如果你只是完成课程大作业，按上面的流程跑通已经足够；如果想深入研究，建议用多个角度（比如 0°、30°、60°、90°）拍四幅以上图像，再用最小二乘拟合斯托克斯参数，稳噪效果会更好。先动手跑通三幅图的方案，再慢慢扩展，你会理解得比我当初快很多。

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