|
|
马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。
您需要 登录 才可以下载或查看，没有账号？我要加入
x
单步走，不报错，可是一直接运行就死循环，恳请高手帮忙看看！
可能是两个子函数有问题，他们的终止条件是不是有问题?
function myconvexhull
clear all
clc
hold on
L_x=20;L_y=20;a_min=1;a_max=2;
S_min=4;S_max=10;
Area_all=0; % 生成的所有凸多边形的总面积
MYAREA=[]; % 存储每个凸多边形的面积
n=0;
x0=[];y0=[];
theta0=[];
tic
warning off MATLAB:divideByZero
while Area_all<=0.2.*L_x.*L_y
n=n+1;
x=[];y=[];R=[];
U1=rand;
S(n)=fix(U1*(S_max-S_min)+S_min);
U2=rand;
x0(n)=U2*(L_x-a_max)+1;
U3=rand;
y0(n)=U3*(L_y-a_max)+1;
U4=rand;
theta0(n)=U4*2*pi;
for i=1:S(n)
U5=rand;
R(i)=U5*(a_max/2-a_min/2)+a_min/2; %vector radius
x(i)=x0(n)+R(i).*cos(theta0(n)+(i-1).*2.*pi/S(n));
y(i)=y0(n)+R(i).*sin(theta0(n)+(i-1).*2.*pi/S(n));
end
% 避免畸形(狭长)
C=[];
for i=1:length(x)
eval(['A_x' num2str(i) '=repmat(x(i),1,length(x)-1)']);
eval(['A_y' num2str(i) '=repmat(y(i),1,length(x)-1)']);
j=find(x~=x(i));
eval(['B_x' num2str(i) '=x(j)']);eval(['B_y' num2str(i) '=y(j)']);
a=(eval(['A_x' num2str(i)])-eval(['B_x' num2str(i)]));
b=(eval(['A_y' num2str(i)])-eval(['B_y' num2str(i)]));
c=sqrt(a.^2+b.^2);
c=max(c);
C=[C,c];
end
Mymaxdistance(n)=max(C(:)); %第n个多边形的最大直径
Mymindistance(n)=min(C(:)); %第n个多边形的最小直径
if Mymaxdistance(n)>2.*Mymindistance(n)
n=n-1;
continue
end
% 得到二维凸壳
[ii,AREA_n]=convhull(x,y);
x=x(ii);y=y(ii);
S(n)=length(x);
x=[x,x(1)];y=[y,y(1)];
MYAREA=[MYAREA,AREA_n]; % 前n个凸多边形的面积
B{1,n}=[x',y'];
R=R(ii);R=[R,R(1)];
%-----------------------------
% 重叠判断
if n==1
plot(x,y)
axis([0 20 0 20])
end
if n>1
N=myoverlapping_determination(B,x,y,n)
if N~=n
continue
end
end
%-------------
% 侵入判断
if n==1
plot(x,y)
axis([0 20 0 20])
end
if n>1
M=myoverlay_determination(B,x,y,n)
if M~=n
continue
end
end
%---------------
% 面积率判断
Area_all=sum(MYAREA);
if Area_all>0.2.*L_x.*L_y
break
end
% 防止死循环
if n>100
disp('Maybe there is something wrong about this programm.')
lasterr
lastwarn
return
end
end
toc
%--------------
function N=myoverlapping_determination(B,x,y,n)
for i=1:n-1
XYi1=[];xi1=[];yi1=[];
XYi1=B{1,i};
xi1=XYi1(:,1)';yi1=XYi1(:,2)';
in=inpolygon(xi1,yi1,x,y)
if all(in==0)
plot(x,y)
axis([0 20 0 20])
N=n;
return
else
n=n-1;
N=n;
return
end
end
%--------------
function M=myoverlay_determintation(B,x,y,n)
nn=n;
k_n=[];
for i=1:length(x)-1
k_n(i)=(y(i+1)-y(i))./(x(i+1)-x(i));
end
a1=k_n;b1=repmat(-1,1,length(a1));
c1=k_n.*x(1:end-1)-y(1:end-1);
for i=1:n-1
%XYi2=[];xi2=[];yi2=[];
XYi2=B{1,i};
xi2=XYi2(:,1)';yi2=XYi2(:,2)';
k_i=[];
for j=1:length(xi2)-1
k_i(j)=(yi2(j+1)-yi2(j))./(xi2(j+1)-xi2(j));
end
a2=k_i;b2=repmat(-1,1,length(a2));
c2=k_i.*xi2(1:end-1)-yi2(1:end-1);
for k=1:length(x)-1
for l=1:length(xi2)-1
A1=[a1(k),b1(k);a2(l),b2(l)];
B1=[c1(k);c2(l)];
myintersection=A1\B1;
if isempty(myintersection)==0
n=n-1;
M=n;
return
else
plot(x,y)
axis([0 20 0 20])
M=n;
return
end
end
end
end
[ 本帖最后由 eight 于 2007-5-26 12:32 编辑 ] |
|
发表于 2007-5-25 17:04
提升卡
变色卡
显身卡
发表于 2007-5-25 17:07
楼主