利用插值函数求三维图像的二分点及到定曲面等距离的曲面.docx

二元函数的分片段线性插值
一函数构造
1・线性三角插值
处的函数值分别为Zl=f(xl,yl), Z2=f(x2,y2), Z3二f(x3,y3)。

我们需要求出线性函 数
P(x,y)二ax+by+c,使其满足 P(xl,yl)=Zl, P(x2,y2)=Z2, P(x3,y3)二Z3。

三角形的面积的两倍我们可以用行列式来表示，
西X 1
D= x 2 y 2 1,我们假设P(x,y) = 1,(x,y)Zj +12(x,y)Z 2 +13(x,y)Z 3, 么可以
兀3力1
构造基
函数如
2•相邻节点的二分点，三分点的x, y 坐标
所给的数据x,y,z 结构均为矩阵，如下面的图形所示，我们用矩形表示x, y 的矩阵，矩形的顶点为Pl P2 P3 P4,我们计算四个顶点的二分点。

如果我们连 接P2 P3,会得到两个三角形。

丿2(兀,刃=万
y >3
兀
1
x 2
X
显然，三个基函数分别满足插值条件:
如上图所示的三角形，已知三角形的三个顶点P1,
P2, P3
Pl Ml P2
M4 M5 M2
P3 M3 P4
二分点
3.二分点坐标的mat lab代码实现
function xxl=inter2xy(x) %计算相邻点的：分点的坐标
[m,n]=size(x);
xx=zeros(2*m-l,2*n-l); %扩人矩阵来存储节点的坐标
fori=l:m %将原矩阵在矩形顶点的x,y放在扩大后的矩阵相应处for j=l:n
xx(2*i-l ,2*j-l)=x(i,j);
end
end
for i=l:m %计算原始行相邻点屮点的x,y坐标
for j=l:n-l
xx(2*i-l ,2*j)= l/2*(xx(2*i-l ,2*j-l )+xx(2*i・ 1,2*j+1));
end
end
for i=l:m-l %计算所有列的相邻点屮点的x,y坐标
for j=l:2*n-l
xx(2*i,j)= l/2*(xx(2*i-l ,j)+xx(2*i+1 ,j));
end
end
xxl=xx; %得到二分点插值后的人矩阵
end
4.二分点插值函数的matlab代码实现
然后进行插值计算出二分点处的Z的值。

由上面所说的
P(x, y) = I】(x, y)Z, +l2(x, y)Z2 +13 (x, y)Z3即为插值函数。

Matlab代码实现如下：
function zzz=Linear(x,y,z,xx,yy) % 计算二分点插值后的Z的值
[m,n]=size(xx);
zzz=zeros(m,n); %用来存储Z值的矩阵
DO二ones(3,3); %三阶矩阵，用于基函数计算
[m,n]=size(x);
for j=l:n-l
for i=l:m-l
D0(l,l)=x(i,j); % DO即上面所说的第一种三角形
D0(l,2)二y(i,j);
D0(2,l)二x(i,j+l);
D0(2,2)=y(i,j+l);
D0(3,l)=x(i+l,j);
D0(3,2)=y(i+l,j);
D1=DO;
Dl(l,l)=xx(2*i-l,2*j); %D1 用于计算ll(x,y)
Dl(l,2)=yy(2*i-l,2*j);
D2=D0;
D2(2,1 )=xx(2*i・ 1,2*j); % D2用于计算12(x,y)
D2(2,2)=yy(2*i-l,2*j);
D3=D0;
D3(3,l)=xx(2*i-l,2*j); % D3用于计算13(x,y)
D3(3,2)=yy(2*i-l,2*j);
ll=det(Dl)/det(DO);
12=det(D2)/det(D0);
13=det(D3)/det(D0);
zzz(2*i-l ,2*j)=ll *z(i,j)+12*z(i,j+1 )+13*z(i+ l,j); %M1 的Z值
Dl(l,l)=xx(2*i,2*j);
Dl(l,2)=yy(2*i,2*j);
D2(2,l)=xx(2*i,2*j);
D2(2,2)=yy(2*i,2*j);
D3(3,l)=xx(2*i,2*j);
D3(3,2)=yy(2*i,2*j);
ll=det(Dl)/det(DO);
12=det(D2)/det(D0);
13=det(D3)/det(DO);
zzz(2*i,2*j)=Il *z(i,j)+12*z(i,j+1 )+13*z(i+1 ,j); % M5 的Z值
Dl(l,l)=xx(2*i,2*j-1);
Dl(l,2)=yy(2*i,2*j-1);
D2(2,l)=xx(2*i,2*j・l);
D2(2,2)=yy(2*i,2*j-1);
D3(3,l)=xx(2*i,2*j-1);
D3(3,2)=yy(2*i,2*j-1);
ll=det(Dl)/det(DO);
12=det(D2)/det(D0);
13=det(D3)/det(D0); zzz(2*t2*j-l)=ll*z(ij)+12*z(i,j+l)+13*z(i+l,j); % M4的Z
值 end
end
D0=ones(3,3)； for j=l:n-l for i=l:m-l
D0(l,l)=x(i,j+l); % DO 为第二种三角形
13=det(D3)/det(D0);
zzz(2*i,2*j)=l 1 *z(i,j+1 )+12*z(i+1 ,j+1 )+13*z(i+1 j); % M5 的Z 值 Dl(l,l)=xx(2*i+l,2*j); Dl(l,2)=yy(2*i+l,2*j); D2(2,l)=xx(2*i+l,2*j); D2(2,2)=yy(2*i+l,2*j); D3(3,l )=xx(2*i+1,2*j); D3(3,2)=yy(2*i+l,2*j); ll=det(Dl)/det(DO); 12=det(D2)/det(D0); 13=det(D3)/det(D0);
zzz(2*i+l,2*j)=ll*z(i,j+l)+12*z(i+l,j+l)+13*z(i+l,j); % M3的Z 值 Dl(l,l)=xx(2 气 2*j+l); Dl(l,2)=yy(2*i,2*j+1)； D2(2,l)=xx(2*i,2*j+1); D2(2,2)=yy(2*i,2*j+1); D3(3,l)=xx(2*i,2*j+1); D3(3,2)=yy(2*i,2*j+1); ll=det(Dl)/det(DO); 12=det(D2)/det(D0);
D0(l,2)=y(i,j+l); D0(2,l)=x(i+l,j+l); D0(2,2)=y(i+l,j+l); D0(3,l)二x(i+l,j);
D0(3,2)=y(i+lj); D1=DO;
Dl(l,l)=xx(2*i,2*j);
Dl(l,2)=yy(2*i,2*j); D2=D0;
D2(2,l)=xx(2*i,2*j);
D2(2,2)=yy(2*i,2*j); D3=D0;
D3(3,l)=xx(2*i,2*j);
D3(3,2)=yy(2*i,2*j); ll=det(Dl)/det(DO); 12=det(D2)/det(D0); % DI 用于计算ll(x,y) %D2用于计算12(x,y) %D3用于计算13(x,y)
13=det(D3)/det(DO);
zzz(2衍,2*j+l)=ll*z(i,j+l)+12*z(i+l,j+l)+13*z(i+l,j); % M2的Z值end
end
for i=l:m
for j=l:n zzz(2*i・l,2*j・l)=z(i,j); %原始的Z值放到扩大后的矩阵中相应的点end end
end
三利用插值点来作图分析
1 •二分点作图讨论
先利用原有数据作图，load inter2;mesh(x, y, z)。

得结果如下:
然后把我们的插值节点和原始数据得到的图结合起来，如下操作:
» load inter2
xx二inter2xy (x);
yv=inter2xy(y);
zzz=Linear (x, v, z, xx, vv);
mesh (x, v, z) ; hold on; plot3 (xx, vv, zzz, ' o')
得结果如下:
1
5、
,.
-2
0.5 .1 -°5
然后我们直接利用插值后的点来作图，mesh(xx, yy, zzz),结果如下:
-2 -2
分析插值得到的数据ZZZ,我们发现第16行的数据与上下相差较大
为此，我们再作图时，可以去除第16行，继续操作如下:
» XX(16,:)二[];
» yy(16, :) = □；
» zzz (16,:)二]];
>> mesh (xx, vv, zzz)
得到结果为：
发现与原始数据做的图比较接近，达到了插值的效果。

四到插值曲面定距离的两个相邻曲面的研究
1 •简介
一个曲而到另一个曲而的距离相等，即是说这个曲而上所有的点到另一曲而的距离相等，另外，点到面的距离是垂直的距离，而不是简单的竖直距离。

为此我们可以在已知曲面上分片，然后求出小片段曲面的法线向量，进一步得到所求曲面的相应点处的坐标，最后利用所求的点，作图即可。

需要说明的一点是，法线向量有两个，因此会有两个曲面，这在编程时需要注意。

2.算法实现
我们利用己经得到的插值节点，可以得到一个大矩阵，在这个大矩阵中，我们还是选择三角形，如P1M1M4为一个三角形，Ml P2 M5为一个三角形，同理M4 M5 P3是一个三角形，M5 M2 P4也是，利用Pl M4这个向量和Pl Ml这个向量的叉积，就可以得到这个分片得到后的小曲面的法线方向。

然后将其单位化，再乘以给定距离h,最后加上Pl处的值，就得到了所求点的坐标值。

3. matlab编程代码
function [vx,vy，vz]=vec(xx,yy,zzz) % 求相邻曲面的坐标
vl⑴=xx(i+l,j)-xx(i,j); %分片段曲面上的--个向量
vl⑵二yy(i+l,j)・yy(i,j);
v 1 (3)=zzz(i+1 ,j)-zzz(i,j);
v2(l)=xx(i,j+ l)-xx(i,j); %分片段曲面上的另一个向量
v2(2)=yy(ij+l)-yy(ij);
v2(3)=zzz(i,j+l)-zzz(i,j);
temp=cross(v 1 ,v2); %求出分片段曲面的法向量
temp 1 =sqrt(temp( 1 )A2+temp(2)A2+temp(3)A2); %单位化法向量
vx(i,j)=0.1 *temp( 1 )/templ +xx(i,j); % 相邻曲面的x坐标
vy(i,j)=O. 1 *temp(2)/temp 1 +yy(i,j); % 相邻曲面的y坐标vz(i,j)=O. 1
*temp(3)/temp 1 +zzz(i,j); % 相邻曲面的z坐标end
[m,n]=size(xx);
vx=zeros(m-l,n-l);
vy=zeros(m-l,n-l);
vz=zeros(m-1 ,n-1);
for j=l:n-l
for i=l:m-l
%相邻曲面的X坐标
%相邻曲面的y坐标
%相邻曲面的z坐标
end
end
4.作图分析
为了准确性，我们还是去掉了第16行，先利用二分点计算，操作如下:
>> load inter2
xx=inter2xy(x);
yv=inte:r2xy(y);
zzz二Linear (x, y, z, xx, yy);
» xx(16,:)二[];
» yy(16, :) = [];
» zzz (16, :) = [];
>> [vx, vv, vz] =vec (xx, vv, zzz);
>> mesh(x, v, z) ;hold on;mesh(vx, vv, vz)
得到结果如下;
利用三分点插值画图，结果如下:
得到的结果都比较理想。

最后，我们同吋画出三个曲面。

结果如下:
可见，结果比较理想。

再次证明了算法和代码的正确性。