|  | 
沙发
 
 
 楼主|
发表于 2018-1-11 12:50:24
|
只看该作者 
| 更新: 3D版本的
 
 
     
     
 主程序
 
 复制代码clear all
clc
clf
%% outline
% gen rand circle, sphere in 3D 
% according : http://www.ilovematlab.cn/thread-526875-1-1.html
% gen cen
% gen rad with cen
% 
%% main
% % parameters
% n =200;
% b =10;
% rmin=0.1;
% rmax=0.3;
% 
% % gen initial center
% [ cen ] = fun_gen_cen_3d(n,b);
% 
% % distance
% D=pdist(cen);
% D_mat = squareform(D);
% 
% 
% % generate r ind for cal ri+rj distance 
% x1 = [1:n];              % ind
% [ind1,ind2]=ndgrid(x1,x1);
% out_ind=[ind2(:),ind1(:)];
% ind_cho = out_ind(:,1)<out_ind(:,2);
% out_ind(~ind_cho,:)=[];
% 
% % gen initial r 
% r = fun_gen_r(n,D_mat,rmin,rmax);
% R_mat=squareform(sum(r(out_ind)')');
% 
% % judge
% dis=(D_mat-R_mat)>0;
% n_flg = n^2-n;
% 
% n_cont = 1;
% while sum(dis(:))~= n_flg       | ... % dis
%         sum(cen(:,1)-r < 0 )> 0 | ... % left bound
%         sum(cen(:,2)-r < 0 )> 0 | ... % bottom bound
%         sum(cen(:,3)-r < 0 )> 0 | ... % z dir
%         sum(cen(:,1)+r > b )> 0 | ... % right bound
%         sum(cen(:,2)+r > b )> 0 | ... % up bound
%         sum(cen(:,3)+r > b )> 0
%         
%     n_cont=n_cont+1
%     [ cen ] = fun_gen_cen_3d(n,b);              % new centers
%     r  = fun_gen_r(n,D_mat,rmin,rmax);           % new r
%     R_mat=squareform(sum(r(out_ind)')');  % ri+rj matrix
%     D = pdist(cen);
%     D_mat = squareform(D);          % center matrix
%     dis=(D_mat-R_mat)>0;
% end
% 
% 
% save Mat_4plot.mat
load Mat_4plot.mat
%% output result
hold on
plot3(cen(:,1),cen(:,2),cen(:,3),'k.','markersize',3)
axis([0 b 0 b 0 b])
axis square
ax = gca;
grid on
ax.XTick = [0:1:10];
ax.YTick = [0:1:10];
ax.ZTick = [0:1:10];
% the=linspace(0,2*pi,100);
% the=[the,0];
% 
for k = 1: n
    [x_sp,y_sp,z_sp]=sphere(20);
    x_sp = r(k)*x_sp + cen(k,1);
    y_sp = r(k)*y_sp + cen(k,2);
    z_sp = r(k)*z_sp + cen(k,3);
    surf(x_sp,y_sp,z_sp,r(k)*ones(size(x_sp)),'edgecolor','none')
%     x_r = r(k)*x_sp+cen(k,1);
%     y_r = r(k)*sin(the)+cen(k,2);
%     plot(x_r,y_r,'-','linewidth',2)
end
box on
axis equal
colormap(jet)
alpha(0.5)
xlabel('x')
ylabel('y')
zlabel('z')
 3D的中心程序
 
 复制代码function [ cen ] = fun_gen_cen_3d(n,b)
tm=randperm(b*b*b);
ind  = sort(tm(1:n));
tm_0 = zeros(b,b,b);
tm_0(tm(1:n))=1;
[ind_x,ind_y,ind_z]=findND(tm_0);
% cen_x = ceil(ind/b)+0.5-rand(1,n)*0.35;
% cen_y = mod(ind,b) +0.5-rand(1,n)*0.35;
cen_x = ind_x' -0.5+rand(1,n)*0.35;
cen_y = ind_y' -0.5+rand(1,n)*0.35;
cen_z = ind_z' -0.5+rand(1,n)*0.35;
cen=[cen_x;
     cen_y;
     cen_z;]';
 
 cen = abs(cen);
 
 | 
 |