大概描述下算法,为了看得清楚以sphere(5)为例,毕竟点少。n=5;
=sphere(n);
data=y+x.*z;
surf(x,y,z,'cdata',data,'facecolor','interp','facelighting','phong');
axis equal
hold on
colormap(jet(10));效果如图:
这里x,y,z,data都是(n+1)*(n+1)的矩阵。
就最一般的情况而言,点阵中相邻的四点P(i,j),P(i,j+1),P(i+1,j+1),P(i+1,j)会组成一个四边形。所以我们可以通过在矩阵的行方向或列方向上进行相关操作而对边进行插值。
比如说我们要画data=0.5的等高线,而
data(1,2)= 0.8 data(2,2)=0.4
data(1,1)=0.4 data(2,1)=0.3
那我们就可以用定比分公式,在这个四边形的上边和左边分别得到定比分值lambda,并用这个值和两个端点的坐标值计算出插值点的x,y,z.
这个方法也是contourc的算法,在matlab的doc文件中有描述。
扫描所有的横边(改变矩阵的列数)和纵边(改变矩阵的行数),如果某条边的端点正好对应目标值0.5又或者两端点的data数值跨越了0.5,我们就可以在该边上得到相应的数据点。这里的纵横是针对矩阵下标来说的。
扫描这些边的方法不止一种,用gradient filter或者基本的矩阵操作都没问题。比如说扫描横边大致可以这样: gtc= data>ciso; %csio=0.5
leftlist=logical();
left=data(leftlist);
rightlist=logical();
right=data(rightlist);
lambda_lr=(repmat(ciso,size(left))-left)./(right-left); %mm*1
x_lr=x(leftlist)+lambda_lr.*(x(rightlist)-x(leftlist));
y_lr=y(leftlist)+lambda_lr.*(y(rightlist)-y(leftlist));
z_lr=z(leftlist)+lambda_lr.*(z(rightlist)-z(leftlist));x_lr, y_lr, z_lr就会是在横边上得到的插值点。纵边上插值点的寻找方法类似。
sphere(5),ciso=0.5对应的插值点如下图。显示方便纵边上的点标为蓝色,横边的标为红色。
当sphere的点比较多的时候当然更准确,看起来也舒服许多:
最后把所有点排序,用plot3画出来再加上文字标注就行了。
当然这只是比较理想的情况。
从外部导来的数据往往还需要进行其他处理,比如说散点的储存顺序往往和sphere的不同,球面上相邻的点可能在矩阵中并不相邻。这就得先把那些点fit对应到sphere上。
所以应该还是TECPLOT之类的软件更方便些。 回复 22# nwcwww
测试这组数据存在问题%% data generate
clear
n=100;
=sphere(n);
data=y+x.*z;
data(30:80,30:60)=0.1;
surf(x,y,z,'cdata',data,'facecolor','interp','facelighting','phong','edgecolor','none');
axis equal
hold on
colormap(jet(10));
gtc= data>0.4
%%horizon
leftlist=logical();
left=data(leftlist);
rightlist=logical();
right=data(rightlist);
lambda_lr=(repmat(ciso,size(left))-left)./(right-left); %mm*1
x_lr=x(leftlist)+lambda_lr.*(x(rightlist)-x(leftlist));
y_lr=y(leftlist)+lambda_lr.*(y(rightlist)-y(leftlist));
z_lr=z(leftlist)+lambda_lr.*(z(rightlist)-z(leftlist));
%%vertical
uplist=logical()
bottomlist=logical()
up=data(uplist)
bottom=data(bottomlist)
lambda_ub=(repmat(ciso,size(up))-up)./(bottom-up); %mm*1
x_ub=x(uplist)+lambda_ub.*(x(bottomlist)-x(uplist));
y_ub=y(uplist)+lambda_ub.*(y(bottomlist)-y(uplist));
z_ub=z(uplist)+lambda_ub.*(z(bottomlist)-z(uplist));
%%plot
plot3(x_ub,y_ub,z_ub,'b*')
plot3(x_lr,y_lr,z_lr,'ro')
回复 23# liuyalong008
我现在手头没MATLAB。具体问题是什么? 回复 24# nwcwww
sorry, 我自己的问题,:handshake 回复 25# liuyalong008
:handshake
刚才试了下,用isosurface也可以画,而且省事不少。以前还真没用过这个命令。clear
n=100;
=sphere(n);
data=y+x.*z;
data(30:80,30:60)=0.1;
surf(x,y,z,'cdata',data,'facecolor','interp','facelighting','phong','edgecolor','none');
axis equal
hold on
colormap(jet(10));
xx=reshape(,);
yy=reshape(,);
zz=reshape(,);
vv=reshape(,);
hold on;
patch(isosurface(xx,yy,zz,vv,0.4));
patch(isosurface(xx,yy,zz,vv,0.9)); 回复 26# nwcwww
赞一个!perfect! 最好能等高线上显示等高线数据
页:
1
[2]