ModuleGeometry

18def Normalize(Vector):
19    """Normalize Vector."""
20    return Vector / np.linalg.norm(Vector)

24def VectorPerpendicular(Vector):
25    """
26    Find a perpendicular 3D vector in some arbitrary direction
27    """
28
29    if abs(Vector[0]) < 1e-15:
30        return np.array([1, 0, 0])
31    if abs(Vector[1]) < 1e-15:
32        return np.array([0, 1, 0])
33    if abs(Vector[2]) < 1e-15:
34        return np.array([0, 0, 1])
35
36    # set arbitrarily a = b =1
37    return Normalize(np.array([1, 1, -1.0 * (Vector[0] + Vector[1]) / Vector[2]]))

41def AngleBetweenTwoVectors(U, V):
42    """Return the angle in radians between the vectors U and V ; formula from W.Kahan"""
43    u = np.linalg.norm(U)
44    v = np.linalg.norm(V)
45    return 2 * np.arctan2(np.linalg.norm(U * v - V * u), np.linalg.norm(U * v + V * u))

49def IntersectionLinePlane(A, u, P, n):
50    """
51    Return the intersection point between a line and A plane.
52    A is a point of the line, u a vector of the line ; P is a point of the plane, n a normal vector
53    Line's equation : OM = u*t + OA , t a real
54    Plane's equation : n.OP - n.OM = 0
55    """
56    t = np.dot(n, -A + P) / np.dot(u, n)
57    I = u * t + A
58    return I

62def SpiralVogel(NbPoint, Radius):
63    """
64    Return a NbPoint x 2 matrix of 2D points representative of Vogel's spiral with radius Radius
65    """
66
67    GoldenAngle = np.pi * (3 - np.sqrt(5))
68    r = np.sqrt(np.arange(NbPoint) / NbPoint) * Radius
69
70    theta = GoldenAngle * np.arange(NbPoint)
71
72    Matrix = np.zeros((NbPoint, 2))
73    Matrix[:, 0] = np.cos(theta)
74    Matrix[:, 1] = np.sin(theta)
75    Matrix = Matrix * r.reshape((NbPoint, 1))
76
77    return Matrix

81def SolverQuadratic(a, b, c):
82    """
83    Solve the quadratic equation a*x^2 + b*x +c = 0 ; keep only real solutions
84    """
85    Solution = np.roots([a, b, c])
86    RealSolution = []
87
88    for k in range(len(Solution)):
89        if abs(Solution[k].imag) < 1e-15:
90            RealSolution.append(Solution[k].real)
91
92    return RealSolution

 96def SolverQuartic(a, b, c, d, e):
 97    """
 98    Solve the quartic equation a*x^4 + b*x^3 +c*x^2 + d*x + e = 0 ; keep only real solutions
 99    """
100    Solution = np.roots([a, b, c, d, e])
101    RealSolution = []
102
103    for k in range(len(Solution)):
104        if abs(Solution[k].imag) < 1e-15:
105            RealSolution.append(Solution[k].real)
106
107    return RealSolution

111def KeepPositiveSolution(SolutionList):
112    """
113    Keep only positive solution (numbers) in the list
114    """
115    PositiveSolutionList = []
116    epsilon = 1e-12
117    for k in SolutionList:
118        if k > epsilon:
119            PositiveSolutionList.append(k)
120
121    return PositiveSolutionList

125def KeepNegativeSolution(SolutionList):
126    """
127    Keep only positive solution (numbers) in the list
128    """
129    NegativeSolutionList = []
130    epsilon = -1e-12
131    for k in SolutionList:
132        if k < epsilon:
133            NegativeSolutionList.append(k)
134
135    return NegativeSolutionList

139def ClosestPoint(A, I1, I2):
140    """
141    Return the closest point I1 or I2 from A
142    """
143    DistanceSquare1 = np.dot(I1 - A, I1 - A)
144    DistanceSquare2 = np.dot(I2 - A, I2 - A)
145    if DistanceSquare1 < DistanceSquare2:
146        return I1
147    else:
148        return I2

152def FarestPoint(A, I1, I2):
153    """
154    Return the farthest point I1 or I2 from A
155    """
156    DistanceSquare1 = np.dot(I1 - A, I1 - A)
157    DistanceSquare2 = np.dot(I2 - A, I2 - A)
158    if DistanceSquare1 > DistanceSquare2:
159        return I1
160    else:
161        return I2

165def DiameterPointList(PointList):
166    """
167    Return the diameter of the smallest circle (for 2D points) or sphere (3D points) including all the points
168    """
169    if len(PointList) == 0:
170        return None
171
172    elif len(PointList[0]) == 2:
173        Xmax = PointList[0][0]
174        Xmin = PointList[0][0]
175        Ymax = PointList[0][1]
176        Ymin = PointList[0][1]
177
178        for k in PointList:
179            x = k[0]
180            y = k[1]
181
182            Xmax = max(x, Xmax)
183            Xmin = min(x, Xmin)
184            Ymax = max(y, Ymax)
185            Ymin = min(y, Ymin)
186
187        X = abs(Xmax - Xmin)
188        Y = abs(Ymax - Ymin)
189        Diameter = max(X, Y)
190
191        return Diameter
192
193    elif len(PointList[0]) == 3:
194        Xmax = PointList[0][0]
195        Xmin = PointList[0][0]
196        Ymax = PointList[0][1]
197        Ymin = PointList[0][1]
198        Zmax = PointList[0][2]
199        Zmin = PointList[0][2]
200
201        for k in PointList:
202            x = k[0]
203            y = k[1]
204            z = k[2]
205
206            Xmax = max(x, Xmax)
207            Xmin = min(x, Xmin)
208            Ymax = max(y, Ymax)
209            Ymin = min(y, Ymin)
210            Zmax = max(z, Zmax)
211            Zmin = min(z, Zmin)
212
213        X = abs(Xmax - Xmin)
214        Y = abs(Ymax - Ymin)
215        Z = abs(Zmax - Zmin)
216        Diameter = max(X, Y)
217        Diameter = max(Diameter, Z)
218
219        return Diameter

223def CentrePointList(PointList):
224    """
225    Shift all 2D-points in PointList so as to center the point-cloud on the origin [0,0].
226    """
227    ListCoordX = []
228    ListCoordY = []
229    appendX = ListCoordX.append
230    appendY = ListCoordY.append
231
232    for k in PointList:
233        appendX(k[0])
234        appendY(k[1])
235
236    CentreX = (np.amax(ListCoordX) + np.amin(ListCoordX)) * 0.5
237    CentreY = (np.amax(ListCoordY) + np.amin(ListCoordY)) * 0.5
238
239    PointListCentered = []
240    append = PointListCentered.append
241    for k in PointList:
242        x = k[0] - CentreX
243        y = k[1] - CentreY
244        append(np.array([x, y]))
245
246    return PointListCentered

250def IncludeRectangle(X, Y, Point):
251    """
252    Check if a 2D Point is a element of the rectangle X x Y
253    """
254    x = Point[0]
255    y = Point[1]
256    return abs(x) <= abs(X / 2) and abs(y) <= abs(Y / 2)

260def IncludeDisk(R, Point):
261    """
262    Check if a 2D Point is a element of the disk of radius R
263    """
264    x = Point[0]
265    y = Point[1]
266    if (x**2 + y**2) <= R**2:
267        return True
268    else:
269        return False

273def SymmetricalVector(V, SymmetryAxis):
274    """
275    Return the symmetrical vector to V
276    """
277    return RotationAroundAxis(SymmetryAxis, np.pi, V)

284def TranslationPoint(Point, T):
285    """Translate Point by Vector T"""
286    Pointprime = Point + T
287    return Pointprime

291def TranslationPointList(PointList, T):
292    """Translate all points in PointList by Vector T"""
293    PointListPrime = []
294    append = PointListPrime.append
295    for k in PointList:
296        append(TranslationPoint(k, T))
297    return PointListPrime

301def TranslationRay(Ray, T):
302    """Translate a Ray by vector T"""
303    Rayprime = Ray.copy_ray()
304    Rayprime.point = Rayprime.point + T
305    return Rayprime

309def TranslationRayList(RayList, T):
310    """Translate all rays of RayList by vector T"""
311    RayListPrime = []
312    append = RayListPrime.append
313    for k in RayList:
314        append(TranslationRay(k, T))
315    return RayListPrime

322def RotationAroundAxis(Axis, Angle, Vector):
323    """Rotate Vector by Angle (in rad) around Axis"""
324    rot_axis = Normalize(np.array([0.0] + Axis))
325    axis_angle = (Angle * 0.5) * rot_axis
326    vec = quaternion(*Vector)
327    qlog = quaternion(*axis_angle)
328    q = np.exp(qlog)
329    VectorPrime = q * vec * np.conjugate(q)
330    return VectorPrime.imag

334def RotationPoint(Point, Axis1, Axis2):
335    """Transform vector Point such that Axis1 in the old coordinate frame becomes vector Axis 2 in the new coordinate system"""
336    Angle = AngleBetweenTwoVectors(Axis1, Axis2)
337    if abs(Angle) < 1e-10:
338        Pointprime = Point
339    elif abs(Angle - np.pi) < 1e-10:
340        Pointprime = -Point
341    else:
342        N = np.cross(Axis1, Axis2)
343        Pointprime = RotationAroundAxis(N, Angle, Point)
344    return Pointprime

348def RotationPointList(PointList, Axis1, Axis2):
349    """Transform all vectors Point such that Axis1 in the old coordinate frame becomes vector Axis 2 in the new coordinate system"""
350    PointListPrime = []
351    append = PointListPrime.append
352    for k in PointList:
353        append(RotationPoint(k, Axis1, Axis2))
354    return PointListPrime

358def RotationRay(Ray, Axis1, Axis2):
359    """Like RotationPoint but with a Ray object"""
360    Rayprime = Ray.copy_ray()
361    OA = Ray.point
362    u = Ray.vector
363    OB = u + OA
364    OAprime = RotationPoint(OA, Axis1, Axis2)
365    OBprime = RotationPoint(OB, Axis1, Axis2)
366    uprime = OBprime - OAprime
367    Rayprime.point = OAprime
368    Rayprime.vector = uprime
369    return Rayprime

373def RotationRayList(ListeRay, Axis1, Axis2):
374    """Like RotationPointList but with a list of Ray objects"""
375    RayListPrime = []
376    append = RayListPrime.append
377    for k in ListeRay:
378        append(RotationRay(k, Axis1, Axis2))
379    return RayListPrime

383def RotationAroundAxisRayList(ListeRay, Axis, Angle):
384    """Like RotationAroundAxis but with a list of Ray objects. Angle in rad."""
385    RayListPrime = []
386    append = RayListPrime.append
387    for k in ListeRay:
388        Ray = k.copy_ray()
389        Ray.vector = RotationAroundAxis(Axis, Angle, Ray.vector)
390        append(Ray)
391
392    return RayListPrime