Coverage for /usr/local/lib/python3.8/site-packages/spam/orientations/generate.py: 100%
103 statements
« prev ^ index » next coverage.py v7.2.3, created at 2023-11-22 13:26 +0000
1import math
3import numpy
4import scipy
5import scipy.special
8def generateIsotropic(N):
9 """
10 There is no analytical solution for putting equally-spaced points on a unit sphere.
11 This Saff and Kuijlaars spiral algorithm gets close.
13 Parameters
14 ----------
15 N : integer
16 Number of points to generate
18 Returns
19 -------
20 orientations : Nx3 numpy array
21 Z,Y,X unit vectors of orientations for each point on sphere
23 Note
24 ----------
25 For references, see:
26 http://www.cgafaq.info/wiki/Evenly_distributed_points_on_sphere
28 Which in turn was based on:
29 http://sitemason.vanderbilt.edu/page/hmbADS
31 From:
32 Rakhmanov, Saff and Zhou: **Minimal Discrete Energy on the Sphere**, Mathematical Research Letters, Vol. 1 (1994), pp. 647-662:
33 https://www.math.vanderbilt.edu/~esaff/texts/155.pdf
35 Also see discussion here:
36 http://groups.google.com/group/sci.math/browse_thread/thread/983105fb1ced42c/e803d9e3e9ba3d23#e803d9e3e9ba3d23%22%22
37 """
39 # Check that it is an integer
40 assert isinstance(N, int), "\n spam.orientations.generateIsotropic: Number of vectors should be an integer"
41 # Check value of number of vectors
42 assert N > 0, "\n spam.orientations.generateIsotropic: Number of vectors should be > 0"
44 M = int(N) * 2
46 s = 3.6 / math.sqrt(M)
48 delta_z = 2 / float(M)
49 z = 1 - delta_z / 2
51 longitude = 0
53 points = numpy.zeros((N, 3))
55 for k in range(N):
56 r = math.sqrt(1 - z * z)
57 points[k, 2] = math.cos(longitude) * r
58 points[k, 1] = math.sin(longitude) * r
59 points[k, 0] = z
60 z = z - delta_z
61 longitude = longitude + s / r
62 return points
65def generateIcosphere(subDiv):
66 """
67 This function creates an unit icosphere (convex polyhedron made from triangles) starting from an icosahedron (polyhedron with 20 faces) and then making subdivision on each triangle.
68 The number of faces is 20*(4**subDiv).
70 Parameters
71 ----------
72 subDiv : integer
73 Number of times that the initial icosahedron is divided.
74 Suggested value: 3
76 Returns
77 -------
78 icoVerts: numberOfVerticesx3 numpy array
79 Coordinates of the vertices of the icosphere
81 icoFaces: numberOfFacesx3 numpy array
82 Indeces of the vertices that compose each face
84 icoVectors: numberOfFacesx3
85 Vectors normal to each face
87 Note
88 ----------
89 From: https://sinestesia.co/blog/tutorials/python-icospheres/
90 """
91 # Chech that it is an integer
92 assert isinstance(subDiv, int), "\n spam.orientations.generateIcosphere: Number of subDiv should be an integer"
93 assert subDiv > 0, print("\n spam.orientations.generateIcosphere: Number of subDiv should be > 0")
95 # 1. Internal functions
97 middle_point_cache = {}
99 def vertex(x, y, z):
100 """Return vertex coordinates fixed to the unit sphere"""
102 length = numpy.sqrt(x**2 + y**2 + z**2)
104 return [i / length for i in (x, y, z)]
106 def middle_point(point_1, point_2):
107 """Find a middle point and project to the unit sphere"""
109 # We check if we have already cut this edge first
110 # to avoid duplicated verts
111 smaller_index = min(point_1, point_2)
112 greater_index = max(point_1, point_2)
114 key = "{}-{}".format(smaller_index, greater_index)
116 if key in middle_point_cache:
117 return middle_point_cache[key]
118 # If it's not in cache, then we can cut it
119 vert_1 = icoVerts[point_1]
120 vert_2 = icoVerts[point_2]
121 middle = [sum(i) / 2 for i in zip(vert_1, vert_2)]
122 icoVerts.append(vertex(middle[0], middle[1], middle[2]))
123 index = len(icoVerts) - 1
124 middle_point_cache[key] = index
125 return index
127 # 2. Create the initial icosahedron
128 # Golden ratio
129 PHI = (1 + numpy.sqrt(5)) / 2
130 icoVerts = [
131 vertex(-1, PHI, 0),
132 vertex(1, PHI, 0),
133 vertex(-1, -PHI, 0),
134 vertex(1, -PHI, 0),
135 vertex(0, -1, PHI),
136 vertex(0, 1, PHI),
137 vertex(0, -1, -PHI),
138 vertex(0, 1, -PHI),
139 vertex(PHI, 0, -1),
140 vertex(PHI, 0, 1),
141 vertex(-PHI, 0, -1),
142 vertex(-PHI, 0, 1),
143 ]
145 icoFaces = [
146 # 5 faces around point 0
147 [0, 11, 5],
148 [0, 5, 1],
149 [0, 1, 7],
150 [0, 7, 10],
151 [0, 10, 11],
152 # Adjacent faces
153 [1, 5, 9],
154 [5, 11, 4],
155 [11, 10, 2],
156 [10, 7, 6],
157 [7, 1, 8],
158 # 5 faces around 3
159 [3, 9, 4],
160 [3, 4, 2],
161 [3, 2, 6],
162 [3, 6, 8],
163 [3, 8, 9],
164 # Adjacent faces
165 [4, 9, 5],
166 [2, 4, 11],
167 [6, 2, 10],
168 [8, 6, 7],
169 [9, 8, 1],
170 ]
172 # 3. Work on the subdivisions
173 for i in range(subDiv):
174 faces_subDiv = []
175 for tri in icoFaces:
176 v1 = middle_point(tri[0], tri[1])
177 v2 = middle_point(tri[1], tri[2])
178 v3 = middle_point(tri[2], tri[0])
179 faces_subDiv.append([tri[0], v1, v3])
180 faces_subDiv.append([tri[1], v2, v1])
181 faces_subDiv.append([tri[2], v3, v2])
182 faces_subDiv.append([v1, v2, v3])
183 icoFaces = faces_subDiv
185 # 4. Compute the normal vector to each face
186 icoVectors = []
187 for tri in icoFaces:
188 # Get the points
189 P1 = numpy.array(icoVerts[tri[0]])
190 P2 = numpy.array(icoVerts[tri[1]])
191 P3 = numpy.array(icoVerts[tri[2]])
192 # Create two vector
193 v1 = P2 - P1
194 v2 = P2 - P3
195 v3 = numpy.cross(v1, v2)
196 norm = vertex(*v3)
197 icoVectors.append(norm)
199 return icoVerts, icoFaces, icoVectors
202def generateVonMisesFisher(mu, kappa, N=1):
203 """
204 This function generates a set of N 3D unit vectors following a vonMises-Fisher distribution, centered at a mean orientation mu and with a spread K.
206 Parameters
207 -----------
208 mu : 1x3 array of floats
209 Z, Y and X components of mean orientation.
210 Non-unit vectors are normalised.
212 kappa : int
213 Spread of the distribution, must be > 0.
214 Higher values of kappa mean a higher concentration along the main orientation
216 N : int
217 Number of vectors to generate
219 Returns
220 --------
221 orientations : Nx3 array of floats
222 Z, Y and X components of each vector.
224 Notes
225 -----
226 Sampling method taken from https://github.com/dlwhittenbury/von-Mises-Fisher-Sampling
228 """
230 def randUniformCircle(N):
231 # N number of orientations
232 v = numpy.random.normal(0, 1, (N, 2))
233 v = numpy.divide(v, numpy.linalg.norm(v, axis=1, keepdims=True))
234 return v
236 def randTmarginal(kappa, N=1):
237 # Start of algorithm
238 b = 2 / (2.0 * kappa + numpy.sqrt(4.0 * kappa**2 + 2**2))
239 x0 = (1.0 - b) / (1.0 + b)
240 c = kappa * x0 + 2 * numpy.log(1.0 - x0**2)
241 orientations = numpy.zeros((N, 1))
242 # Loop over number of orientations
243 for i in range(N):
244 # Continue unil you have an acceptable sample
245 while True:
246 # Sample Beta distribution
247 Z = numpy.random.beta(1, 1)
248 # Sample Uniform distributionNR
249 U = numpy.random.uniform(low=0.0, high=1.0)
250 # W is essentially t
251 W = (1.0 - (1.0 + b) * Z) / (1.0 - (1.0 - b) * Z)
252 # Check whether to accept or reject
253 if kappa * W + 2 * numpy.log(1.0 - x0 * W) - c >= numpy.log(U):
254 # Accept sample
255 orientations[i] = W
256 break
257 return orientations
259 # Check for non-scalar value of kappa
260 assert numpy.isscalar(kappa), "\n spam.orientations.generateVonMisesFisher: kappa should a scalar"
261 assert kappa > 0, "\n spam.orientations.generateVonMisesFisher: kappa should be > 0"
262 assert N > 1, "\n spam.orientations.generateVonMisesFisher: The number of vectors should be > 1"
264 try:
265 mu = numpy.reshape(mu, (1, 3))
266 except Exception:
267 print("\n spam.orientations.generateVonMisesFisher: The main orientation vector must be an array of 1x3")
268 return
269 # Normalize mu
270 mu = mu / numpy.linalg.norm(mu)
271 # check that N > 0!
272 # Array to store orientations
273 orientations = numpy.zeros((N, 3))
274 # Component in the direction of mu (Nx1)
275 t = randTmarginal(kappa, N)
276 # Component orthogonal to mu (Nx(p-1))
277 xi = randUniformCircle(N)
278 # Component in the direction of mu (Nx1).
279 orientations[:, [0]] = t
280 # Component orthogonal to mu (Nx(p-1))
281 Eddy = numpy.sqrt(1 - t**2)
282 orientations[:, 1:] = numpy.hstack((Eddy, Eddy)) * xi
283 # Rotation of orientations to desired mu
284 Olga = scipy.linalg.null_space(mu)
285 Roubin = numpy.concatenate((mu.T, Olga), axis=1)
286 orientations = numpy.dot(Roubin, orientations.T).T
288 return orientations