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