MATLAB符号数学解非线性方程组 -- 眼球反射模型

1.利用符号函数求解非线性方程组

1 定义符号变量 syms

2 方程求解

　　solve：求解出来的是个解析解，意味着其实不用matlab，自己一直算也能算出一个通式；

　　vpasolve：数值解，可理解为是逼近值，自己算是不能算出通式的。

　　而实际我碰到的问题求解非线性方程组没有通式，所以考虑采用vpasolve

2.实例分析 球面反射点glint三维坐标求解

　　 红外LED灯照射眼球，然后反射，摄像头捕捉到反光点。由于眼球角膜是球面反射，没有通式，所以需要用数值方法求解。

　　假设反射点为G(xg,yg,zg)，包含三个未知数，需要找到三个独立的约束方程求解，包括：1，反射点glint在角膜球面上；2，角膜球心O，摄像机C，LED灯位置L，以及反射点G，四点共面；3，发射定律，入射角等于反射角。使用vpasolve求解如下：

solG = vpasolve([norm(OG) == R1,OG*subs(n,soln)',a == b],[xg,yg,zg],[nan,nan;nan,nan;0,inf]);

其中a,b为入射角以及发射角余弦公式，即

a = OG * GL'/(norm(OG)*norm(GL));
b = OG * GC'/(norm(OG)*norm(GC));

　　求解出来的solG，仍然是符号变量：

　　solG =

　　　　xg: [1x1 sym]

　　　　yg: [1x1 sym]

　　　　zg: [1x1 sym]

　　需要带入原始反射点G中，得到：

G2 = subs(G,solG)
G2 =
[ -6.6928147380739451481396476913942, 0.80366754303004353771951008713299, 4.3081723923376420662103160847784]

　　属性仍然是符号变量sym，将其转换为double，得到：

G2 = double(subs(G,solG))
G2 =
-6.6928    0.8037    4.3082

　　MATLAB源代码如下，增加了一个球体sclera巩膜：

%% 1 建立眼球模型，求相机和LED灯任意位置的glint
% 1.1 原始眼球模型
% 人眼正前方为Z轴，朝上为X轴，右手定则得Y轴为人眼右边
clc,clear
O1 = [0,0,0];           %cornea
R1 = 8;
[alpha1,beta1] = meshgrid(0:0.05:pi,0:0.05:pi);
O1O2 = 7;               %两球间隔设置为7mm
O2 = [0,0,-O1O2];       %sclera巩膜
R2 = 12;
[alpha2,beta2] = meshgrid(0:0.05:pi,0:0.05:2*pi);
P = [0,0,R1];           %pupil中心
% 设定camera和led的位置，头戴式，桌面式
% L = [-250,0,500];
% C = [-200,0,500];
L = [-40,10,30];
C = [-55,0,30];
% 1.2 求glint空间坐标
% 不存在解析解，只能用寻找数值解
syms xn yn zn xg yg zg
n = [xn,yn,zn];
G = [xg,yg,zg];
GL = L - G;
GC = C - G;
% 求平面法向量，约束方程--四点共面
OL = L - O1;
OC = C - O1;
soln = vpasolve([n*OL',n*OC',norm(n) == 1],[xn,yn,zn]);
% 求glint的位置，3个未知数，3个约束方程
% 约束方程：glint在球面上，四点共面，反射定律
OG = G - O1;
a = OG * GL'/(norm(OG)*norm(GL));
b = OG * GC'/(norm(OG)*norm(GC));
solG = vpasolve([norm(OG) == R1,OG*subs(n,soln)',a == b],[xg,yg,zg],[nan,nan;nan,nan;0,inf]);
G2 = double(subs(G,solG));
% 绘制
hs1 = mesh(O1(1)+R1*cos(beta1).*cos(alpha1),O1(2)+R1*cos(beta1).*sin(alpha1),O1(3)+R1*sin(beta1));hold on
hs2 = mesh(O2(1)+R2*cos(beta2).*cos(alpha2),O2(2)+R2*cos(beta2).*sin(alpha2),O2(3)+R2*sin(beta2));hold on
axis equal
plot3(C(1),C(2),C(3),'s',L(1),L(2),L(3),'*','markersize',10),hold on
hg = plot3(G2(1),G2(2),G2(3),'+','markersize',10);hold on
hp = plot3(P(1),P(2),P(3),'o','markersize',10,'markerFaceColor','red');hold on
hl = plot3([L(1),G2(1),C(1)],[L(2),G2(2),C(2)],[L(3),G2(3),C(3)]);hold on

　　运行结果如下图，星号*表示LED，矩形表示相机，红色圆形是瞳孔。

　　

3 眼球模型知识补充

　　cornea角膜半径8mm，sclera巩膜半径为12mm。

　　两个球体之间距离范围，根据前面的球体和后面球体相交得到范围为（4,12），但是由于前面最多露出一半，进一步约束上界大概9mm，故范围（4,9）之间。取6mm。