Mathematica中文论坛-非官方

标题: 豆粑粑 matlab contour3D 散点的点密度 3D等值面 [打印本页]

---

作者: meatball1982    时间: 2017-3-16 15:25
标题: 豆粑粑 matlab contour3D 散点的点密度 3D等值面


空间散点，根据点密度画图。在3维空间中等高值曲面显示。
这个问题经常被问到，其它软件也有相似的实现。
记录下来，给以后的我\footnote{如果能活到再次使用这类图的时候}以帮助。
数据,m文件附件上传。


cv_357.txt.tar (271 KB, 下载次数: 2)

2017-3-16 15:23 上传

点击文件名下载附件
data

work_contour_3D.m.tar (4.5 KB, 下载次数: 1)

2017-3-16 15:24 上传

点击文件名下载附件
mfile

复制代码

1. clear all
   
2. clc
   
3. close all
   
4. 
   
5. %% outline
   
6. % plot scatter points, surf, countour3D
   
7. 
   
8. %% main
   
9. % load data ---------------------------------------------------------------
   
10. dat=load('cv_357.txt');
    
11. n=length(dat(:,1));
    
12. 
    
13. x=dat(:,1);
    
14. y=dat(:,2);
    
15. z=dat(:,3);
    
16. 
    
17. % grid size ---------------------------------------------------------------
    
18. n_cv=30;
    
19. cv_g=linspace(-pi,pi,n_cv);
    
20. 
    
21. cv_int=0.5; % 2 calculate point denisty, half of cubic size
    
22. cv_x=cv_g;
    
23. cv_y=cv_g;
    
24. cv_z=cv_g;
    
25. 
    
26. [cv_X,cv_Y,cv_Z]=meshgrid(cv_g,cv_g,cv_g);
    
27. 
    
28. % points denisty ----------------------------------------------------------
    
29. for i=1:n_cv
    
30.     ind_x= (x>=cv_x(i)-cv_int) & (x<cv_x(i)+cv_int);
    
31.     for j=1:n_cv
    
32.         ind_y= (y>=cv_y(j)-cv_int) & (y<cv_y(j)+cv_int);
    
33.         for k=1:n_cv
    
34.             ind_z= (z>=cv_z(k)-cv_int) & (z<cv_z(k)+cv_int);
    
35.             p_n = sum(ind_x & ind_y & ind_z);
    
36.             cv_V(j,i,k)= p_n;
    
37.         end
    
38.     end
    
39. end
    
40. 
    
41. % %% 4 conv
    
42. % save Mat_post_3D_contour.mat
    
43. %
    
44. % load Mat_post_3D_contour.mat
    
45. % 2 smooth the points density, change 1 to 3, or 5 or 7 according to your
    
46. % data
    
47. cv_V=smooth3(cv_V,'box',1);
    
48. 
    
49. %% original points --------------------------------------------------------
    
50. figure(1)
    
51. plot3(x,y,z,'.')
    
52. axis([-pi pi -pi pi -pi pi])
    
53. 
    
54. axis equal tight
    
55. view(20,20)
    
56. box on
    
57. ax=gca;
    
58. ax.BoxStyle = 'full';
    
59. 
    
60. %% surf version -----------------------------------------------------------
    
61. 
    
62. figure(2)
    
63. hold on
    
64. for i=1:3:n_cv
    
65.     surf( cv_X(:,:,i),cv_Y(:,:,i),cv_Z(:,:,i),cv_V(:,:,i),...
    
66.         'edgecolor','none','facealpha',0.05+0.01*i)
    
67. end
    
68. axis([-pi pi -pi pi -pi pi])
    
69. colormap(flipud(jet))
    
70. 
    
71. axis equal tight
    
72. view(20,20)
    
73. box on
    
74. ax=gca;
    
75. ax.BoxStyle = 'full';
    
76. 
    
77. %% contour3D --------------------------------------------------------------
    
78. figure(3)
    
79. n_layer=9; % layer number
    
80. col_mm=flipud(jet(n_layer)); % color map
    
81. 
    
82. lev_add=2.3.^[2:n_layer+1]; % which layer you want to plot
    
83. lev_beg=80;                 % begin layer
    
84. 
    
85. % contour 3D with each layer
    
86. for i = 1: n_layer
    
87.     mm_lev=lev_beg+lev_add(i);
    
88.     p1 = patch(isosurface(cv_V,mm_lev),'FaceColor',col_mm(i,:),...
    
89.     'EdgeColor','none','FaceAlpha',0.1+0.015*i);
    
90.     isonormals(cv_V,p1);
    
91. end
    
92. % other setting
    
93. axis equal
    
94. view(20,20);
    
95. axis vis3d equal tight
    
96. camlight headlight;
    
97. lighting phong
    
98. set(gca,'GridLineStyle',':')
    
99. box on
    
100. ax=gca;
     
101. ax.BoxStyle = 'full';
     
102. grid on
     
103. 
     
104. 
     
105. 
     
106. %% logs
     
107. % typed by : mm
     
108. % mod : 2017年 03月 16日 星期四 15:23:22 CST
     
109. % contact me : meatball1982@163.com
     
110. %

复制代码

---

欢迎光临 Mathematica中文论坛-非官方 (http://ilovemathematica.com/)

Powered by Discuz! X3.2