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

二元函数的分片段线性插值
一 简介
由一元函数如y=f(x)的分段线性插值，我们知道，插值的关键就是要找到基函数，对于二元函数也一样，我们需要找到基函数，但此时的基函数，应该是x ，y 的函数，而不再是简单的x 的函数，因此我们需要找到另外的方法来将三维的曲面分为小片段，从而构造基函数。

二 基函数和插值函数的构造
1. 三维曲面的基函数
由曲面的分片段构造的基函数方法有分片段线性插值基函数，分片段样条插值基函数两大类。

其中分片段线性插值又分为三角形，矩形等等，分片段样条插值又分为B 样条，三次样条等等。

这里根据给出来的数据的特点，我们选择了分片段线性插值的三角形插值来计算相邻节点的二分点，三分点的值。

2. 线性三角插值
如上图所示的三角形，已知三角形的三个顶点P1， P2，P3处的函数值分别为Z1=f(x1,y1)，Z2=f(x2,y2)，Z3=f(x3,y3)。

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

由此得线性方程组问题：
⎪⎩⎪⎨⎧=++=++=++333
222111Z c by ax Z c by ax Z c by ax 显然该方程组有唯一解的充要条件是系数矩阵行列式不等于0，该行列式恰好是三角形面积两倍。

即只要当P1，P2，P3不共线时，三角形区域上的二元线性插值函数是存在的也是唯一的。

三角形的面积的两倍我们可以用行列式来表示，
332211321321y)Z (x,l y)Z (x,l y)Z (x,l y)P(x,,1
11
++==我们假设y y y x x x D ，那么可以
利用行列式的性质来构造基函数如下：111
1),(,1111),(,1111),(212133131232321y y y x x x D y x l y y y x x x D y x l y y y x x x D y x l === 显然，三个基函数分别满足插值条件：)3,2,1,(,1,0),(=⎩⎨⎧=≠=j i j
i j i y x l j j i 三角形区域上线性插值函数的图形是空间上的一个三角形，此空间三角形的顶点恰好是(x1,y1,Z1)，(x2,y2,Z2），(x3,y3,Z3)，由此我们就构造出来了三角形插值的基函数。

3.相邻节点的二分点，三分点的x ，y 坐标
由于所给的数据比较有特点，它的数据结构上的相邻也是空间结构上的相邻，因此我们在数据结构中求解相邻节点的二分点，三分点。

所给的数据x,y,z 结构均为矩阵，因此如下面的图形所示，我们用矩形表示x ，y 的矩阵，矩形的顶点为P1 P2 P3 P4,我们计算四个顶点的二分点，即中点的坐标值。

如果我们连接P2 P3，会有两种形状的三角形，且M5位于三角形斜边的中点处。

二分点 三分点 依照同样的道理，我们可以做出三分点。

M4,M8为M1，M11的三分点，M5,M9
为M2，M12的三分点，那么在计算时，我们需要首先计算矩形四条边上的三分点，然后才能计算矩形中的三分点。

连接P2 P3我们依然可以得到两种形状的三角形，此时对角线上有两个元素M5 M8，而且分成的两个三角形内部各有一个点,如第一种三角形内部点为M4，第二种三角形内部点为M9。

当然上面给出的矩形相当于一个二阶矩阵，类似可得[m,n]阶矩阵的情形。

从上面我们可以看出，插值完之后矩阵的阶数要变大。

因此需要一个大矩阵来存储插值节点的坐标。

要注意在大矩阵中相应点处的值还是原来小矩阵中的元素的值。

4. 二分点，三分点坐标的matlab代码实现
function xx1=inter2xy(x) %计算相邻点的二分点的坐标
[m,n]=size(x);
xx=zeros(2*m-1,2*n-1); %扩大矩阵来存储节点的坐标
for i=1:m %将原矩阵在矩形顶点的x,y放在扩大后的矩阵相应处for j=1:n
xx(2*i-1,2*j-1)=x(i,j);
end
end
for i=1:m %计算原始行相邻点中点的x,y坐标
for j=1:n-1
xx(2*i-1,2*j)=1/2*(xx(2*i-1,2*j-1)+xx(2*i-1,2*j+1));
end
end
for i=1:m-1 %计算所有列的相邻点中点的x,y坐标
for j=1:2*n-1
xx(2*i,j)=1/2*(xx(2*i-1,j)+xx(2*i+1,j));
end
end
xx1=xx; %得到二分点插值后的大矩阵
end
function xx1=inter3xy(x) %计算相邻点的三分点的坐标
[m,n]=size(x);
xx=zeros(3*m-2,3*n-2); %扩大矩阵来存储节点的坐标
for i=1:m %原始位置
for j=1:n
xx(3*i-2,3*j-2)=x(i,j);
end
end
for i=1:m %原始行位置
for j=1:n-1
xx(3*i-2,3*j-1)=(2*x(i,j)+x(i,j+1))/3;
xx(3*i-2,3*j)=(x(i,j)+2*x(i,j+1))/3;
end
end
for i=1:m-1 %原始列位置
for j=1:n
xx(3*i-1,3*j-2)=(2*x(i,j)+x(i+1,j))/3;
xx(3*i,3*j-2)=(x(i,j)+2*x(i+1,j))/3;
end
end
for i=1:m-1 %对角上元素位置
for j=1:n-1
xx(3*i-1,3*j)=(2*xx(3*i-2,3*j)+xx(3*i+1,3*j))/3;
xx(3*i,3*j-1)=(2*xx(3*i+1,3*j-1)+xx(3*i-2,3*j-1))/3;
end
end
for i=1:m-1 % 第一种三角形里的一个点
for j=1:n-1
xx(3*i-1,3*j-1)=(xx(3*i+1,3*j-1)+2*xx(3*i-2,3*j-1))/3;
end
end
for i=1:m-1 %第二种三角形里的一个点
for j=1:n-1
xx(3*i,3*j)=(2*xx(3*i+1,3*j)+xx(3*i-2,3*j))/3;
end
end
xx1=xx;
end
下面我们来举例验证上面的算法和代码正确。

如下图所示：
得到a 和aa 分别为：
这就证明了我们的算法和代码的正确性。

5. 二分点，三分点插值函数的matlab 代码实现
在计算出来相邻点的二分点和三分点的坐标值后，我们要进行插值计算出二分，三分点处的Z 的值。

由上面所说的332211y)Z (x,l y)Z (x,l y)Z (x,l y)P(x,++=即为插值函数。

Matlab 代码实现如下：
function zzz=Linear(x,y,z,xx,yy) % 计算二分点插值后的Z 的值
[m,n]=size(xx);
zzz=zeros(m,n); % 用来存储Z 值的矩阵
D0=ones(3,3); % 三阶矩阵，用于基函数计算
[m,n]=size(x);
for j=1:n-1
for i=1:m-1
D0(1,1)=x(i,j); % D0即上面所说的第一种三角形
D0(1,2)=y(i,j);
D0(2,1)=x(i,j+1);
D0(2,2)=y(i,j+1);
D0(3,1)=x(i+1,j);
D0(3,2)=y(i+1,j);
D1=D0;
D1(1,1)=xx(2*i-1,2*j); %D1用于计算l1(x,y)
D1(1,2)=yy(2*i-1,2*j);
D2(2,1)=xx(2*i-1,2*j); % D2用于计算l2(x,y)
D2(2,2)=yy(2*i-1,2*j);
D3=D0;
D3(3,1)=xx(2*i-1,2*j); % D3用于计算l3(x,y)
D3(3,2)=yy(2*i-1,2*j);
l1=det(D1)/det(D0);
l2=det(D2)/det(D0);
l3=det(D3)/det(D0);
zzz(2*i-1,2*j)=l1*z(i,j)+l2*z(i,j+1)+l3*z(i+1,j); % M1的Z值
D1(1,1)=xx(2*i,2*j);
D1(1,2)=yy(2*i,2*j);
D2(2,1)=xx(2*i,2*j);
D2(2,2)=yy(2*i,2*j);
D3(3,1)=xx(2*i,2*j);
D3(3,2)=yy(2*i,2*j);
l1=det(D1)/det(D0);
l2=det(D2)/det(D0);
l3=det(D3)/det(D0);
zzz(2*i,2*j)=l1*z(i,j)+l2*z(i,j+1)+l3*z(i+1,j); % M5的Z值
D1(1,1)=xx(2*i,2*j-1);
D1(1,2)=yy(2*i,2*j-1);
D2(2,1)=xx(2*i,2*j-1);
D2(2,2)=yy(2*i,2*j-1);
D3(3,1)=xx(2*i,2*j-1);
D3(3,2)=yy(2*i,2*j-1);
l1=det(D1)/det(D0);
l2=det(D2)/det(D0);
l3=det(D3)/det(D0);
zzz(2*i,2*j-1)=l1*z(i,j)+l2*z(i,j+1)+l3*z(i+1,j); % M4的Z值end
end
D0=ones(3,3);
for j=1:n-1
for i=1:m-1
D0(1,1)=x(i,j+1); % D0为第二种三角形
D0(1,2)=y(i,j+1);
D0(2,1)=x(i+1,j+1);
D0(2,2)=y(i+1,j+1);
D0(3,1)=x(i+1,j);
D0(3,2)=y(i+1,j);
D1(1,1)=xx(2*i,2*j); % D1用于计算l1(x,y)
D1(1,2)=yy(2*i,2*j);
D2=D0;
D2(2,1)=xx(2*i,2*j); % D2用于计算l2(x,y)
D2(2,2)=yy(2*i,2*j);
D3=D0;
D3(3,1)=xx(2*i,2*j); % D3用于计算l3(x,y)
D3(3,2)=yy(2*i,2*j);
l1=det(D1)/det(D0);
l2=det(D2)/det(D0);
l3=det(D3)/det(D0);
zzz(2*i,2*j)=l1*z(i,j+1)+l2*z(i+1,j+1)+l3*z(i+1,j); % M5的Z值
D1(1,1)=xx(2*i+1,2*j);
D1(1,2)=yy(2*i+1,2*j);
D2(2,1)=xx(2*i+1,2*j);
D2(2,2)=yy(2*i+1,2*j);
D3(3,1)=xx(2*i+1,2*j);
D3(3,2)=yy(2*i+1,2*j);
l1=det(D1)/det(D0);
l2=det(D2)/det(D0);
l3=det(D3)/det(D0);
zzz(2*i+1,2*j)=l1*z(i,j+1)+l2*z(i+1,j+1)+l3*z(i+1,j); % M3的Z值
D1(1,1)=xx(2*i,2*j+1);
D1(1,2)=yy(2*i,2*j+1);
D2(2,1)=xx(2*i,2*j+1);
D2(2,2)=yy(2*i,2*j+1);
D3(3,1)=xx(2*i,2*j+1);
D3(3,2)=yy(2*i,2*j+1);
l1=det(D1)/det(D0);
l2=det(D2)/det(D0);
l3=det(D3)/det(D0);
zzz(2*i,2*j+1)=l1*z(i,j+1)+l2*z(i+1,j+1)+l3*z(i+1,j); % M2的Z值end
end
for i=1:m
for j=1:n
zzz(2*i-1,2*j-1)=z(i,j); % 原始的Z值放到扩大后的矩阵中相应的点end
end
end
同理，三分点处的Z值代码如下：
function zzz=Linear3(x,y,z,xx,yy) %计算三分点插值后的Z的值
[m,n]=size(xx);
zzz=zeros(m,n); %用来存储Z值的矩阵
D0=ones(3,3); %三阶矩阵，用于基函数计算
[m,n]=size(x);
for j=1:n-1
for i=1:m-1
D0(1,1)=x(i,j); % D0即上面所说的第一种三角形
D0(1,2)=y(i,j);
D0(2,1)=x(i,j+1);
D0(2,2)=y(i,j+1);
D0(3,1)=x(i+1,j);
D0(3,2)=y(i+1,j);
D1=D0;
D1(1,1)=xx(3*i-2,3*j-1); %D1用于计算l1(x,y)
D1(1,2)=yy(3*i-2,3*j-1);
D2=D0;
D2(2,1)=xx(3*i-2,3*j-1); % D2用于计算l2(x,y)
D2(2,2)=yy(3*i-2,3*j-1);
D3=D0;
D3(3,1)=xx(3*i-2,3*j-1); % D3用于计算l3(x,y)
D3(3,2)=yy(3*i-2,3*j-1);
l1=det(D1)/det(D0);
l2=det(D2)/det(D0);
l3=det(D3)/det(D0);
zzz(3*i-2,3*j-1)=l1*z(i,j)+l2*z(i,j+1)+l3*z(i+1,j); % M1的Z值
D1(1,1)=xx(3*i-2,3*j);
D1(1,2)=yy(3*i-2,3*j);
D2(2,1)=xx(3*i-2,3*j);
D2(2,2)=yy(3*i-2,3*j);
D3(3,1)=xx(3*i-2,3*j);
D3(3,2)=yy(3*i-2,3*j);
l1=det(D1)/det(D0);
l2=det(D2)/det(D0);
l3=det(D3)/det(D0);
zzz(3*i-2,3*j)=l1*z(i,j)+l2*z(i,j+1)+l3*z(i+1,j); % M2的Z值
D1(1,1)=xx(3*i-1,3*j-2);
D1(1,2)=yy(3*i-1,3*j-2);
D2(2,1)=xx(3*i-1,3*j-2);
D2(2,2)=yy(3*i-1,3*j-2);
D3(3,1)=xx(3*i-1,3*j-2);
D3(3,2)=yy(3*i-1,3*j-2);
l2=det(D2)/det(D0);
l3=det(D3)/det(D0);
zzz(3*i-1,3*j-2)=l1*z(i,j)+l2*z(i,j+1)+l3*z(i+1,j); % M3的Z值D1(1,1)=xx(3*i-1,3*j-1);
D1(1,2)=yy(3*i-1,3*j-1);
D2(2,1)=xx(3*i-1,3*j-1);
D2(2,2)=yy(3*i-1,3*j-1);
D3(3,1)=xx(3*i-1,3*j-1);
D3(3,2)=yy(3*i-1,3*j-1);
l1=det(D1)/det(D0);
l2=det(D2)/det(D0);
l3=det(D3)/det(D0);
zzz(3*i-1,3*j-1)=l1*z(i,j)+l2*z(i,j+1)+l3*z(i+1,j); % M4的Z值D1(1,1)=xx(3*i-1,3*j);
D1(1,2)=yy(3*i-1,3*j);
D2(2,1)=xx(3*i-1,3*j);
D2(2,2)=yy(3*i-1,3*j);
D3(3,1)=xx(3*i-1,3*j);
D3(3,2)=yy(3*i-1,3*j);
l1=det(D1)/det(D0);
l2=det(D2)/det(D0);
l3=det(D3)/det(D0);
zzz(3*i-1,3*j)=l1*z(i,j)+l2*z(i,j+1)+l3*z(i+1,j); % M5的Z值D1(1,1)=xx(3*i,3*j-2);
D1(1,2)=yy(3*i,3*j-2);
D2(2,1)=xx(3*i,3*j-2);
D2(2,2)=yy(3*i,3*j-2);
D3(3,1)=xx(3*i,3*j-2);
D3(3,2)=yy(3*i,3*j-2);
l1=det(D1)/det(D0);
l2=det(D2)/det(D0);
l3=det(D3)/det(D0);
zzz(3*i,3*j-2)=l1*z(i,j)+l2*z(i,j+1)+l3*z(i+1,j); % M7的Z值D1(1,1)=xx(3*i,3*j-1);
D1(1,2)=yy(3*i,3*j-1);
D2(2,1)=xx(3*i,3*j-1);
D2(2,2)=yy(3*i,3*j-1);
D3(3,1)=xx(3*i,3*j-1);
D3(3,2)=yy(3*i,3*j-1);
l2=det(D2)/det(D0);
l3=det(D3)/det(D0);
zzz(3*i,3*j-1)=l1*z(i,j)+l2*z(i,j+1)+l3*z(i+1,j); % M8的Z值end
end
D0=ones(3,3);
for j=1:n-1
for i=1:m-1
D0(1,1)=x(i,j+1); % D0即上面所说的第一种三角形
D0(1,2)=y(i,j+1);
D0(2,1)=x(i+1,j+1);
D0(2,2)=y(i+1,j+1);
D0(3,1)=x(i+1,j);
D0(3,2)=y(i+1,j);
D1=D0;
D1(1,1)=xx(3*i,3*j);
D1(1,2)=yy(3*i,3*j);
D2=D0;
D2(2,1)=xx(3*i,3*j);
D2(2,2)=yy(3*i,3*j);
D3=D0;
D3(3,1)=xx(3*i,3*j);
D3(3,2)=yy(3*i,3*j);
l1=det(D1)/det(D0);
l2=det(D2)/det(D0);
l3=det(D3)/det(D0);
zzz(3*i,3*j)=l1*z(i,j+1)+l2*z(i+1,j+1)+l3*z(i+1,j); % M9的Z值end
end
for i=1:m
for j=1:n
zzz(3*i-2,3*j-2)=z(i,j); % 原始的Z值放到扩大后的矩阵中相应的点end
end
end
算法的正确性通过下一步验证。

三利用插值点来作图分析
1.二分点作图讨论
先利用原有数据作图，load inter2;mesh(x,y,z)。

得结果如下：
然后把我们的插值节点和原始数据得到的图结合起来，如下操作：
得结果如下：
插值得到的值基本正确，有个别点插值可能不正确，如上图我们发现有几个点跑到了所给曲面的外侧，没在曲面上。

我们直接利用插值后的点来作图，mesh(xx,yy,zzz)，结果如下：
我们发现它只有在某一列插值较差，其他点处插值较好。

我们分析可能是因为在马鞍面的中心处由于旋转度较大，导致P1 P2 P3三点几乎共线，所以得到结果不准确。

分析插值得到的数据zzz，我们发现第16行的数据与上下相差较大
为此，我们再作图时，可以去除第16行，继续操作如下：
得到结果为：
发现与原始数据做的图比较接近，达到了插值的效果。

2.三分点作图讨论
先操作如下：
这次我们先观察数据结构，
我们发现得到的zzz的第24,26,27,46行和第最后一列的数据和上下相差太大，这可能是因为利用的插值方法插值没有成功所导致，也可能是因为所给的图形的特点所导致结果如下。

为此我们去掉上述的四行三列，并作图，操作如下：
得到了理想的结果，我们将三分点的值和原始数据得到的图结合起来观察：
得到结果：
去除一些行和列后发现插值点都落在了曲面上，说明了算法和代码的正确性。

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

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

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

2.算法实现
我们利用已经得到的插值节点，可以得到一个大矩阵，在这个大矩阵中，我们还是选择三角形，如P1 M1 M4为一个三角形，M1 P2 M5为一个三角形，同理M4 M5 P3是一个三角形，M5 M2 P4也是，它们都是同一种形状的三角形，因此我们就选取这样的三角形。

利用
P1 M4这个向量和P1 M1这个向量的叉积，就可以得到这个分片得到后的小曲面的法线方向。

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

对于三分点插值后的结果同样可以利用此法得到相邻曲面。

3.matlab编程代码
function [vx,vy,vz]=vec(xx,yy,zzz) % 求相邻曲面的坐标
[m,n]=size(xx);
vx=zeros(m-1,n-1); % 相邻曲面的x坐标
vy=zeros(m-1,n-1); % 相邻曲面的y坐标
vz=zeros(m-1,n-1); % 相邻曲面的z坐标
for j=1:n-1
for i=1:m-1
v1(1)=xx(i+1,j)-xx(i,j); % 分片段曲面上的一个向量
v1(2)=yy(i+1,j)-yy(i,j);
v1(3)=zzz(i+1,j)-zzz(i,j);
v2(1)=xx(i,j+1)-xx(i,j); % 分片段曲面上的另一个向量
v2(2)=yy(i,j+1)-yy(i,j);
v2(3)=zzz(i,j+1)-zzz(i,j);
temp=cross(v1,v2); % 求出分片段曲面的法向量
temp1=sqrt(temp(1)^2+temp(2)^2+temp(3)^2); % 单位化法向量
vx(i,j)=0.1*temp(1)/temp1+xx(i,j); % 相邻曲面的x坐标
vy(i,j)=0.1*temp(2)/temp1+yy(i,j); % 相邻曲面的y坐标
vz(i,j)=0.1*temp(3)/temp1+zzz(i,j); % 相邻曲面的z坐标
end
end
end
5.作图分析
为了准确性，我们还是去掉了第16行，先利用二分点计算，操作如下：
得到结果如下；
利用三分点插值画图，代码略，结果如下：
得到的结果都比较理想。

最后，我们改变法向量的方向，(添一个负号就可以了），先保留着原来的两个曲面，改完后，再操作如下：
可以得到另一个相邻曲面。

这样，就得到了三个曲面的图，结果如下：
可见，结果比较理想。

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