曲面拟合之最小二乘法（矩形域）

矩形阈的最小二乘曲面拟合

算法实现：

#include "stdlib.h"#include "math.h"void pir2(x,y,z,n,m,a,p,q,dt)int n,m,p,q;double x[],y[],z[],a[],dt[];{ int i,j,k,l,kk;double apx[20],apy[20],bx[20],by[20],u[20][20];double t[20],t1[20],t2[20],xx,yy,d1,d2,g,g1,g2;double x2,dd,y1,x1,*v;v=malloc(20*m*sizeof(double));for (i=0; i<=p-1; i++){ l=i*q;for (j=0; j<=q-1; j++) a[l+j]=0.0;}if (p>n) p=n;if (p>20) p=20;if (q>m) q=m;if (q>20) q=20;xx=0.0;for (i=0; i<=n-1; i++)xx=xx+x[i]/(1.0*n);yy=0.0;for (i=0; i<=m-1; i++)yy=yy+y[i]/(1.0*m);d1=1.0*n; apx[0]=0.0;for (i=0; i<=n-1; i++)apx[0]=apx[0]+x[i]-xx;apx[0]=apx[0]/d1;for (j=0; j<=m-1; j++){ v[j]=0.0;for (i=0; i<=n-1; i++)v[j]=v[j]+z[i*m+j];v[j]=v[j]/d1;}if (p>1){ d2=0.0; apx[1]=0.0;for (i=0; i<=n-1; i++){ g=x[i]-xx-apx[0];d2=d2+g*g;apx[1]=apx[1]+(x[i]-xx)*g*g;}apx[1]=apx[1]/d2;bx[1]=d2/d1;for (j=0; j<=m-1; j++){ v[m+j]=0.0;for (i=0; i<=n-1; i++){ g=x[i]-xx-apx[0];v[m+j]=v[m+j]+z[i*m+j]*g;}v[m+j]=v[m+j]/d2;}d1=d2;}for (k=2; k<=p-1; k++){ d2=0.0; apx[k]=0.0;for (j=0; j<=m-1; j++) v[k*m+j]=0.0;for (i=0; i<=n-1; i++){ g1=1.0; g2=x[i]-xx-apx[0];for (j=2; j<=k; j++){ g=(x[i]-xx-apx[j-1])*g2-bx[j-1]*g1;g1=g2; g2=g;}d2=d2+g*g;apx[k]=apx[k]+(x[i]-xx)*g*g;for (j=0; j<=m-1; j++)v[k*m+j]=v[k*m+j]+z[i*m+j]*g;}for (j=0; j<=m-1; j++)v[k*m+j]=v[k*m+j]/d2;apx[k]=apx[k]/d2;bx[k]=d2/d1;d1=d2;}d1=m; apy[0]=0.0;for (i=0; i<=m-1; i++)apy[0]=apy[0]+y[i]-yy;apy[0]=apy[0]/d1;for (j=0; j<=p-1; j++){ u[j][0]=0.0;for (i=0; i<=m-1; i++)u[j][0]=u[j][0]+v[j*m+i];u[j][0]=u[j][0]/d1;}if (q>1){ d2=0.0; apy[1]=0.0;for (i=0; i<=m-1; i++){ g=y[i]-yy-apy[0];d2=d2+g*g;apy[1]=apy[1]+(y[i]-yy)*g*g;}apy[1]=apy[1]/d2;by[1]=d2/d1;for (j=0; j<=p-1; j++){ u[j][1]=0.0;for (i=0; i<=m-1; i++){ g=y[i]-yy-apy[0];u[j][1]=u[j][1]+v[j*m+i]*g;}u[j][1]=u[j][1]/d2;}d1=d2;}for (k=2; k<=q-1; k++){ d2=0.0; apy[k]=0.0;for (j=0; j<=p-1; j++) u[j][k]=0.0;for (i=0; i<=m-1; i++){ g1=1.0;g2=y[i]-yy-apy[0];for (j=2; j<=k; j++){ g=(y[i]-yy-apy[j-1])*g2-by[j-1]*g1;g1=g2; g2=g;}d2=d2+g*g;apy[k]=apy[k]+(y[i]-yy)*g*g;for (j=0; j<=p-1; j++)u[j][k]=u[j][k]+v[j*m+i]*g;}for (j=0; j<=p-1; j++)u[j][k]=u[j][k]/d2;apy[k]=apy[k]/d2;by[k]=d2/d1;d1=d2;}v[0]=1.0; v[m]=-apy[0]; v[m+1]=1.0;for (i=0; i<=p-1; i++)for (j=0; j<=q-1; j++)a[i*q+j]=0.0;for (i=2; i<=q-1; i++){ v[i*m+i]=v[(i-1)*m+(i-1)];v[i*m+i-1]=-apy[i-1]*v[(i-1)*m+i-1]+v[(i-1)*m+i-2];if (i>=3)for (k=i-2; k>=1; k--)v[i*m+k]=-apy[i-1]*v[(i-1)*m+k]+v[(i-1)*m+k-1]-by[i-1]*v[(i-2)*m+k];v[i*m]=-apy[i-1]*v[(i-1)*m]-by[i-1]*v[(i-2)*m];}for (i=0; i<=p-1; i++){ if (i==0) { t[0]=1.0; t1[0]=1.0;}else{ if (i==1){ t[0]=-apx[0]; t[1]=1.0;t2[0]=t[0]; t2[1]=t[1];}else{ t[i]=t2[i-1];t[i-1]=-apx[i-1]*t2[i-1]+t2[i-2];if (i>=3)for (k=i-2; k>=1; k--)t[k]=-apx[i-1]*t2[k]+t2[k-1]-bx[i-1]*t1[k];t[0]=-apx[i-1]*t2[0]-bx[i-1]*t1[0];t2[i]=t[i];for (k=i-1; k>=0; k--){ t1[k]=t2[k]; t2[k]=t[k];}}}for (j=0; j<=q-1; j++)for (k=i; k>=0; k--)for (l=j; l>=0; l--)a[k*q+l]=a[k*q+l]+u[i][j]*t[k]*v[j*m+l];}dt[0]=0.0; dt[1]=0.0; dt[2]=0.0;for (i=0; i<=n-1; i++){ x1=x[i]-xx;for (j=0; j<=m-1; j++){ y1=y[j]-yy;x2=1.0; dd=0.0;for (k=0; k<=p-1; k++){ g=a[k*q+q-1];for (kk=q-2; kk>=0; kk--)g=g*y1+a[k*q+kk];g=g*x2; dd=dd+g; x2=x2*x1;}dd=dd-z[i*m+j];if (fabs(dd)>dt[2]) dt[2]=fabs(dd);dt[0]=dt[0]+dd*dd;dt[1]=dt[1]+fabs(dd);}}free(v);return;}

例子：

解答：

#include "math.h"#include "stdio.h"#include "8pir2.c"main(){ int i,j;double x[11],y[21],z[11][21],a[6][5],dt[3];for (i=0; i<=10; i++) x[i]=0.2*i;for (i=0; i<=20; i++) y[i]=0.1*i;for (i=0; i<=10; i++)for (j=0; j<=20; j++)z[i][j]=exp(x[i]*x[i]-y[j]*y[j]);pir2(x,y,z,11,21,a,6,5,dt);printf("\n");printf("MAT A(i,j) IS:\n");for (i=0; i<=5; i++){ for (j=0; j<=4; j++)printf("%e  ",a[i][j]);printf("\n");}printf("\n");for (i=0; i<=2; i++)printf("dt(%2d)=%e  ",i,dt[i]);printf("\n");}

结果：

曲面拟合之最小二乘法（矩形域）相关推荐

matlab 二维矩形函数,科学网—利用MATLAB对非矩形域实现二维插值 - 张乐乐的博文...
>> load('x1.mat'); >> load('y1.mat') >> load('T.mat'); >> load('Lon.mat'); & ...
matlab 二维插值验证,科学网-利用MATLAB对非矩形域实现二维插值-张乐乐的博文...
>> load('x1.mat'); >> load('y1.mat') >> load('T.mat'); >> load('Lon.mat'); & ...
移动最小二乘法（MLS）曲线曲面拟合C++代码实现
移动最小二乘法(MLS)曲线曲面拟合曲线曲面拟合有很多种方法,Beizer,B样条等,曲面拟合移动最小二乘法是一个很好的选择,本文会详细讲解一下移动最小二乘法方法拟合曲面,并给出C++代码实现. 本 ...
matlab mpt工具箱帮助文档_替代 Matlab 的国产软件出现？
文章来源:ithome.雷锋网亡羊补牢,犹未晚也.近日,哈工大.哈工程 Matlab 被禁一事引起了各方科研人员的注意.不少专业人士表示 Matlab 被禁是意料之中,但如何找到一款替代品却是一大难 ...
知乎热议：替代 Matlab 的国产软件出现，开发商称半年内实现 Matlab 功能的70%
点击上方"3D视觉工坊",选择"星标" 干货第一时间送达整理:公众号@雷锋网本文仅做学术分享,如有侵权,请联系删除. 谁能替代 Matlab? " ...
号称可替代MATLAB的国产软件来了！
本文转载自:新智元来源:知乎 | 编辑:舒婷.白峰 [导读]「MATLAB热」整整持续了两周,关于MATLAB被禁之后中国科研人员如何应对引发了热烈的讨论.近日,国产软件TRUFFE ...
三维扫描原理及精度控制
三维扫描学习目录一.理论基础 1. 三维扫描原理及精度控制二.边缘定位(原理) 2. 边缘细定位边缘(求解普遍亚像素边缘) 3. 针对圆型标志点曲率滤波三.求解标志点圆心 4. 三种基于矩的亚像 ...
Matlab中griddata函数拟合三维散点
griddata可以插入二维或三维散点数据 griddata有以下三种形式: ①vq = griddata(x,y,v,xq,yq) ②vq = griddata(x,y,z,v,xq,yq,zq) ...
ASP.NET网站建设基本常用代码
1.为按钮添加确认对话框 Button.Attributes.Add("onclick","return confirm('确认?')"); Button.At ...

曲面拟合之最小二乘法（矩形域）

矩形阈的最小二乘曲面拟合

曲面拟合之最小二乘法（矩形域）相关推荐

最新文章

热门文章