请教：三重循环优化
问题：判断两凸多面体是否相交方法：检测第一个多面体的每个面是否与第二个多面体的面有交点(r0=(abs(z1-z2)<=.5);)
分析：多面体每个面为三角形，三顶点M1,M2,M3确定一平面。该平面的法线向量n与向量M1M2,M1M3都垂直，
所以n=M1M2叉乘M1M3，提取系数a,b,c.根据平面的点法式方程，得所求平面的方程：
a(x-x1)+b(y-y1)+c(z-z1)=0;
效率极其低，大家有没有什么好的方法，不妨分享以下！
小弟先谢谢了！
程序：function N=myoverlaptest
% 2008 4 5 4:34
n=50;
X1=rand(n,3).*4;
c1=convhulln(X1);
X2=X1+3.5;
c2=convhulln(X2);
number=2;
len1=size(c1,1)
figure,hold on,view(3)
for i=1:size(c1,1)
j=c1(i,);
patch(X1(j,1),X1(j,2),X1(j,3),[.5 .5 .5])
end
for i=1:size(c2,1)
j=c2(i,);
patch(X2(j,1),X2(j,2),X2(j,3),[.5 .5 .5]);
end
N=0; % 不相交
m=size(c1,1);
mm=0;
for i=1:m
A=X1(c1(i,:),1)'; % x坐标 (1*3)
B=X1(c1(i,:),2)'; % y坐标
C=X1(c1(i,:),3)'; % z坐标
x=min(A):.1:max(A);
y=min(B):.1:max(B);
=meshgrid(x,y);
M1M2=;
M1M3=;
syms I J K
n1=det(); % 求平面的法线向量
a=double(maple('coeff',n1,I));
b=double(maple('coeff',n1,J));
c=double(maple('coeff',n1,K));
d=-(a*A(1)+b*B(1)+c*C(1));
z1=-(a/c*x+b/c*y+d/c);
for j=1:number-1
cj=c2;
X_convex1=X2;
for k=1:size(c2,1)
mm=mm+1;
A=X_convex1(cj(k,:),1)'
B=X_convex1(cj(k,:),2)'
C=X_convex1(cj(k,:),3)'
M1M2=
M1M3=
n2=det()
a=double(maple('coeff',n2,I))
b=double(maple('coeff',n2,J))
c=double(maple('coeff',n2,K))
d=-(a*A(1)+b*B(1)+c*C(1))
z2=-(a/c*x+b/c*y+d/c)
r0=(abs(z1-z2)<=.5);
zz=r0.*z1;yy=r0.*y;xx=r0.*x;
x2=;y2=
in=inpolygon(xx(r0~=0),yy(r0~=0),x2',y2');
if all(in==0)
continue % 不相交,则进入下轮循环
else
N=1
return
end
if mm>len1^2
disp('死循环');
return
end
end
end
end[
[ 本帖最后由 chengweifeng 于 2008-5-22 18:01 编辑 ] 两个多面体相交，重合部分的三维坐标肯定相等，先将x轴坐标相等的数据a选出来，再将a中y轴坐标相等的数据b选出来，最后将b中z轴坐标相等的数据c选出来，c就是两个多面体相交的数据。 Suppose there are two convex bodies A and B whose vertices are A1, A2,...,An and B1,B2,...Bm, respectively.I think the following statements are equivalent:
(1) A \cap B = \emptyset
(2) for all i=1,2,...,n,Ai \notin B
(3) for all i=1,2,...,m,Bi \notin A
So instead of checking the intersection of faces between A and B, I will proceed with checking whether all vertices of A are out of B.If it is true, A and B have no intersection. to:hhwuai007
你的方法让我有种焕然一新的感觉，但是个人暂时的感觉 ，你的方法在程序上实现恐怕更为不易。因为“重合部分的三维坐标肯定相等”----〉最终还是归结到求彼此相交面的交点，好像又回到了我的方法上来。那么，效率问题仍然是没有跨越。 考虑体B的所有线段是否与体A面相交，同理考虑A
得到体A某面方程F(X,Y,Z)=0
带入体B某线段两坐标
G=F(点一）×F（点二）
是否存在
G<0线面相交
G>0不相交
G=0线有点在面上
没仔细考虑，直接根据直线和直线判断相交类比过来
[ 本帖最后由 shunfly 于 2008-4-5 13:55 编辑 ] To:pcqsl
Actually,checking the intersection of two convex hulln,two conditions must be considered.The following pictures show whatshould be taken into account:
I think the ways you supplied are just available for the second conditon.Maybe the following function show what your mean,I call it Volume Add text Method
function N=AddVolume
% 体积和法测试两多面体是否相交，相交则N=1,反之N=0
% 但是此法不能判断另一类特殊的多面体相交情况
% 2008 3 23cheng wf
% chengwf@yahoo.cn
clear all
clc
=sphere(20);
i=1+round(rand.*4);
x1=x(11,i);y1=y(11,i);z1=z(11,i);
x2=x(11,i+5);y2=y(11,i+5);z2=z(11,i+5);
x3=x(11,i+10);y3=y(11,i+10);z3=z(11,i+10);
x4=x(11,i+15);y4=y(11,i+15);z4=z(11,i+15);
i1_row=1;i1_coloun=1+round(rand.*19);
j1_row=21;j1_coloun=1+round(rand.*19);
z5=z(i1_row,i1_coloun);
x5=x(i1_row,i1_coloun);
y5=y(i1_row,i1_coloun);
z6=z(j1_row,j1_coloun);
x6=x(j1_row,j1_coloun);
y6=y(j1_row,j1_coloun);
if i<=3
j=i+2;
elseif i>3&i<=5
j=i-2;
end
x7=x(16,j);y7=y(16,j);z7=z(16,j);
x8=x(16,j+5);y8=y(16,j+5);z8=z(16,j+5);
x9=x(16,j+10);y9=y(16,j+10);z9=z(16,j+10);
x10=x(16,j+15);y10=y(16,j+15);z10=z(16,j+15);
x11=x(6,j);y11=y(6,j);z11=z(6,j);
x12=x(6,j+5);y12=y(6,j+5);z12=z(6,j+5);
x13=x(6,j+10);y13=y(6,j+10);z13=z(6,j+10);
x14=x(6,j+15);y14=y(6,j+15);z14=z(6,j+15);
xi=;
yi=;
zi=;
X=;
=convhulln(X);
hold on
for i=1:size(c,1)
j=c(i,);
patch(X(j,1),X(j,2),X(j,3),[.5 .5 .5]);
end
view(3)
% X2=; % 不相交 N=0
X2=; % 相交 N=1
=convhulln(X2);
for i=1:size(c2,1);
j=c2(i,);
patch(X2(j,1),X2(j,2),X2(j,3),[.5 .3 .2]);
end
% ======相交判断=====
% 体积和测试法
i=0;
j=0;
% N=0;% 不相交
for i=1:length(X) % 对第1个多面体的每个顶点进行检测
v=0;
for j=1:size(c2,1)
k=c2(j,;
x_convexj=;
y_convexj=;
z_convexj=;
convexj=;
=convhulln(convexj);
v=v+vj;
end
if abs(v-V2)<5*10e-5
% if v==V2
N=1.0% 相交
return
else
N=0 % 是否相交，待定
% 调用其他函数判断另一种相交情况
end
end
return
[ 本帖最后由 chengweifeng 于 2008-4-5 21:24 编辑 ] To: shunfly 兄
你的方法理论上可行，但是就我#1的帖子看，其实很难提取出每个多面体的棱.另外，此法同样难逃三重循环的限制。你试一下我#1的帖子就知道效率是何等的低。:'( 原帖由 chengweifeng 于 2008-4-5 14:18 发表 http://www.simwe.com/forum/images/common/back.gif
To:pcqsl
Actually,checking the intersection of two convex hulln,two conditions must be considered.The following pictures show whatshould be taken into account:
You are correct.I made mistakes.
By giving more thought to your problem, I would like to proceed with Matlab function fmincon(). The fundamental reason is that *ANY POINT* inside a *CONVEX HULL* can be expressed with a linear function of the vertices. Let us define convex bodies A and B with vertices A1...An, and B1...Bm.Then
for all a \in A, exists p1,p2,...,pn where for i=1...n, 0<= pi<=1, and p1+p2+...+pn=1,such that a=A1*p1+A2*p2+...+An*pn.With these constrains on p1,...pn, the express never give points outside of A.There is a similar result on B.
Then your problem is equivalent to find the minimum distance between two points, where one is from A and the other is from B. Compare to your algorithm, this program should have betterperformance.However, since it is based on the "Medium-Scale Algorithm", it will be slow if the number of vertices is large.%%%%%main program
n=50;
X1=rand(n,3).*4;
c1=convhulln(X1);
X2=X1+3.5;
c2=convhulln(X2);
a=X1(unique(c1),:); %only vertices involved
b=X2(unique(c2),:);
N=checkintersection(a,b)%N=1 means there exists an intersection
%%%%%sub program
function r=checkintersection(a,b)
a=a';
b=b';
Na=size(a,2);
Nb=size(b,2);
dim=size(a,1);
Acon=zeros(1,Na+Nb);
bcon=1;
Aeq=];
beq=;
dist=@(p)(distance(p,a,b));
p0=;
options= optimset('display','off','LargeScale','off','outputfcn',@outfun);
=fmincon(dist,p0,Acon,bcon,Aeq,beq,zeros(Na+Nb,1),ones(Na+Nb,1),[],options);
if exitflag>0 || exitflag==-1
r=true;
else
r=false;
end;
%%%%%calculate the distance,points obtained from coefficients x and vertices a and b
function d=distance(x,a,b)
v1=a*x(1:size(a,2));
v2=b*x(size(a,2)+1:end);
d=norm(v1-v2);
%%%%%To check whether a small distance presents, if yes, stop.
%%%%%This is a trick for the case where the intersection is "easy" to find.
function stop=outfun(x, optimValues, state)
stop=false;
delta=1e-16;
if optimValues.fval<delta % distance is small enough
if sum( x>-delta )+ sum(x<1+delta) == 2*length(x)
% coefficients meet constraints ( in other words, the point is inside the convex bodies.)
stop=true;
end
end;
end;
[ 本帖最后由 pcqsl 于 2008-4-6 11:15 编辑 ]
for all a \in A, exists p1,p2,...,pn where for i=1...n, 0<= pi<=1, and p1+p2+...+pn=1,such that a=A1*p1+A2*p2+...+An*pn.With these constrains on p1,...pn, the express never give points outside of A.
pn is a weighted coefficient,right?!You are Great!
Out of regard for your thought,I will take it into account seriously.
The result will be posted here after my test on your function
Thanks again!
[ 本帖最后由 chengweifeng 于 2008-4-6 16:22 编辑 ] To:pcqsl
*ANY POINT* inside a *CONVEX HULL* can be expressed with a linear function of the vertices.
The code you supplied is more efficient than mine.Thank you very much!
but it still works only for the second condition!
I think the method supplied by hhwuai007 may be a good way! 原帖由 chengweifeng 于 2008-4-26 21:47 发表 http://www.simwe.com/forum/images/common/back.gif
but it still works only for the second condition!
I don't think you understand this algorithm fully.My code apparently works for both cases.I give an example for the first situation:a=
b=
figure;
hold off
plot3(a(:,1),a(:,2),a(:,3),'r.','markersize',20)
hold on;
plot3(b(:,1),b(:,2),b(:,3),'m.','markersize',20)
grid on
ha=convhulln(a)
for i=1:size(ha,1)
j=ha(i,);
patch(a(j,1),a(j,2),a(j,3),[.5 .5 .5])
end
hold on
hb=convhulln(b)
for i=1:size(hb,1)
j=hb(i,);
patch(b(j,1),b(j,2),b(j,3),[.5 .5 .5]*1.5)
end
alpha(.6)
N=checkintersection(a,b)%N=1 means there exists an intersectionN returns 1
Can you show a counter example when my code failed?
[ 本帖最后由 pcqsl 于 2008-4-27 01:11 编辑 ] To pcqsl:
是的，我错了，程序运行良好！
感谢pcqsl的鼎力相助，帮了我的大忙！:victory:
同时也向该论坛所有帮助和鼓励过我的兄弟们表示感谢，因为你们，这段孤独的旅程才会走的那么坚决！ 惩罚函数法：有约束的最优化问题 ——〉无约束最优化问题
惩罚函数法的基本思想：借助惩罚函数把约束问题转化为无约问题，进而用无约束最优化方法来求解。由于约束的非线性，不能用消元法将问题转化为无约束问题，因此，在求解时必须同时兼顾到既使目标函数值下降，又要满足约束条件这两个方面。实现这一点的途径是由目标函数和约束函数组成辅助函数，把原来的约束问题转化为极小化辅助问题的无约束问题。
惩罚函数法包括：内部惩罚函数法(内点法)、外部惩罚函数法（外点法）、乘子法
内点法：对企图从内部穿越可行域边界的点在目标函数中加入相应的“障碍”，距边界越近，障碍越大，在边界上给以无穷大的障碍，从而保障迭代一直在可行域内部移动。
外点法：对违反约束的点在目标函数中加入相应的“惩罚”，而对可行点不予惩罚。此法的迭代点一般在可行域外部移动。
乘子法：？
***：在通常遇到的实际问题中，如果根据问题的性质，知道函数f(x,y)的最小值（最大值）一定在D的内部取得，而函数在D内 只有一个驻点，那么可以肯定该驻点处的函数值就是函数f(x,y)在D上的最小/最大值。
***：极值问题，对于函数的自变量，除了限制在函数的定义域内以外，并无其他条件-------〉无约束优化
在实际问题中，有时会遇到对函数的自变量还有附加条件的极值问题-------〉有约束优化问题 最近重新关注这个问题，并对pcqsl提出的方法非常感兴趣，只是对8#程序中的
“Acon=zeros(1,Na+Nb); bcon=1;”两句始终理解不透，不知道其理论意义及其在
fmincon函数中的作用/意义何在？恳请路过高手指点迷津。不胜感激！ 14# chengweifeng 呵呵，还能翻出三年之前的帖子，不容易啊，
页:
[1]