[已解决]【高手指点】如何通过matlab画出这样的球面等高线图

nizp1982

三维需要x、y、z的坐标，还需要个函数值，即第四变量

nwcwww

大概描述下算法，为了看得清楚以sphere(5)为例，毕竟点少。

1. n=5;
   
2. [x,y,z]=sphere(n);
   
3. data=y+x.*z;
   
4. surf(x,y,z,'cdata',data,'facecolor','interp','facelighting','phong');
   
5. axis equal
   
6. hold on
   
7. 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或者基本的矩阵操作都没问题。比如说扫描横边大致可以这样：

1.     gtc= data>ciso; %csio=0.5
   
2.     leftlist=logical([diff(gtc,1,2),zeros(n+1,1)]);
   
3.     left=data(leftlist);
   
4. 
   
5.     rightlist=logical([zeros(n+1,1),diff(gtc,1,2)]);
   
6.     right=data(rightlist);
   
7.    
   
8.     lambda_lr=(repmat(ciso,size(left))-left)./(right-left); %mm*1
   
9.     x_lr=x(leftlist)+lambda_lr.*(x(rightlist)-x(leftlist));
   
10.     y_lr=y(leftlist)+lambda_lr.*(y(rightlist)-y(leftlist));
    
11.     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之类的软件更方便些。

liuyalong008

回复 22# nwcwww

测试这组数据存在问题

1. %% data generate
   
2. clear
   
3. n=100;
   
4. [x,y,z]=sphere(n);
   
5. data=y+x.*z;
   
6. data(30:80,30:60)=0.1;
   
7. surf(x,y,z,'cdata',data,'facecolor','interp','facelighting','phong','edgecolor','none');
   
8. axis equal
   
9. hold on
   
10. colormap(jet(10));
    
11. gtc= data>0.4
    
12. %%horizon
    
13. leftlist=logical([diff(gtc,1,2),zeros(n+1,1)]);
    
14. left=data(leftlist);
    
15. rightlist=logical([zeros(n+1,1),diff(gtc,1,2)]);
    
16. right=data(rightlist);
    
17. lambda_lr=(repmat(ciso,size(left))-left)./(right-left); %mm*1
    
18. x_lr=x(leftlist)+lambda_lr.*(x(rightlist)-x(leftlist));
    
19. y_lr=y(leftlist)+lambda_lr.*(y(rightlist)-y(leftlist));
    
20. z_lr=z(leftlist)+lambda_lr.*(z(rightlist)-z(leftlist));
    
21. %%vertical
    
22. uplist=logical([diff(gtc,1,1);zeros(1,n+1)])
    
23. bottomlist=logical([zeros(1,n+1);diff(gtc,1,1)])
    
24. up=data(uplist)
    
25. bottom=data(bottomlist)
    
26. lambda_ub=(repmat(ciso,size(up))-up)./(bottom-up); %mm*1
    
27. x_ub=x(uplist)+lambda_ub.*(x(bottomlist)-x(uplist));
    
28. y_ub=y(uplist)+lambda_ub.*(y(bottomlist)-y(uplist));
    
29. z_ub=z(uplist)+lambda_ub.*(z(bottomlist)-z(uplist));
    
30. %%plot
    
31. plot3(x_ub,y_ub,z_ub,'b*')
    
32. plot3(x_lr,y_lr,z_lr,'ro')
    
33. 
    
34. 
    
35. 
    

复制代码

nwcwww

回复 23# liuyalong008


我现在手头没MATLAB。具体问题是什么？

liuyalong008

回复 24# nwcwww

sorry, 我自己的问题，:handshake

nwcwww

回复 25# liuyalong008


:handshake

刚才试了下，用isosurface也可以画，而且省事不少。以前还真没用过这个命令。

1. clear
   
2. n=100;
   
3. [x,y,z]=sphere(n);
   
4. data=y+x.*z;
   
5. data(30:80,30:60)=0.1;
   
6. surf(x,y,z,'cdata',data,'facecolor','interp','facelighting','phong','edgecolor','none');
   
7. axis equal
   
8. hold on
   
9. colormap(jet(10));
   
10. xx=reshape([x,x],[size(x),2]);
    
11. yy=reshape([y,y],[size(y),2]);
    
12. zz=reshape([z,z],[size(z),2]);
    
13. vv=reshape([data,data],[size(data),2]);
    
14. hold on;
    
15. patch(isosurface(xx,yy,zz,vv,0.4));
    
16. patch(isosurface(xx,yy,zz,vv,0.9));

复制代码

liuyalong008

回复 26# nwcwww
赞一个！perfect！

nizp1982

最好能等高线上显示等高线数据