matlab画任意多面体,[转载]matlab判断一个点是否在多面体内

要在空位区域随机放置一定数量的原子，这些原子在空位区域任何一处存在的概念是相同的。空位区域是由包围这个空位周边的一些原子定义的。

如果这个空位区域是一个标准的长方体，那么问题就比较简单，只需要产生随机数，然后再将随机数沿着基矢方向进行相应的缩放。

对于不规则的空间区域，也可以采用类似的思想：将空位区域(多面体)扩大到一个长方体，即长方体刚好是多面体的包络。然后在长方体内部随机产生点，如果点在多面体内部就保留；不在多面体内部就舍去，重新产生。

这其中就出现一个基础问题：如何判断一个点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点在四面体外部。

具体代码如下：

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

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