我遇到的一个实际问题是:要在空位区域随机放置一定数量的原子,这些原子在空位区域任何一处存在的概念是相同的。空位区域是由包围这个空位周边的一些原子定义的。
如果这个空位区域是一个标准的长方体,那么问题就比较简单,只需要产生随机数,然后再将随机数沿着基矢方向进行相应的缩放。
对于不规则的空间区域,也可以采用类似的思想:将空位区域(多面体)扩大到一个长方体,即长方体刚好是多面体的包络。然后在长方体内部随机产生点,如果点在多面体内部就保留;不在多面体内部就舍去,重新产生。
这其中就出现一个基础问题:如何判断一个点P是否在一个多面体内?多面体由空间三维点定义。
我的初步想法:1. 将多面体划分成四面体
2. 判断这个点是否在这些四面体内部
3. 如果这个点在任何一个四面体内部,那么这个点就在这个多面体内部;
4. 如果这个点不在任何一个四面体内部,那么这个点就不在这个多面体内部。
那么核心问题就转换为:1. 如何将多面体划分成许多的四面体?
2. 如何判断一个点是否在四面体内部?
对于第一个问题,matlab提供直接使用的函数 DelaunayTri 可以实现。
对于第二个问题,matlab也提供一个函数 tsearchn ,但是这个函数的健壮性比较差,至少我的问题没法解决。
没办法,在网上找到了有关的算法,自己写了代码。算法如下:
四面体由四个顶点定义:V1 = (x1, y1, z1)
V2 = (x2, y2, z2)
V3 = (x3, y3, z3)
V4 = (x4, y4, z4)
要检测的点定义为:P = (x, y, z)
可以根据这些点组成的4个行列式来进行判断:D0
|x1 y1 z1 1|
|x2 y2 z2 1|
|x3 y3 z3 1|
|x4 y4 z4 1|D1
|x y z 1|
|x2 y2 z2 1|
|x3 y3 z3 1|
|x4 y4 z4 1|D2
|x1 y1 z1 1|
|x y z 1|
|x3 y3 z3 1|
|x4 y4 z4 1|D3
|x1 y1 z1 1|
|x2 y2 z2 1|
|x y z 1|
|x4 y4 z4 1|D4
|x1 y1 z1 1|
|x2 y2 z2 1|
|x3 y3 z3 1|
|x y z 1|判据:
如果Di (i=1,2,3,4)与D0的符号都相同,即同为正,或同为负,那么P点就在四面体内部;
否则P点在四面体外部。function inflag = inpolyhedron(point_set,p_detected)
% point_set: a set of points stores the coordinates
% p_detected: point to be detected
% inflag:
% flag = 1: the point is in the polyhedron.
% flag = 0: the point is not in the polyhedron.
% Powered by: Xianbao Duan xianbao.d@gmail.com
% stores the coordinates of the convexes.
tri = DelaunayTri(point_set);
% number of the tetrahedrons decomposed from the polyhedron
num_tet = size(tri,1);
t_inflag = zeros(1,11);
for i = 1:num_tet
v1_coord = point_set(tri(i,1),:);
v2_coord = point_set(tri(i,2),:);
v3_coord = point_set(tri(i,3),:);
v4_coord = point_set(tri(i,4),:);
D0 =det( [v1_coord,1;v2_coord,1;v3_coord,1;v4_coord,1]);
D1 = det([p_detected,1;v2_coord,1;v3_coord,1;v4_coord,1]);
D2 = det([v1_coord,1;p_detected,1;v3_coord,1;v4_coord,1]);
D3 = det([v1_coord,1;v2_coord,1;p_detected,1;v4_coord,1]);
D4 = det([v1_coord,1;v2_coord,1;v3_coord,1;p_detected,1]);
if D0*D1 > 0 && D0*D2>0 && D0*D3>0 && D0*D4 > 0
t_inflag(i) = 1;
break;
end
end
if sum(t_inflag) > 0
inflag = 1;
% disp('The point is in the polyhedron.');
else
inflag = 0;
% disp('The point is not in the polyhedron.');
end