请教：三重循环优化

chengweifeng · 发表于 2008-4-5 07:07:38

问题：判断两凸多面体是否相交
方法：检测第一个多面体的每个面是否与第二个多面体的面有交点(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,[1 2 3 1]);
patch(X1(j,1),X1(j,2),X1(j,3),[.5 .5 .5])
end
for i=1:size(c2,1)
j=c2(i,[1 2 3 1]);
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);
[x,y]=meshgrid(x,y);
M1M2=[A(2)-A(1),B(2)-B(1),C(2)-C(1)];
M1M3=[A(3)-A(1),B(3)-B(1),C(3)-C(1)];
syms I J K
n1=det([I,J,K;M1M2;M1M3]); % 求平面的法线向量
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=[A(2)-A(1),B(2)-B(1),C(2)-C(1)]
M1M3=[A(3)-A(1),B(3)-B(1),C(3)-C(1)]
n2=det([I,J,K;M1M2;M1M3])
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=[A,A(1)];y2=[B,B(1)]
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 编辑 ]

hhwuai007 · 发表于 2008-4-5 10:48:14

两个多面体相交，重合部分的三维坐标肯定相等，先将x轴坐标相等的数据a选出来，再将a中y轴坐标相等的数据b选出来，最后将b中z轴坐标相等的数据c选出来，c就是两个多面体相交的数据。

pcqsl · 发表于 2008-4-5 10:57:44

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.

chengweifeng · 发表于 2008-4-5 13:10:38

to:hhwuai007
你的方法让我有种焕然一新的感觉，但是个人暂时的感觉，你的方法在程序上实现恐怕更为不易。因为“重合部分的三维坐标肯定相等”----〉最终还是归结到求彼此相交面的交点，好像又回到了我的方法上来。那么，效率问题仍然是没有跨越。

shunfly · 发表于 2008-4-5 13:54:20

考虑体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 编辑 ]

chengweifeng · 发表于 2008-4-5 14:18:16

To:pcqsl
Actually,checking the intersection of two convex hulln,two conditions must be considered.The following pictures show what should 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 23 cheng wf
% chengwf@yahoo.cn
clear all
clc
[x,y,z]=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=[x1;x2;x3;x4;x5;x6;x7;x8;x9;x10;x11;x12;x13;x14];
yi=[y1;y2;y3;y4;y5;y6;y7;y8;y9;y10;y11;y12;y13;y14];
zi=[z1;z2;z3;z4;z5;z6;z7;z8;z9;z10;z11;z12;z13;z14];
X=[xi,yi,zi];
[c,V1]=convhulln(X);
hold on
for i=1:size(c,1)
j=c(i,[1 2 3 1]);
patch(X(j,1),X(j,2),X(j,3),[.5 .5 .5]);
end
view(3)
% X2=[xi+1.5,yi,zi]; % 不相交 N=0
X2=[xi+.5,yi,zi]; % 相交 N=1
[c2,V2]=convhulln(X2);
for i=1:size(c2,1);
j=c2(i,[1 2 3 1]);
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=[X2(k,1);X(i,1)];
y_convexj=[X2(k,2);X(i,2)];
z_convexj=[X2(k,3);X(i,3)];
convexj=[x_convexj,y_convexj,z_convexj];
[L,vj]=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 编辑 ]

chengweifeng · 发表于 2008-4-5 14:25:52

To: shunfly 兄
你的方法理论上可行，但是就我#1的帖子看，其实很难提取出每个多面体的棱.另外，此法同样难逃三重循环的限制。你试一下我#1的帖子就知道效率是何等的低。:'(

pcqsl · 发表于 2008-4-6 08:17:45

原帖由 chengweifeng 于 2008-4-5 14:18 发表
To:pcqsl
Actually,checking the intersection of two convex hulln,two conditions must be considered.The following pictures show what should 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 better  performance.  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=[ones(1,Na) zeros(1,Nb);zeros(1,Na) ones(1,Nb);[a -b]];
beq=[1;1;zeros(dim,1)];
dist=@(p)(distance(p,a,b));
p0=[1;zeros(Na-1,1);1;zeros(Nb-1,1)];
options= optimset('display','off','LargeScale','off','outputfcn',@outfun);
[p,fval,exitflag]=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 编辑 ]

chengweifeng · 发表于 2008-4-6 16:00:22

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 编辑 ]

chengweifeng · 发表于 2008-4-26 21:47:34

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!

pcqsl · 发表于 2008-4-27 01:06:19

原帖由 chengweifeng 于 2008-4-26 21:47 发表

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=[0 0 0;1 0 0;1 1 0;1 0 1]
b=[0.5 1 0.3;1 1.5 0.3;1 2 .5;.75 -1 .5]
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,[1 2 3 1]);
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,[1 2 3 1]);
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 intersection

复制代码

N returns 1

Can you show a counter example when my code failed?

[ 本帖最后由 pcqsl 于 2008-4-27 01:11 编辑 ]

chengweifeng · 发表于 2008-5-6 09:27:44

To pcqsl:
是的，我错了，程序运行良好！
感谢pcqsl的鼎力相助，帮了我的大忙！:victory:
同时也向该论坛所有帮助和鼓励过我的兄弟们表示感谢，因为你们，这段孤独的旅程才会走的那么坚决！

chengweifeng · 发表于 2008-5-22 00:15:35

惩罚函数法：有约束的最优化问题 ——〉无约束最优化问题

惩罚函数法的基本思想：借助惩罚函数把约束问题转化为无约问题，进而用无约束最优化方法来求解。由于约束的非线性，不能用消元法将问题转化为无约束问题，因此，在求解时必须同时兼顾到既使目标函数值下降，又要满足约束条件这两个方面。实现这一点的途径是由目标函数和约束函数组成辅助函数，把原来的约束问题转化为极小化辅助问题的无约束问题。

惩罚函数法包括：内部惩罚函数法(内点法)、外部惩罚函数法（外点法）、乘子法

内点法：对企图从内部穿越可行域边界的点在目标函数中加入相应的“障碍”，距边界越近，障碍越大，在边界上给以无穷大的障碍，从而保障迭代一直在可行域内部移动。

外点法：对违反约束的点在目标函数中加入相应的“惩罚”，而对可行点不予惩罚。此法的迭代点一般在可行域外部移动。

乘子法：？

***：在通常遇到的实际问题中，如果根据问题的性质，知道函数f(x,y)的最小值（最大值）一定在D的内部取得，而函数在D内只有一个驻点，那么可以肯定该驻点处的函数值就是函数f(x,y)在D上的最小/最大值。

***：极值问题，对于函数的自变量，除了限制在函数的定义域内以外，并无其他条件-------〉无约束优化
在实际问题中，有时会遇到对函数的自变量还有附加条件的极值问题-------〉有约束优化问题

chengweifeng · 发表于 2011-4-2 09:20:10

最近重新关注这个问题，并对pcqsl提出的方法非常感兴趣，只是对8#程序中的
“Acon=zeros(1,Na+Nb); bcon=1;”两句始终理解不透，不知道其理论意义及其在
fmincon函数中的作用/意义何在？恳请路过高手指点迷津。不胜感激！

liuyalong008 · 发表于 2011-4-2 10:01:42

14# chengweifeng 呵呵，还能翻出三年之前的帖子，不容易啊，

账号		自动登录	找回密码
密码			注册

账号		自动登录	找回密码
密码			立即注册

请教：三重循环优化

相关帖子

本帖子中包含更多资源

本帖子中包含更多资源