Hide keyboard shortcuts

Hot-keys on this page

r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

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

38

39

40

41

42

43

44

45

46

47

48

49

50

51

52

53

54

55

56

57

58

59

60

61

62

63

64

65

66

67

68

69

70

71

72

73

74

75

76

77

78

79

80

81

82

83

84

85

86

87

88

89

90

91

92

93

94

95

96

97

98

99

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

#***************************************************************************** 

# Copyright (C) 2007 Robert Bradshaw <robertwb@math.washington.edu> 

# 

# Distributed under the terms of the GNU General Public License (GPL) 

# 

# This code is distributed in the hope that it will be useful, 

# but WITHOUT ANY WARRANTY; without even the implied warranty of 

# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 

# General Public License for more details. 

# 

# The full text of the GPL is available at: 

# 

# http://www.gnu.org/licenses/ 

#***************************************************************************** 

  

  

# Utility functions for operating on the point_c struct. 

from libc cimport math 

  

  

cdef inline bint point_c_set(point_c* res, P) except -2: 

res.x, res.y, res.z = P[0], P[1], P[2] 

  

cdef inline bint point_c_eq(point_c P, point_c Q): 

return P.x == Q.x and P.y == Q.y and P.z == Q.z 

  

cdef inline int point_c_cmp(point_c P, point_c Q): 

""" 

Lexographic order 

""" 

if P.x == Q.x: 

if P.y == Q.y: 

if P.z == Q.z: 

return 0 

elif P.z < Q.z: 

return -1 

else: 

return 1 

elif P.y < Q.y: 

return -1 

else: 

return 1 

elif P.x < Q.x: 

return -1 

else: 

return 1 

  

cdef inline void point_c_update_finite_lower_bound(point_c* res, point_c P): 

# We use the condition "not P.x > res.x" so that the condition returns True if res.x is NaN 

if math.isfinite(P.x) and not P.x > res.x: 

res.x = P.x 

if math.isfinite(P.y) and not P.y > res.y: 

res.y = P.y 

if math.isfinite(P.z) and not P.z > res.z: 

res.z = P.z 

  

cdef inline void point_c_update_finite_upper_bound(point_c* res, point_c P): 

# We use the condition "not P.x < res.x" so that the condition returns True if res.x is NaN 

if math.isfinite(P.x) and not P.x < res.x: 

res.x = P.x 

if math.isfinite(P.y) and not P.y < res.y: 

res.y = P.y 

if math.isfinite(P.z) and not P.z < res.z: 

res.z = P.z 

  

cdef inline void point_c_lower_bound(point_c* res, point_c P, point_c Q): 

res.x = P.x if P.x < Q.x else Q.x 

res.y = P.y if P.y < Q.y else Q.y 

res.z = P.z if P.z < Q.z else Q.z 

  

cdef inline void point_c_upper_bound(point_c* res, point_c P, point_c Q): 

res.x = P.x if P.x > Q.x else Q.x 

res.y = P.y if P.y > Q.y else Q.y 

res.z = P.z if P.z > Q.z else Q.z 

  

cdef inline void point_c_add(point_c* res, point_c P, point_c Q): 

res.x = P.x + Q.x 

res.y = P.y + Q.y 

res.z = P.z + Q.z 

  

cdef inline void point_c_sub(point_c* res, point_c P, point_c Q): 

res.x = P.x - Q.x 

res.y = P.y - Q.y 

res.z = P.z - Q.z 

  

cdef inline void point_c_mul(point_c* res, point_c P, double a): 

res.x = a * P.x 

res.y = a * P.y 

res.z = a * P.z 

  

cdef inline double point_c_dot(point_c P, point_c Q): 

return P.x*Q.x + P.y*Q.y + P.z*Q.z 

  

cdef inline void point_c_cross(point_c* res, point_c P, point_c Q): 

res.x = P.y * Q.z - P.z * Q.y 

res.y = P.z * Q.x - P.x * Q.z 

res.z = P.x * Q.y - P.y * Q.x 

  

cdef inline double point_c_len(point_c P): 

return math.sqrt(point_c_dot(P, P)) 

  

cdef inline void point_c_transform(point_c* res, double* M, point_c P): 

""" 

M is a flattened 4x4 matrix, row major, representing an Euclidean Transformation. 

Operate on P as a point. 

""" 

res.x = M[0]*P.x + M[1]*P.y + M[2]*P.z + M[3] 

res.y = M[4]*P.x + M[5]*P.y + M[6]*P.z + M[7] 

res.z = M[8]*P.x + M[9]*P.y + M[10]*P.z + M[11] 

  

cdef inline void point_c_stretch(point_c* res, double* M, point_c P): 

""" 

M is a flattened 4x4 matrix, row major, representing an Euclidean Transformation. 

Operate on P as a vector (i.e. ignore the translation component) 

""" 

res.x = M[0]*P.x + M[1]*P.y + M[2]*P.z 

res.y = M[4]*P.x + M[5]*P.y + M[6]*P.z 

res.z = M[8]*P.x + M[9]*P.y + M[10]*P.z 

  

cdef inline void face_c_normal(point_c* res, face_c face, point_c* vlist): 

cdef point_c e1, e2 

point_c_sub(&e1, vlist[face.vertices[0]], vlist[face.vertices[1]]) 

point_c_sub(&e2, vlist[face.vertices[2]], vlist[face.vertices[1]]) 

point_c_cross(res, e1, e2) 

  

cdef inline double cos_face_angle(face_c F, face_c E, point_c* vlist): 

cdef point_c nF, nE 

face_c_normal(&nF, F, vlist) 

face_c_normal(&nE, E, vlist) 

cdef double dot = point_c_dot(nF, nE) 

return dot / math.sqrt(point_c_dot(nF, nF)*point_c_dot(nE, nE)) 

  

cdef inline double sin_face_angle(face_c F, face_c E, point_c* vlist): 

cdef point_c nF, nE 

face_c_normal(&nF, F, vlist) 

face_c_normal(&nE, E, vlist) 

cdef double dot = point_c_dot(nF, nE) 

return math.sqrt(1-(dot*dot)/(point_c_dot(nF, nF)*point_c_dot(nE, nE)))