偏振成像与MATLAB实现：从三角度图像到DoP/AoP参数提取

简介：本资源是一套基于MATLAB实现偏振成像分析的完整实验代码与数据集，面向光学工程、机器视觉及图像处理方向的本科生、研究生和科研人员，解决从原始偏振图像获取斯托克斯参量、偏振度（DOLP）与偏振角（AOP）图像的核心计算问题。压缩包共9个文件（540KB），含8幅BMP格式原始偏振图像（对应偏振片在0°、45°、90°等三组角度采集）及1个主运行脚本pianzhenxinxiyuanshi.m，脚本封装了斯托克斯矢量求解、DOLP与AOP像素级计算及可视化流程，可直接运行复现实验结果。已有1356人学习下载，资源结构简洁明确：原始图像按物理意义命名（如hh1.bmp、hp1.bmp、hv1.bmp），输出图像（I.bmp、U.bmp、DOLP.bmp、AOP.bmp）与中间结果一一对应，便于理解偏振参量的物理含义与算法逻辑，是开展偏振光信息提取教学实验或算法验证的理想入门材料。 偏振成像这事儿，听起来像是实验室里才会碰到的课题，但只要你手里有一台普通相机、一片偏振片，再会一点MATLAB，就能从三张看似“只是变暗了一点”的照片里，把物体表面人眼看不见的物理信息挖出来。偏振度（DoP）告诉你反射光中有多少是规则偏振的，偏振角（AoP）则能反推反射面的朝向和材质特性。这篇博文就围绕“三角度偏振图像到偏振参数图”这条主线，把从物理原理、实验采集到MATLAB实现的完整链路拆开揉碎讲清楚，适合正在做图像处理大作业、搞机器视觉项目，或者刚接触偏振成像想快速上手的同学参考。

1. 为什么是三张图？偏振成像的数学基底与测量策略

1.1 从光强的“叠加态”说起：斯托克斯参数到底在描述什么

普通相机记录的灰度值，只是光强这个标量。但光作为横波，它的电场矢量在垂直于传播方向的平面内是可以有特定振动方向的，这就是偏振态。一束光可能是完全非偏振的（自然光），可能是完全线偏振的（比如经过反射或散射后的光），更多时候是两者的混合。

斯托克斯参数就是用来定量描述这种混合状态的四个数，通常写作S0、S1、S2、S3。其中S0是总光强，S1描述水平/垂直线偏振的差异，S2描述45°/135°线偏振的差异，S3涉及圆偏振分量。对于绝大多数被动偏振成像场景——比如天空光、水面反射、塑料或金属表面的反射——圆偏振分量微乎其微，所以我们只需要关注S0、S1、S2这三个参数就够了。

那么怎么从测量中“解出”S0、S1、S2呢？核心思路是：让光通过一个可旋转的线偏振片，改变偏振片的透光轴方向，记录不同角度下的光强I(θ)。理想情况下，透过偏振片后的光强与入射光的斯托克斯参数满足一个余弦关系：I(θ) = 0.5 * (S0 + S1 * cos(2θ) + S2 * sin(2θ))。

这里θ是偏振片透光轴相对参考方向（通常是水平方向）的旋转角。也就是说，只要你在多个角度下测量光强，就能反推出未知数S0、S1、S2。

1.2 三次测量与四次测量：角度选择的博弈

从数学上看，I(θ)公式里有三个未知数，理论上最少需要三个不同角度下的独立测量。于是就有了两种最常见的方案：

• 四强度法：取0°、45°、90°、135°四个角度。此时S0 = I(0) + I(90)，S1 = I(0) - I(90)，S2 = I(45) - I(135)。计算最直观，每个参数都有清晰的物理对应，是教材里最常见的方案。
• 三强度法：取0°、60°、120°三个角度。这也就是你这个项目标题里提到的方案。它利用了余弦函数在特定角度下的对称性，只需要旋转三次，而且旋转步长均匀（60°），机械实现上更简单，更适合连续转动的电动旋转台。

三角度方案的公式推导也很优雅。把θ = 0°、60°、120°分别代入I(θ)的表达式，得到一个三元一次方程组：

I(0) = 0.5 * (S0 + S1) I(60) = 0.5 * (S0 - 0.5S1 + 0.866S2) I(120) = 0.5 * (S0 - 0.5S1 - 0.866S2)

（其中S1和S2前面的系数来自cos和sin在60°、120°处的值）

解这个方程组，就能得到三个斯托克斯参数的显式解。在MATLAB里你可以直接写矩阵求解，也可以手算化简成更简洁的形式。我建议在代码里保留线性方程组的写法，这样以后想换角度组合（比如改成0°、45°、90°），只需要改矩阵系数就行。

1.3 角度设计里的隐藏代价：误差放大与噪声均衡

你可能想问：为什么市面上有些设备用0°、45°、90°、135°，而这里偏偏选0°、60°、120°？除了减少一次拍摄和旋转的机械便利，还有一个原因是误差在各方向上的均衡性。

在三强度方案里，三个角度互相错开60°，它们在偏振态空间的分布是均匀的。当你把这三个强度值当作线性方程组的输入时，系数矩阵的条件数相对较小，数学上意味着噪声不会在某些方向上被过度放大。但代价是：S1和S2的计算不再是简单的加减法，而是会引入系数相乘（比如S2表达式中会有除以√3的项），这会让结果对像素噪声更敏感一点点。

所以，如果你是在静态场景下实验，相机和物体都不动，我其实更推荐用四强度法，因为计算简单、不容易出错；如果场景对拍摄次数敏感（比如有微光漂移），或者你用的旋转台只能按固定步长转动，那三强度法就是更好的选择。这个项目用三角度，说明应用场景大概率是后者——记住这个前提非常重要，因为后文的噪声处理方法完全是围绕它展开的。

2. 图像采集环节的纪律：偏振片角度、相机设置与配准检查

2.1 偏振片安装与角度基准的确定

偏振片不是随便夹在镜头前就完事的。首先要确定“0°基准”。把偏振片装到镜头前时，最好用光轴方向（比如桌面的镜面反射光）来标定：当偏振片透光轴与反射光的偏振方向垂直时，反射光会被大幅消光，画面会变得很暗。以此为基准，把这一位置记为偏振片的机械零点，或者直接在镜头上做个标记，接下来旋转60°、120°就有参照了。

如果用的是手动旋转环，我建议拍三张之后，再补拍一张0°重复，用来检查相机和光源的稳定性。如果场景中某些高光区域在两次0°拍摄中的灰度值差异超过1%—2%，说明光源有闪烁或者相机没固定好，这时候得先解决稳定性问题，否则后面算出来的DoP和AoP都会带条纹状噪点。

2.2 相机参数与图像格式：RAW优先，避免非线性变换

偏振度计算本质上是强度值的比值运算，所以相机采集到的灰度必须与入射光强成线性关系。但相机内部默认的γ校正、自动增益、白平衡，都会把线性关系破坏掉。

我做这个实验时踩过一次坑：用JPG直出，算出来的偏振度图边缘总是怪怪的，后来才反应过来是JPG做了非线性色调映射。解决方法是：有RAW就拍RAW，在MATLAB里用rawread读取，然后自己做一个简单的线性去马赛克，或者直接调用demosaic函数，它会保留线性强度信息。找不到RAW功能的话，至少在相机菜单里把“照片风格”设为“可靠设置”或“中性”，把对比度、饱和度全部拉到最低，同时关闭自动亮度优化。

镜头光圈也不宜开太大。有些镜头的镜片本身会有应力双折射，在边缘区域引入伪偏振。我用50mm f/1.8镜头在f/1.8光圈下测过一次，边缘区域的DoP会比中心高出2%—3%，收两档光圈到f/4以后，整个画面就干净了。光圈收小带来的曝光不足，用延长快门或者提高ISO来补，但ISO不要超过800，否则噪声会直接影响S1、S2的精度。

2.3 像素级对齐：旋转带来的微小位移是重灾区

手动旋转偏振片时，很难保证旋转轴与相机光轴完全重合。哪怕偏移只有几十个微米，在2000万像素的传感器上也会造成几个像素的位移。而偏振度计算是对每个像素独立进行的，如果同一物理点在三次拍摄中分别落在不同的像素上，算出来的DoP和AoP就会出现边缘重影和伪纹理。

稳妥的做法是在镜头后端加一个转接环，让偏振片在镜头前端的滤镜螺纹处旋转，这样轴心基本稳定。如果你用的偏振片是插入式滤镜，那就要在软件层面做配准。常用的方法是在场景里放一个特征点丰富的标靶，用MATLAB的detectSURFFeatures和estimateGeometricTransform2D做配准，以第一张0°图作为参考，把60°和120°图通过仿射变换对齐。注意，配准时要用外插还是内插？我建议用imwarp时把'OutputView'设为'same'，并且用'cubic'插值，双线性插值会在边缘引入模糊，对偏振度影响虽然不是特别大，但做精细研究时会看出来。

3. MATLAB计算流程：从像素值到偏振参数图的完整实现

3.1 数据准备与预处理

下面假设你已经把三张图像读入MATLAB，分别命名为I0、I60、I120。这里有个细节：如果图像是RGB彩色图，你有两个选择：一是转成灰度图再算；二是分别对R、G、B三个通道独立计算偏振参数，最后再做通道融合。前者的好处是速度快，后者能保留光谱信息——偏振度在彩色图上其实比灰度图更丰富，因为不同波长的光在物体表面的反射偏振特性并不相同。对单次实验而言，直接用灰度图就能验证整个算法流程。

第一步，把图像统一转成double类型，并把值域放缩到0到1之间。这一步极其重要，因为uint8的数据在做减法和除法时会溢出，导致S1、S2算出来全是0或255。

I0 = im2double(imread('polar_0.jpg')); I60 = im2double(imread('polar_60.jpg')); I120 = im2double(imread('polar_120.jpg'));

如果图像有坏点或热噪声，可以先做一个3x3的中值滤波，但滤波窗口不要超过5x5，否则会抹掉细小的偏振结构。我通常的做法是只对明显异常的孤立点做修复，用medfilt2配合imflatfield做背景平坦化，看情况使用。

3.2 构建线性方程组并解算斯托克斯参数

前面提到，测量强度I(θ)与斯托克斯参数的关系是：I(θ) = 0.5 * (S0 + S1cos(2θ) + S2sin(2θ))。三个角度对应三个方程，写成矩阵就是：

theta = [0, 60, 120] * pi/180; A = 0.5 * [1, cos(2*theta(1)), sin(2*theta(1)); 1, cos(2*theta(2)), sin(2*theta(2)); 1, cos(2*theta(3)), sin(2*theta(3))];

然后用矩阵除法直接求出S0、S1、S2：

% 对每个像素独立解算，使用reshape优化为矩阵运算 sz = size(I0); M = [I0(:), I60(:), I120(:)]'; % 3 x N Stokes = A \ M; % 3 x N S0 = reshape(Stokes(1,:), sz); S1 = reshape(Stokes(2,:), sz); S2 = reshape(Stokes(3,:), sz);

这里用A \ M而不是手动写出解析解公式，好处是当时如果你临时改用了0°、45°、90°三个角度，只需要重新定义theta向量和A矩阵，后面代码一行都不用改。

从三次测量推导出S0、S1、S2后，可以画一张三通道伪彩色图来看效果。把S0作为亮度通道，把S1、S2分别映射到色度通道，就能直观看到物体表面偏振态的连续变化。这是调试阶段最直观的手段。

3.3 偏振度（DoP）与偏振角（AoP）的计算

偏振度定义为线偏振分量占总光强的比例：DoP = sqrt(S1^2 + S2^2) / S0。用MATLAB写，一行向量化就能完成：

DoP = sqrt(S1.^2 + S2.^2) ./ (S0 + eps);

这里加eps是为了防止S0为0的像素处除零。理论上DoP的值域是[0, 1]：0表示完全非偏振，1表示完全线偏振。实际数据里，由于噪声和S0估计误差，DoP可能会出现略微超出[0, 1]的情况，尤其是暗区。处理方法是直接截断到[0, 1]：

DoP = min(max(DoP, 0), 1);

偏振角表示偏振椭圆长轴的方位角，公式是AoP = 0.5 * atan2(S2, S1)。注意这里必须用atan2而不是atan，因为atan的值域只有[-π/2, π/2]，无法区分S1和S2同号或异号的情况，会导致偏振角方向在45°处发生跳变。用atan2直接得到的是弧度值，值域是[-π/2, π/2]——不对，atan2(S2, S1)的值域是[-π, π]，除以2后是[-π/2, π/2]，但物理上偏振角应该在[0°, 180°)区间。所以要做一个角度包装：

AoP = 0.5 * atan2(S2, S1); % 弧度 AoP = AoP * 180 / pi; % 转角度 AoP = mod(AoP, 180); % 映射到[0, 180)

这段代码里的mod操作是关键，我见过好几个人因为漏了这一步，画出来的偏振角图在0°和180°交界处出现很突兀的锯齿。实际上偏振角物理上就是180°周期的，0°和180°是同一个方向，不包装的话可视化会非常让人困惑。

3.4 把三个斯托克斯分量图像化的技巧

有些同学会直接把S0、S1、S2当作灰度图显示，结果发现S1、S2图里有一半像素是负值，灰度图没法直接显示，于是就得做线性映射。

我自己的习惯是：S0直接用原图灰度显示；S1、S2则用imagesc加colormap(jet)显示，并让色标中心对应0值，这样就能清晰看到正负区域的空间分布。比如S1正值表示该区域水平偏振分量大于垂直偏振分量，负值则反之。这个在分析材料的反射特性时很有用。

如果你想把S1、S2映射到0—255的灰度图输出，推荐用：

imagesc(S1); colorbar; caxis([-max(abs(S1(:))) max(abs(S1(:)))]);

归一化的时候不要用mat2gray，因为mat2gray会把最小值映射到0，这样就把正负信息破坏了。记住，S1、S2是有符号物理量，可视化时必须保留正负。

4. 我调试过程中踩过的几个关键坑：值域、角度绕圈与噪声放大

4.1 uint8陷阱：亮度减出“雪花”纹

第一次跑通这个流程时，我用的是imread读进来的uint8类型直接做减法，结果S1图像上一片雪花状的噪点。原因是uint8减uint8，在MATLAB里默认是饱和运算：如果被减数小于减数，结果直接变成0而不是负数。这等于把S1的所有负半轴信息全部截断了，偏振角当然就乱了。

解决办法就是在一切计算之前先im2double。这个操作花不了多少时间，但能避免至少半小时的抓狂排查。我甚至建议不要用double(imread(...))/255这种写法，因为如果原始图像是16位TIFF，除以255就会把值域弄得不对。im2double会自动判断图像类型并正确归一化，省心很多。

4.2 atan2还是atan：偏振角的四象限歧义

我在前面已经强调过了，这里再说一个更具体的例子。假设某个像素的S1 = -1，S2 = 1，用atan(S2/S1)计算，得到的是-45°，然后你把它映射到[0°—180°)，得到135°。但用atan2(S2, S1)直接得到的是135°。这说明atan会丢符号信息，导致45°的偏差。如果你的数据里既有正S1又有负S1的区域，用atan算出来的AoP图会像地图一样出现对角方向的错位。所以别偷懒，老老实实用atan2。

4.3 S0太小导致的偏振度爆炸

在阴影或黑色物体的区域，S0趋近于0，但S1和S2并不一定也趋近于0——它们可能是小的正值或负值。这时DoP = sqrt(S1^2+S2^2)/S0，本质上是用两个小噪声除以另一个小噪声，结果会有很大的随机波动，甚至超过10。虽然我已经在代码里进行了min(max())截断，但截断只是掩盖问题，会生成一片像素值全是1的区域，看起来像“过曝”。

更稳妥的办法是做一个阈值掩膜：当S0小于某个阈值（比如全图S0最大值的5%）时，把DoP和AoP设置为无效值（NaN），在可视化时用黑色或透明色标出。这样做既保留了有效区域的偏振信息，也避免了视觉上的误导。

threshold = 0.05 * max(S0(:)); mask = S0 > threshold; DoP(~mask) = NaN; AoP(~mask) = NaN;

4.4 边缘配准残影：绕不开的对齐问题

如果你拍的是静态场景，用三脚架固定相机，只旋转偏振片，但三张图之间仍然可能有一两个像素的位移。这种残影对DoP的直接影响是边缘处出现亮暗交替的“镶边”。特别是景物边缘，反射率变化剧烈，哪怕一个像素的错位，都会让S1、S2在边缘产生很大的正负跳变。

我在实际项目中用过两种解决方式。第一种是用图像配准，上文已经提过；第二种更粗暴但很有效——对最终算出的DoP图做一个2到3像素的形态学腐蚀，把边缘那一圈异常值剥掉。腐蚀会损失一点空间分辨率，但能让结果图非常干净。学术论文里这个叫“边界效应去除”，没有人会说你做得不对。

4.5 光源稳定性与背景扣除

这个坑不太容易注意到。如果你在室内用钨丝灯或LED灯照明，这些光源本身是部分偏振的。尤其是LED，其荧光粉散射可能导致光源本身带有几个百分点的偏振度。光源偏振会叠加在物体反射的偏振信号上，导致你测的DoP并不是物体的真实属性。

如果你做的是高精度定量测量，建议在拍摄前把偏振片放在光源前，旋转偏振片看光强变化：如果变化幅度超过1%，说明光源本身偏振明显，需要加一个退偏振片。对于大作业或工程验证级别的精度，这个因素通常可以忽略，但心里要有数。

5. 结果可视化的进阶玩法：伪彩色映射与数据验证

5.1 偏振度图的灰度与伪彩色输出

DoP图本身是一个灰度图像，直接用imshow显示即可。但灰度图对细微偏振差异的辨识度不高，所以我一般会用colormap('jet')或parula做伪彩色映射。特别是当DoP值集中在0到0.3这个低动态范围时，用线性灰度显示会显得非常平淡，伪彩色能瞬间拉开层次。

要注意的是，显示用伪彩色，但保存数据时一定要保存原始的DoP矩阵，最好存成MAT文件或GeoTIFF格式（如果你有坐标信息）。不要直接把伪彩色图存成JPG当作结果数据，因为JPG压缩会损失量化精度，而且读取的时候你已经拿不回原始的浮点DoP值了。

5.2 偏振角图的HSV映射：方向信息如何优雅展示

偏振角AoP是一个周期量，用普通线性色标会出现一条从白色突变到黑色的分割线，视觉上非常不自然。我推荐用HSV色彩空间来表示AoP：把色调H映射为AoP的角度值，把饱和度S设置为1，把亮度V设置为S0归一化的结果。这样在任何区域，亮度高的地方表示信号强，色调表示偏振方向，可以在一张图里同时表达两个维度的信息。

MATLAB里可以这么写：

H = AoP / 180; % 归一化到[0,1] S = ones(size(AoP)); % 最大饱和度 V = mat2gray(S0); % 亮度用光强 hsvImage = cat(3, H, S, V); rgbImage = hsv2rgb(hsvImage); imshow(rgbImage);

生成的效果很惊艳：不同朝向的表面呈现出不同的颜色，高光区域颜色鲜明，阴影区域暗淡。这是我在项目汇报PPT里最喜欢放的一张图。

5.3 用已知目标验证算得对不对

算法写完别急着上真实验数据，先用一个已知偏振特性的目标做验证。最简单的方法：找一片偏振片，放在一个均匀光源前面（比如手机屏幕显示纯白），然后旋转这个偏振片到已知角度，用你的系统去测。这时候得到的AoP应该就是那个已知角度，误差应该在±3°以内。测几个角度，如果都准，说明你的系统标定和算法流程是可信的。

没有偏振片的话，也可以用天空作为验证目标。晴朗日子里的天空光，在很多角度下是高度线偏振的，特别是在与太阳成90°的方向。用你的系统拍三张天空图，算出来的DoP应该在合理范围（通常0.3—0.7），AoP的分布也应该是平滑过渡的。我刚开始做偏振成像时，就是用天空做验证，发现算法和相位对齐都没问题后才放心跑实际数据的。

6. 算法封装与效率优化：把脚本变成可复用函数

6.1 函数化设计：输入三张图，输出五个参数

把核心计算封装成一个函数，输入是I0、I60、I120三张尺寸相同的图像，输出是S0、S1、S2、DoP、AoP五个矩阵。函数内部按前文的流程走一遍。这样以后再拿到不同场景下的偏振数据，一行代码就能出结果，不用再复制粘贴脚本。

function [S0, S1, S2, DoP, AoP] = calcStokesFrom3Angles(I0, I60, I120, thetaDeg) % thetaDeg 是三个测量角度的向量，默认 [0, 60, 120] if nargin < 4 thetaDeg = [0, 60, 120]; end % ... 实现见上文 end

6.2 内存与速度的平衡

三张2000万像素的图像，如果用double类型存，每张约160MB，三张就是约480MB，加上计算过程中的临时变量，有些低配电脑会卡顿。我建议在读入图像后立即转成double，但计算过程中尽量用向量化操作。如果内存还是吃紧，可以分块处理：把图像切成若干个512x512的块，逐块计算，最后拼回去。我通常用blockproc实现这个，它允许你用自定义函数处理每个块，非常方便。

fun = @(block) calcStokesBlock(block.data, A_matrix); result = blockproc(cat(3, I0, I60, I120), [512 512], fun, 'BorderSize', [2 2]);

注意，分块处理会引入边界效应问题——每个块边缘的像素在滤波或插值时，需要相邻块的数据。blockproc的'BorderSize'参数就是干这个的，设成2的时候，每个块会包含2像素宽的重叠区域，拼接结果不会出现块状痕迹。

6.3 与后续应用的衔接：偏振特征如何接入视觉任务

算出了DoP和AoP，这还不是终点。偏振成像的价值在于，它能提供比普通灰度图像更丰富的信息。两个灰度值相同的物体，可能在DoP图像上有明显差异——比如水面和路面、塑料和金属、疤痕组织和正常组织。你可以把DoP和AoP当作额外的通道，和RGB图像拼在一起，形成一个五通道或六通道的数据立方体，再送入深度学习网络做分割或识别。

我在一个表面缺陷检测项目里就是这么干的：在灰度图上几乎无法分辨的细微划痕，在AoP图像上因为改变了局部反射的偏振方向，变得非常显眼。当时我用了一个非常朴素的两层卷积网络，只输入DoP和AoP两个通道，检测率就比用灰度图提升了将近20个百分点。偏振信息对纹理弱、几何特征不明显的表面缺陷，几乎是一种天然的增强。

7. 精度评估与误差来源的定量分析

7.1 偏振片角度误差带来的后果

机械旋转偏振片时，角度读数可能误差±0.5°。这个误差对DoP的影响有多大？我做过一个简单仿真：在用理想公式生成I(0)、I(60)、I(120)的基础上，给每个角度加入随机高斯噪声（标准差0.5°代表旋转台回程误差），然后比较解算出的DoP与真实DoP的误差。结果是在DoP=0.5附近，角度误差导致的DoP相对误差约为2%—3%。换句话说，如果你不是做超高精度的定量测量，0.5°的旋转精度完全够用，没必要花大价钱买精密旋转台。

7.2 噪声对偏振角的影响及信噪比门槛

偏振角对噪声更敏感。在S0比较低的区域，S1和S2的信噪比下降，atan2的结果会发生随机抖动。做一个简单的估计：假设S0=0.2，S1和S2都是0.05，噪声标准差0.01，那么AoP的不确定度大约在5°左右。如果你的应用需要偏振角的精度在1°以内，那你必须保证暗区的光强不能太低，否则就得延长曝光时间或使用更高的ISO。

从实践角度来说，我一般会在拍摄前看直方图：确保绝大多数像素的亮度值都落在满量程的20%—80%区间。如果高光区域过曝，S0会被截断，偏振度偏小；暗部太暗，噪声主导，偏振角乱跳。中间调是最均衡的。

7.3 偏振片消光比的影响

市面上便宜的线偏振片消光比可能只有100:1到300:1，好的比如消光比在10000:1以上。消光比不够高意味着：即使透光轴与偏振方向垂直，仍然有一小部分光漏过来。这个漏光对偏振数据的影响是系统性偏差：实测的DoP会偏低，尤其是高偏振区域，因为泄漏光是非偏振的，会把“纯偏振信号”稀释掉。

在你用三角度方案做实验时，如果发现测一个已知的完全偏振光源，DoP只有0.85而不是接近1，那多半就是偏振片消光比不足。这时候可以在结果上做一个修正：DoP_true = (DoP_measured - 1/ER) / (1 - 1/ER)，其中ER是消光比。虽然这个修正模型是近似处理，但能把系统偏差拉回来不少。

8. 从三角度到四角度：算法换一个输入怎么改

8.1 四角度方案的矩阵表达

如果你的实验条件允许拍摄四张图，也就是0°、45°、90°、135°，那计算会更简单：S0 = I0 + I90，S1 = I0 - I90，S2 = I45 - I135。不需要解方程组，直接用加减法就行。这一方案特别适合那些需要实时处理的场景，因为排除除法运算后，计算量大幅下降。

需要改动的只有矩阵A和读入图像的数量。你可以在函数里加一个参数控制“测量模式”，如果是三角度，就用A矩阵解方程；如果是四角度，就走加减法分支。

8.2 圆偏振分量S3与旋转偏振片的局限

偏振片只能测量强度在某个线偏振方向上的投影，也就是说它只能测到S0、S1、S2三个分量，对S3（圆偏振分量）完全不敏感。想测S3，需要额外加一个四分之一波片，改变相位延迟后才能分离出S3。如果你真的需要完整的斯托克斯矢量（比如做水下主动偏振成像），那就要考虑在光路里加入波片，MATLAB算法层面对应的就是多几个输入图，多解几个未知数。

对大多数课程项目和工业应用，线偏振部分的DoP和AoP已经提供了足够的信息，先不用急着上全套斯托克斯测量。

9. 个人经验与后续扩展

讲真，偏振成像这个领域，入门最大的障碍不是数学，也不是MATLAB代码，而是采集数据时一丝不苟的实验习惯。三个角度的图像只要有一张出现运动模糊、焦距漂移或者光源闪烁，后面整个计算流程做得再漂亮也白搭。我在做这个项目时，每次采集前都会先拍三张测试图，跑一遍全流程，看到DoP和AoP分布符合预期后，才开始正式采集多组数据。这个“预扫”习惯帮我省了大量返工时间。

后续可以做的扩展很多。比如把采集过程做成自动化的：用Arduino控制步进电机旋转偏振片，同时给MATLAB发送触发信号，实现多角度的自动采集和处理。我曾经用这种方式在一个小时内连续采集了20组不同场景的数据，大幅提高了效率。

另外一个非常值得尝试的方向是动态偏振成像：把偏振相机装在无人机或者机械臂上，对同一场景从不同视角拍摄，然后结合三维重建算法，同时恢复物体的几何形状和表面材质偏振特性。这在遥感、考古、文物保护、刑侦取证里都有应用潜力。不过这个方向的技术栈就比单纯算DoP要复杂得多了，属于另一个话题了。

最后，代码里最重要的一个习惯：每算出一个中间量，就立刻figure; imshow()看一眼数据是否合理。偏振成像里，你对着公式写代码容易，但数据是肉眼可见的，一旦某个环节错了，图像上立刻会暴露出来。调试效率最高的路，永远是“算一步、看一步、验证一步”。

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