matlab 散点边界，外表面, boundary, alphashape 儿子的papa

meatball1982

用matlab 找散点的表面。3D或是2D.有boundary函数和 alphaShape函数。

1. clear all
   
2. clc
   
3. clf
   
4. 
   
5. %% outline
   
6. % find the out surface of vdw points
   
7. 
   
8. %% main
   
9. % % % load position of LYS
   
10. % % pos=load('../data/LYS/LYS_xyz');
    
11. % %
    
12. % % n_at=length(pos(:,1));
    
13. % %
    
14. % % tm=load('../data/LYS/LYS_vdw.txt');
    
15. % % r_mat=tm(:,3)/2;
    
16. % %
    
17. % % save Mat_post_pos_rand.mat pos r_mat
    
18. 
    
19. %% load data
    
20. load Mat_post_pos_rand.mat
    
21. n_at=length(pos(:,1));
    
22. 
    
23. % sphere surface ---------------------------------
    
24. [sp_x,sp_y,sp_z]=sphere(30);
    
25. [n_sp_x,n_sp_y]=size(sp_x);
    
26. 
    
27. % for store points pos
    
28. x_mat=[];
    
29. y_mat=[];
    
30. z_mat=[];
    
31. 
    
32. for i=1:n_at
    
33.     % generate sphere surface points -------------
    
34.     at_x=sp_x*r_mat(i)+pos(i,1);
    
35.     at_y=sp_y*r_mat(i)+pos(i,2);
    
36.     at_z=sp_z*r_mat(i)+pos(i,3);
    
37.    
    
38.     x_mat=[x_mat;at_x(:)];
    
39.     y_mat=[y_mat;at_y(:)];
    
40.     z_mat=[z_mat;at_z(:)];
    
41. end
    
42. 
    
43. 
    
44. 
    
45. 
    
46. % % boundary version ------------------------------------
    
47. % % boundary 2D------------------
    
48. % n_layer=61;
    
49. % z_lin=linspace(min(unique(z_mat)),max(unique(z_mat)),n_layer);
    
50. %
    
51. % for i=1:(n_layer-1)
    
52. %     ind=z_lin(i)<=z_mat & z_mat< z_lin(i+1);
    
53. %     x_lay=x_mat(ind);
    
54. %     y_lay=y_mat(ind);
    
55. %     z_lay=z_mat(ind);
    
56. %     clf
    
57. %     h=gcf;
    
58. %     set(h,'position',[1100,100,400, 300]);
    
59. %     plot(x_lay,y_lay,'.')
    
60. %     axis equal
    
61. %     k=boundary(x_lay,y_lay,0.8);
    
62. %     hold on
    
63. %     plot(x_lay(k),y_lay(k),'r-');
    
64. %     axis equal
    
65. %     axis tight
    
66. %     axis([-5 7 -6 7])
    
67. % %     h=gcf;
    
68. % %     fig_na=['../imgs/LYS/VDW/boundary/fig_points_2D_',mat2str(i)];
    
69. % %     fun_work_li_035_myfig_out(h,fig_na,3);
    
70. % end
    
71. 
    
72. % % boundary 3D ------------------
    
73. % k=boundary([x_mat';y_mat';z_mat']');
    
74. % h=gcf;
    
75. % set(h,'position',[1100,100,800, 600]);
    
76. % plot3(x_mat(k),y_mat(k),z_mat(k),'r.')
    
77. % axis equal
    
78. % axis tight
    
79. % view(-20,30)
    
80. % % h=gcf;
    
81. % % fig_na=['../imgs/LYS/VDW/boundary/fig_boundary_3D'];
    
82. % % fun_work_li_035_myfig_out(h,fig_na,3);
    
83. 
    
84. 
    
85. % % alphaShape version -----------------------------------
    
86. % shp=alphaShape(x_mat,y_mat,z_mat,0.55);
    
87. % % boundary surface points ind
    
88. % ind=unique(shp.boundaryFacets);
    
89. % shp_P=shp.Points;
    
90. %
    
91. % h=gcf;
    
92. % set(h,'position',[1100,100,800, 600]);
    
93. %
    
94. % plot3(x_mat,y_mat,z_mat,'b.')
    
95. % % plot3(shp_P(ind,1),shp_P(ind,2),shp_P(ind,3),'r.')
    
96. % plot(shp,'edgecolor','none','facecolor','r','facealpha',0.4)
    
97. % axis equal
    
98. % axis tight
    
99. % view(-20,30)
    
100. % % h=gcf;
     
101. % % fig_na=['../imgs/LYS/VDW/fig_alphashape_3D'];
     
102. % % fun_work_li_035_myfig_out(h,fig_na,3);
     
103. %
     
104. % plot3(x_mat,y_mat,z_mat,'b.')
     
105. % axis equal
     
106. % axis tight
     
107. % view(-20,30)
     
108. % % h=gcf;
     
109. % % fig_na=['../imgs/LYS/VDW/fig_points_3D'];
     
110. % % fun_work_li_035_myfig_out(h,fig_na,3);
     
111. %
     
112. %
     
113. % plot3(shp_P(ind,1),shp_P(ind,2),shp_P(ind,3),'r.')
     
114. % axis equal
     
115. % axis tight
     
116. % view(-20,30)
     
117. % % h=gcf;
     
118. % % fig_na=['../imgs/LYS/VDW/fig_alphashape_3D_poi'];
     
119. % % fun_work_li_035_myfig_out(h,fig_na,3);
     
120. 
     
121. %% logs
     

复制代码