Coverage for trimesh/ray/ray_util.py: 94%

1import numpy as np 

2 

3from .. import bounds, constants, util 

4from ..typed import ArrayLike 

5 

6 

7@constants.log_time 

8def contains_points(intersector, points: ArrayLike, check_direction: bool | None = None): 

9 """ 

10 Check if a mesh contains a set of points, using ray tests. 

11 

12 If the point is on the surface of the mesh, behavior is 

13 undefined. 

14 

15 Parameters 

16 --------- 

17 mesh: Trimesh object 

18 points: (n,3) points in space 

19 

20 Returns 

21 --------- 

22 contains : (n) bool 

23 Whether point is inside mesh or not 

24 """ 

25 # convert points to float and make sure they are 3D 

26 points = np.asanyarray(points, dtype=np.float64) 

27 if not util.is_shape(points, (-1, 3)): 

28 raise ValueError("points must be (n,3)") 

29 

30 # placeholder result with no hits we'll fill in later 

31 contains = np.zeros(len(points), dtype=bool) 

32 

33 # cull points outside of the axis aligned bounding box 

34 # this avoids running ray tests unless points are close 

35 inside_aabb = bounds.contains(intersector.mesh.bounds, points) 

36 

37 # if everything is outside the AABB, exit early 

38 if not inside_aabb.any(): 

39 return contains 

40 

41 # default ray direction is random, but we are not generating 

42 # uniquely each time so the behavior of this function is easier to debug 

43 default_direction = np.array([0.4395064455, 0.617598629942, 0.652231566745]) 

44 if check_direction is None: 

45 # if no check direction is specified use the default 

46 # stack it only for points inside the AABB 

47 ray_directions = np.tile(default_direction, (inside_aabb.sum(), 1)) 

48 else: 

49 # if a direction is passed use it 

50 ray_directions = np.tile( 

51 np.array(check_direction).reshape(3), (inside_aabb.sum(), 1) 

52 ) 

53 

54 # cast a ray both forwards and backwards 

55 # need multiple_hits=True so the parity test counts cavity walls 

56 _location, index_ray, _c = intersector.intersects_location( 

57 np.vstack((points[inside_aabb], points[inside_aabb])), 

58 np.vstack((ray_directions, -ray_directions)), 

59 multiple_hits=True, 

60 ) 

61 

62 # if we hit nothing in either direction just return with no hits 

63 if len(index_ray) == 0: 

64 return contains 

65 

66 # reshape so bi_hits[0] is the result in the forward direction and 

67 # bi_hits[1] is the result in the backwards directions 

68 bi_hits = np.bincount(index_ray, minlength=len(ray_directions) * 2).reshape((2, -1)) 

69 # a point is probably inside if it hits a surface an odd number of times 

70 bi_contains = np.mod(bi_hits, 2) == 1 

71 

72 # if the mod of the hit count is the same in both 

73 # directions, we can save that result and move on 

74 agree = np.equal(*bi_contains) 

75 

76 # in order to do an assignment we can only have one 

77 # level of boolean indexes, for example this doesn't work: 

78 # contains[inside_aabb][agree] = bi_contains[0][agree] 

79 # no error is thrown, but nothing gets assigned 

80 # to get around that, we create a single mask for assignment 

81 mask = inside_aabb.copy() 

82 mask[mask] = agree 

83 

84 # set contains flags for things inside the AABB and who have 

85 # ray tests that agree in both directions 

86 contains[mask] = bi_contains[0][agree] 

87 

88 # if one of the rays in either direction hit nothing 

89 # it is a very solid indicator we are in free space 

90 # as the edge cases we are working around tend to 

91 # add hits rather than miss hits 

92 one_freespace = (bi_hits == 0).any(axis=0) 

93 

94 # rays where they don't agree and one isn't in free space 

95 # are deemed to be broken 

96 broken = np.logical_and(np.logical_not(agree), np.logical_not(one_freespace)) 

97 

98 # if all rays agree return 

99 if not broken.any(): 

100 return contains 

101 

102 # try to run again with a new random vector 

103 # only do it if check_direction isn't specified 

104 # to avoid infinite recursion 

105 if check_direction is None: 

106 # we're going to run the check again in a random direction 

107 new_direction = util.unitize(np.random.random(3) - 0.5) 

108 # do the mask trick again to be able to assign results 

109 mask = inside_aabb.copy() 

110 mask[mask] = broken 

111 

112 contains[mask] = contains_points( 

113 intersector, points[inside_aabb][broken], check_direction=new_direction 

114 ) 

115 

116 constants.log.debug( 

117 "detected %d broken contains test, attempted to fix", broken.sum() 

118 ) 

119 

120 return contains