【分享】不相交的随机大小的椭圆

funi

此函数产生的椭圆通过调整长径比可以用来模拟纤维混凝土中的纤维，当然也可以模拟类似的其它材料。

1. function myellipse
   
2. % cheng wf 2007 6 11 19:51
   
3. close all
   
4. clear all
   
5. clc
   
6. %format compact
   
7. depth=60;
   
8. a=[3 2 1 0.5];b=[2 1 0.5 0.25];
   
9. number=[10 10 100 20];
   
10. x0=[];y0=[];xlat=[];ylon=[];
    
11. [x0new1,y0new1,xlat1,ylon1]=myellipse_one(depth,a(1),b(1),number(1),x0,y0,xlat,ylon)
    
12. x0=x0new1;y0=y0new1;xlat=xlat1;ylon=ylon1;
    
13. for i=2:length(a)
    
14.     [x0new2,y0new2,xlat2,ylon2]=myellipse_rest(depth,a(i),b(i),number(i),x0,y0,xlat,ylon)
    
15.     x0=x0new2;y0=y0new2;xlat=xlat2;ylon=ylon2;
    
16. end
    
17. %  绘制椭圆图形
    
18. [m,n]=size(xlat)
    
19.     plot(xlat,ylon)
    
20.     axis([0 60 0 60])

复制代码

1. function [x0new,y0new,xlat,ylon]=myellipse_one(depth,a,b,number,x0,y0,xlat,ylon)
   
2. xxlen=length(x0);yylen=length(y0);
   
3. x0(xxlen+1)=a+(depth-2*a)*rand;
   
4. y0(yylen+1)=a+(depth-2*a)*rand;
   
5. ecc=axes2ecc(a,b);
   
6. [lat1,lon1]=ellipse1(x0(xxlen+1),y0(yylen+1),[a,ecc],rand*180);
   
7. xlat=[xlat,lat1];ylon=[ylon,lon1];
   
8. for i=2:number
   
9.     while 1
   
10.     x0(xxlen+i)=a+(depth-2*a)*rand;
    
11.     y0(yylen+i)=a+(depth-2*a)*rand;
    
12.     [lati,loni]=ellipse1(x0(xxlen+i),y0(yylen+i),[a,ecc],rand*180);
    
13.     xlat=[xlat,lati];ylon=[ylon,loni];
    
14. 
    
15.     for j=1:xxlen+i-1
    
16.         in1=[];in2=[];
    
17.         in1=inpolygon(xlat(:,xxlen+i),ylon(:,yylen+i),xlat(:,j),ylon(:,j))
    
18.         in2=inpolygon(xlat(:,j),ylon(:,j),xlat(:,xxlen+i),ylon(:,yylen+i))
    
19.         if all(in1==0)&&all(in2==0)
    
20.             n_j=j;
    
21.             if n_j==xxlen+i-1
    
22.                N=1
    
23.                 break
    
24.             end
    
25.             continue
    
26.         else
    
27.               N=2;
    
28.               xlat(:,xxlen+i)=[];ylon(:,xxlen+i)=[];x0(xxlen+i)=[];y0(yylen+i)=[];
    
29.               break
    
30.          end
    
31.      end
    
32.       if N==1
    
33.         break
    
34.       end
    
35.    end
    
36. end
    
37. x0new=x0;y0new=y0;xlat=xlat;ylon=ylon;

复制代码

1. function [x0new,y0new,xlat,ylon]=myellipse_rest(depth,a,b,number,x0,y0,xlat,ylon)
   
2. xxlen=length(x0);yylen=length(y0);
   
3. ecc=axes2ecc(a,b);N=0;
   
4. for i=1:number
   
5.      while 1
   
6.     x0(xxlen+i)=a+(depth-2*a)*rand;
   
7.     y0(xxlen+i)=a+(depth-2*a)*rand;
   
8.     [lat,lon]=ellipse1(x0(xxlen+i),y0(xxlen+i),[a,ecc],rand*180);
   
9.    xlat=[xlat,lat];ylon=[ylon,lon];
   
10.     for j=1:xxlen+i-1
    
11.         in1=[];in2=[];
    
12.         in1=inpolygon(xlat(:,xxlen+i),ylon(:,xxlen+i),xlat(:,j),ylon(:,j));
    
13.         in2=inpolygon(xlat(:,j),ylon(:,j),xlat(:,xxlen+i),ylon(:,yylen+i));
    
14.         if all(in1==0)&&all(in2==0)
    
15.              n_j=j;
    
16.              if n_j==xxlen+i-1
    
17.                 N=1;
    
18.                 break
    
19.             end
    
20.            continue
    
21.        else
    
22.              N=2;
    
23.              xlat(:,xxlen+i)=[];ylon(:,xxlen+i)=[];x0(xxlen+i)=[];y0(yylen+i)=[];
    
24.               break
    
25.        end
    
26.     end
    
27.        if N==1
    
28.        break
    
29.        end
    
30.   end
    
31. end
    
32. x0new=x0;y0new=y0;xlat=xlat;ylon=ylon;

复制代码