|
本帖最后由 meatball1982 于 2016-6-14 21:06 编辑
用matlab 找散点的表面。3D或是2D.有boundary函数和 alphaShape函数。
- clear all
- clc
- clf
- %% outline
- % find the out surface of vdw points
- %% main
- % % % load position of LYS
- % % pos=load('../data/LYS/LYS_xyz');
- % %
- % % n_at=length(pos(:,1));
- % %
- % % tm=load('../data/LYS/LYS_vdw.txt');
- % % r_mat=tm(:,3)/2;
- % %
- % % save Mat_post_pos_rand.mat pos r_mat
- %% load data
- load Mat_post_pos_rand.mat
- n_at=length(pos(:,1));
- % sphere surface ---------------------------------
- [sp_x,sp_y,sp_z]=sphere(30);
- [n_sp_x,n_sp_y]=size(sp_x);
- % for store points pos
- x_mat=[];
- y_mat=[];
- z_mat=[];
- for i=1:n_at
- % generate sphere surface points -------------
- at_x=sp_x*r_mat(i)+pos(i,1);
- at_y=sp_y*r_mat(i)+pos(i,2);
- at_z=sp_z*r_mat(i)+pos(i,3);
-
- x_mat=[x_mat;at_x(:)];
- y_mat=[y_mat;at_y(:)];
- z_mat=[z_mat;at_z(:)];
- end
- % % boundary version ------------------------------------
- % % boundary 2D------------------
- % n_layer=61;
- % z_lin=linspace(min(unique(z_mat)),max(unique(z_mat)),n_layer);
- %
- % for i=1:(n_layer-1)
- % ind=z_lin(i)<=z_mat & z_mat< z_lin(i+1);
- % x_lay=x_mat(ind);
- % y_lay=y_mat(ind);
- % z_lay=z_mat(ind);
- % clf
- % h=gcf;
- % set(h,'position',[1100,100,400, 300]);
- % plot(x_lay,y_lay,'.')
- % axis equal
- % k=boundary(x_lay,y_lay,0.8);
- % hold on
- % plot(x_lay(k),y_lay(k),'r-');
- % axis equal
- % axis tight
- % axis([-5 7 -6 7])
- % % h=gcf;
- % % fig_na=['../imgs/LYS/VDW/boundary/fig_points_2D_',mat2str(i)];
- % % fun_work_li_035_myfig_out(h,fig_na,3);
- % end
- % % boundary 3D ------------------
- % k=boundary([x_mat';y_mat';z_mat']');
- % h=gcf;
- % set(h,'position',[1100,100,800, 600]);
- % plot3(x_mat(k),y_mat(k),z_mat(k),'r.')
- % axis equal
- % axis tight
- % view(-20,30)
- % % h=gcf;
- % % fig_na=['../imgs/LYS/VDW/boundary/fig_boundary_3D'];
- % % fun_work_li_035_myfig_out(h,fig_na,3);
- % % alphaShape version -----------------------------------
- % shp=alphaShape(x_mat,y_mat,z_mat,0.55);
- % % boundary surface points ind
- % ind=unique(shp.boundaryFacets);
- % shp_P=shp.Points;
- %
- % h=gcf;
- % set(h,'position',[1100,100,800, 600]);
- %
- % plot3(x_mat,y_mat,z_mat,'b.')
- % % plot3(shp_P(ind,1),shp_P(ind,2),shp_P(ind,3),'r.')
- % plot(shp,'edgecolor','none','facecolor','r','facealpha',0.4)
- % axis equal
- % axis tight
- % view(-20,30)
- % % h=gcf;
- % % fig_na=['../imgs/LYS/VDW/fig_alphashape_3D'];
- % % fun_work_li_035_myfig_out(h,fig_na,3);
- %
- % plot3(x_mat,y_mat,z_mat,'b.')
- % axis equal
- % axis tight
- % view(-20,30)
- % % h=gcf;
- % % fig_na=['../imgs/LYS/VDW/fig_points_3D'];
- % % fun_work_li_035_myfig_out(h,fig_na,3);
- %
- %
- % plot3(shp_P(ind,1),shp_P(ind,2),shp_P(ind,3),'r.')
- % axis equal
- % axis tight
- % view(-20,30)
- % % h=gcf;
- % % fig_na=['../imgs/LYS/VDW/fig_alphashape_3D_poi'];
- % % fun_work_li_035_myfig_out(h,fig_na,3);
- %% logs
复制代码
|
|