利用插值函数求三维图像的二分点三分点及其到定曲面等距离的曲面.docx
《利用插值函数求三维图像的二分点三分点及其到定曲面等距离的曲面.docx》由会员分享，可在线阅读，更多相关《利用插值函数求三维图像的二分点三分点及其到定曲面等距离的曲面.docx（23页珍藏版）》请在冰豆网上搜索。
利用插值函数求三维图像的二分点三分点及其到定曲面等距离的曲面
二元函数的分片段线性插值
一简介
由一元函数如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。
由此得线性方程组问题：
显然该方程组有唯一解的充要条件是系数矩阵行列式不等于0，该行列式恰好是三角形面积两倍。
即只要当P1，P2，P3不共线时，三角形区域上的二元线性插值函数是存在的也是唯一的。
三角形的面积的两倍我们可以用行列式来表示，
，那么可以利用行列式的性质来构造基函数如下：
显然，三个基函数分别满足插值条件：
三角形区域上线性插值函数的图形是空间上的一个三角形，此空间三角形的顶点恰好是（x1,y1,Z1），（x2,y2,Z2），（x3,y3,Z3），由此我们就构造出来了三角形插值的基函数。
3.相邻节点的二分点，三分点的x，y坐标
由于所给的数据比较有特点，它的数据结构上的相邻也是空间结构上的相邻，因此我们在数据结构中求解相邻节点的二分点，三分点。
所给的数据x,y,z结构均为矩阵，因此如下面的图形所示，我们用矩形表示x，y的矩阵，矩形的顶点为P1P2P3P4,我们计算四个顶点的二分点，即中点的坐标值。
如果我们连接P2P3，会有两种形状的三角形，且M5位于三角形斜边的中点处。
二分点三分点
依照同样的道理，我们可以做出三分点。
M4,M8为M1，M11的三分点，M5,M9为M2，M12的三分点，那么在计算时，我们需要首先计算矩形四条边上的三分点，然后才能计算矩形中的三分点。
连接P2P3我们依然可以得到两种形状的三角形，此时对角线上有两个元素M5M8，而且分成的两个三角形内部各有一个点,如第一种三角形内部点为M4，第二种三角形内部点为M9。
当然上面给出的矩形相当于一个二阶矩阵，类似可得[m,n]阶矩阵的情形。
从上面我们可以看出，插值完之后矩阵的阶数要变大。
因此需要一个大矩阵来存储插值节点的坐标。
要注意在大矩阵中相应点处的值还是原来小矩阵中的元素的值。
4.二分点，三分点坐标的matlab代码实现
functionxx1=inter2xy（x）%计算相邻点的二分点的坐标
[m,n]=size（x）;
xx=zeros（2*m-1,2*n-1）;%扩大矩阵来存储节点的坐标
fori=1:
m%将原矩阵在矩形顶点的x,y放在扩大后的矩阵相应处
forj=1:
n
xx（2*i-1,2*j-1）=x（i,j）;
end
end
fori=1:
m%计算原始行相邻点中点的x,y坐标
forj=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
fori=1:
m-1%计算所有列的相邻点中点的x,y坐标
forj=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
functionxx1=inter3xy（x）%计算相邻点的三分点的坐标
[m,n]=size（x）;
xx=zeros（3*m-2,3*n-2）;%扩大矩阵来存储节点的坐标
fori=1:
m%原始位置
forj=1:
n
xx（3*i-2,3*j-2）=x（i,j）;
end
end
fori=1:
m%原始行位置
forj=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
fori=1:
m-1%原始列位置
forj=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
fori=1:
m-1%对角上元素位置
forj=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
fori=1:
m-1%第一种三角形里的一个点
forj=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
fori=1:
m-1%第二种三角形里的一个点
forj=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的值。
由上面所说的
即为插值函数。
Matlab代码实现如下：
functionzzz=Linear（x,y,z,xx,yy）%计算二分点插值后的Z的值
[m,n]=size（xx）;
zzz=zeros（m,n）;%用来存储Z值的矩阵
D0=ones（3,3）;%三阶矩阵，用于基函数计算
[m,n]=size（x）;
forj=1:
n-1
fori=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=D0;
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）;
forj=1:
n-1
fori=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（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
fori=1:
m
forj=1:
n
zzz（2*i-1,2*j-1）=z（i,j）;%原始的Z值放到扩大后的矩阵中相应的点
end
end
end
同理，三分点处的Z值代码如下：
functionzzz=Linear3（x,y,z,xx,yy）%计算三分点插值后的Z的值
[m,n]=size（xx）;
zzz=zeros（m,n）;%用来存储Z值的矩阵
D0=ones（3,3）;%三阶矩阵，用于基函数计算
[m,n]=size（x）;
forj=1:
n-1
fori=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）;
l1=det（D1）/det（D0）;
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）;
l1=det（D1）/det（D0）;
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）;
forj=1:
n-1
fori=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
fori=1:
m
forj=1:
n
zzz（3*i-2,3*j-2）=z（i,j）;%原始的Z值放到扩大后的矩阵中相应的点
end
end
end
算法的正确性通过下一步验证。
三利用插值点来作图分析
1.二分点作图讨论
先利用原有数据作图，loadinter2;mesh（x,y,z）。
得结果如下：
可知，这是一个马鞍面。
然后把我们的插值节点和原始数据得到的图结合起来，如下操作：
得结果如下：
插值得到的值基本正确，有个别点插值可能不正确，如上图我们发现有几个点跑到了所给曲面的外侧，没在曲面上。
我们直接利用插值后的点来作图，mesh（xx,yy,zzz），结果如下：
我们发现它只有在某一列插值较差，其他点处插值较好。
我们分析可能是因为在马鞍面的中心处由于旋转度较大，导致P1P2P3三点几乎共线，所以得到结果不准确。
分析插值得到的数据zzz，我们发现第16行的数据与上下相差较大
为此，我们再作图时，可以去除第16行，继续操作如下：
得到结果为：
发现与原始数据做的图比较接近，达到了插值的效果。
2.三分点作图讨论
先操作如下：
这次我们先观察数据结构，
我们发现得到的zzz的第24,26,27,46行和第最后一列的数据和上下相差太大，这可能是因为利用的插值方法插值没有成功所导致，也可能是因为所给的图形的特点所导致结果如下。
为此我们去掉上述的四行三列，并作图，操作如下：
得到了理想的结果，我们将三分点的值和原始数据得到的图结合起来观察：
得到结果：
去除一些行和列后发现插值点都落在了曲面上，说明了算法和代码的正确性。
四到插值曲面定距离的两个相邻曲面的研究
1.简介
一个曲面到另一个曲面的距离相等，即是说这个曲面上所有的点到另一曲面的距离相等，另外，点到面的距离是垂直的距离，而不是简单的竖直距离。
为此我们可以在已知曲面上分片，然后求出小片段曲面的法线向量，进一步得到所求曲面的相应点处的坐标，最后利用所求的点，作图即可。
需要说明的一点是，法线向量有两个，因此会有两个曲面，这在编程时需要注意。
2.算法实现
我们利用已经得到的插值节点，可以得到一个大矩阵，在这个大矩阵中，我们还是选择三角形，如P1M1M4为一个三角形，M1P2M5为一个三角形，同理M4M5P3是一个三角形，M5M2P4也是，它们都是同一种形状的三角形，因此我们就选取这样的三角形。
利用
P1M4这个向量和P1M1这个向量的叉积，就可以得到这个分片得到后的小曲面的法线方向。
然后将其单位化，再乘以给定距离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坐标
forj=1:
n-1
fori=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行，先利用二分点计算，操作如下：
得到结果如下；
利用三分点插值画图，代码略，结果如下：
得到的结果都比较理想。
最后，我们改变法向量的方向，（添一个负号就可以了），先保留着原来的两个曲面，改完后，再操作如下：
可以得到另一个相邻曲面。
这样，就得到了三个曲面的图，结果如下：
可见，结果比较理想。
再次证明了算法和代码的正确性。
升级会员 