Coverage for /usr/local/lib/python3.10/site-packages/spam/orientations/analyse.py: 95%
141 statements
« prev ^ index » next coverage.py v7.8.0, created at 2025-04-24 17:26 +0000
1import numpy
2import scipy
3import scipy.special
4from scipy.stats import chi2
7def fitVonMisesFisher(orientations, confVMF=None, confMu=None, confKappa=None):
8 """
9 This function fits a vonMises-Fisher distribution to a set of N 3D unit vectors.
10 The distribution is characterised by a mean orientation mu and a spread parameter kappa.
12 Parameters
13 -----------
14 orientations : Nx3 array of floats
15 Z, Y and X components of each vector.
17 confVMF : float
18 Confidence interval for the test on the vMF distribution
19 Used for checking wheter the data can be modeled with a vMF distribution
20 Default = 95%
22 confMu : float
23 Confidence interval for the test on the mean orientation mu
24 Used for computing the error on the mean orientation
25 Default = 95%
27 confKappa : float
28 Confidence interval for the test on kappa
29 Used for computing the error on kappa
30 Default = 95%
32 Returns
33 --------
34 Dictionary containing:
36 Keys:
37 orientations : Nx3 array of floats
38 Z, Y and X components of each vector that is located in the same quadrant as the mean orientation
40 mu : 1x3 array of floats
41 Z, Y and X components of mean orientation.
43 theta : float
44 Inclination angle of the mean orientation in degrees - angle with the Z axis
46 alpha : float
47 Azimuth angle of the mean orientation in degrees - angle in the X-Y plane
49 R : float
50 Mean resultant length
51 First order measure of concentration (ranging between 0 and 1)
53 kappa : int
54 Spread of the distribution, must be > 0.
55 Higher values of kappa mean a higher concentration along the main orientation
57 vectorsProj : Nx3 array of floats
58 Z, Y and X components of each vector projected along the mean orientation
60 fisherTest : bool
61 Boolean representing the result of the test on the vMF distribution
62 1 = The data can be modeled with a vMF distribution
63 0 = The data cannot be modeled with a vMF distribution
65 fisherTestVal : float
66 Value to be compared against the critical value, taken from a Chi-squared distrition
68 muTest : float
69 Error associated to the mean orientation
70 Defined as the semi-vertical angle of the cone that comprises the distribution
72 kappaTest : 1x2 list of floats
73 Maximum and minimum value of kappa, given the confidence interval
75 Notes
76 -----
78 The calculation of kappa implemented from Tanabe, A., et al., (2007). Parameter estimation for von Mises_Fisher distributions. doi: 10.1007/s00180-007-0030-7
80 """
82 # Check that the vectors are 3D
83 assert orientations.shape[1] == 3, "\n spam.orientations.fitVonMisesFisher: The vectors must be an array of Nx3"
85 # If needed, assign confidence intervals
86 if confVMF is None:
87 confVMF = 0.95
88 if confMu is None:
89 confMu = 0.95
90 if confKappa is None:
91 confKappa = 0.95
92 # Check the values of the confidence intervals
93 assert confVMF > 0 and confVMF < 1, "\n spam.orientations.fitVonMisesFisher: The confidence interval for confVMF should be between 0 and 1"
94 assert confMu > 0 and confMu < 1, "\n spam.orientations.fitVonMisesFisher: The confidence interval for confMu should be between 0 and 1"
95 assert confKappa > 0 and confKappa < 1, "\n spam.orientations.fitVonMisesFisher: The confidence interval for confKappa should be between 0 and 1"
97 # Create result dictionary
98 res = {}
99 # Remove possible vectors [0, 0, 0]
100 orientations = orientations[numpy.where(numpy.sum(orientations, axis=1) != 0)[0]]
101 # Normalize all the vectors from http://stackoverflow.com/questions/2850743/numpy-how-to-quickly-normalize-many-vectors
102 norms = numpy.apply_along_axis(numpy.linalg.norm, 1, orientations)
103 orientations = orientations / norms.reshape(-1, 1)
104 # Read Number of Points
105 numberOfPoints = orientations.shape[0]
106 # 1. Get raw-orientation with SVD and flip accordingly
107 vectSVD = meanOrientation(orientations)
108 # Flip accordingly to main direction
109 for i in range(numberOfPoints):
110 vect = orientations[i, :]
111 # Compute angle between both vectors
112 delta1 = numpy.degrees(numpy.arccos((numpy.dot(vectSVD, vect)) / (numpy.linalg.norm(vectSVD) * numpy.linalg.norm(vect))))
113 delta2 = numpy.degrees(numpy.arccos((numpy.dot(vectSVD, -1 * vect)) / (numpy.linalg.norm(vectSVD) * numpy.linalg.norm(-1 * vect))))
114 if delta1 < delta2:
115 orientations[i, :] = vect
116 else:
117 orientations[i, :] = -1 * vect
118 res.update({"orientations": orientations})
119 # 2. Compute parameters of vMF
120 # Compute mean orientation
121 mu = numpy.sum(orientations, axis=0) / numpy.linalg.norm(numpy.sum(orientations, axis=0))
122 res.update({"mu": mu})
123 # Decompose mean orientation into polar coordinates
124 thetaR = numpy.arccos(mu[0])
125 alphaR = numpy.arctan2(mu[1], mu[2])
126 if alphaR < 0:
127 alphaR = alphaR + 2 * numpy.pi
128 res.update({"theta": numpy.degrees(thetaR)})
129 res.update({"alpha": numpy.degrees(alphaR)})
130 # Compute mean resultant length
131 R = numpy.linalg.norm(numpy.sum(orientations, axis=0)) / numberOfPoints
132 res.update({"R": R})
133 # Compute rotation matrix - needed for projecting all vector around mu
134 # Taken from pg 194 from MardiaJupp - Fisher book eq. 3.9 is wrong!
135 rotMatrix = numpy.array(
136 [
137 [
138 numpy.cos(thetaR),
139 numpy.sin(thetaR) * numpy.sin(alphaR),
140 numpy.sin(thetaR) * numpy.cos(alphaR),
141 ],
142 [0, numpy.cos(alphaR), -1 * numpy.sin(alphaR)],
143 [
144 -1 * numpy.sin(thetaR),
145 numpy.cos(thetaR) * numpy.sin(alphaR),
146 numpy.cos(thetaR) * numpy.cos(alphaR),
147 ],
148 ]
149 )
150 # Project vectors - needed for computing kappa
151 orientationsProj = numpy.zeros((numberOfPoints, 3))
152 for i in range(numberOfPoints):
153 orientationsProj[i, :] = rotMatrix.dot(orientations[i, :])
154 if orientationsProj[i, 0] < 0:
155 orientationsProj[i, :] = -1 * orientationsProj[i, :]
156 res.update({"vectorsProj": orientationsProj})
157 # Compute Kappa
158 Z_bar = numpy.sum(orientationsProj[:, 0]) / len(orientationsProj)
159 Y_bar = numpy.sum(orientationsProj[:, 1]) / len(orientationsProj)
160 X_bar = numpy.sum(orientationsProj[:, 2]) / len(orientationsProj)
161 R = numpy.sqrt(Z_bar**2 + Y_bar**2 + X_bar**2)
162 # First Kappa guess
163 k_t = R * (3 - 1) / (1 - R**2)
164 error_i = 5
165 # Main Iteration
166 while error_i > 0.001: # t is step i, T is step i+1
167 I_1 = scipy.special.iv(3 / 2 - 1, k_t)
168 I_2 = scipy.special.iv(3 / 2, k_t)
169 k_T = R * k_t * (I_1 / I_2)
170 error_i = 100 * (numpy.abs(k_T - k_t) / k_t)
171 k_t = k_T.copy()
172 # Add results
173 res.update({"kappa": k_t})
174 # 3. Tests
175 # Test for vMF distribution - Can we really model the data with a vMF distribution?
176 valCritic = scipy.stats.chi2.ppf(1 - confVMF, 3)
177 test = 3 * (R**2) / numberOfPoints
178 if test < valCritic:
179 fisherFit = True
180 else:
181 fisherFit = False
182 res.update({"fisherTest": fisherFit})
183 res.update({"fisherTestVal": test})
184 # Test the location of mu - compute the semi-vertical angle of the cone
185 d = 0
186 for vect in orientations:
187 d += (numpy.sum(vect * mu)) ** 2
188 d = 1 - (1 / numberOfPoints) * d
189 sigma = numpy.sqrt(d / (numberOfPoints * R**2))
190 angle = numpy.degrees(numpy.arcsin(numpy.sqrt(-1 * numpy.log(1 - confMu)) * sigma))
191 res.update({"muTest": angle})
192 # Test the value of Kappa - compute interval for Kappa - eq. 5.37 Fisher
193 kappaDown = 0.5 * chi2.ppf(0.5 * (1 - confKappa), 2 * numberOfPoints - 2) / (numberOfPoints - numberOfPoints * R)
194 kappaUp = 0.5 * chi2.ppf(1 - 0.5 * (1 - confKappa), 2 * numberOfPoints - 2) / (numberOfPoints - numberOfPoints * R)
195 res.update({"kappaTest": [kappaDown, kappaUp]})
196 return res
199def meanOrientation(orientations):
200 """
201 This function performs a Singular Value Decomposition (SVD) on a series of 3D unit vectors to find the main direction of the set
203 Parameters
204 -----------
205 orientations : Nx3 numpy array of floats
206 Z, Y and X components of direction vectors.
207 Non-unit vectors are normalised.
209 Returns
210 --------
211 orientationVector : 1x3 numpy arrayRI*numpy.cos(thetaI) of floats
212 Z, Y and X components of direction vector.
214 Notes
215 -----
216 Implementation taken from https://www.ltu.se/cms_fs/1.51590!/svd-fitting.pdf
218 """
220 # Read Number of Points
221 orientations.shape[0]
222 # Remove possible vectors [0, 0, 0]
223 orientations = orientations[numpy.where(numpy.sum(orientations, axis=1) != 0)[0]]
224 # Normalise all the vectors from http://stackoverflow.com/questions/2850743/numpy-how-to-quickly-normalize-many-vectors
225 norms = numpy.apply_along_axis(numpy.linalg.norm, 1, orientations)
226 orientations = orientations / norms.reshape(-1, 1)
227 # Include the negative part
228 orientationsSVD = numpy.concatenate((orientations, -1 * orientations), axis=0)
229 # Compute the centre
230 meanVal = numpy.mean(orientationsSVD, axis=0)
231 # Center array
232 orientationsCenteredSVD = orientationsSVD - meanVal
233 # Run SVD
234 svd = numpy.linalg.svd(orientationsCenteredSVD.T, full_matrices=False)
235 # Principal direction
236 orientationVector = svd[0][:, 0]
237 # Flip (if needed) to have a positive Z value
238 if orientationVector[0] < 0:
239 orientationVector = -1 * orientationVector
240 return orientationVector
243def fabricTensor(orientations):
244 """
245 Calculation of a second order fabric tensor from 3D unit vectors representing orientations
247 Parameters
248 ----------
249 orientations: Nx3 array of floats
250 Z, Y and X components of direction vectors
251 Non-unit vectors are normalised.
253 Returns
254 -------
255 N: 3x3 array of floats
256 normalised second order fabric tensor
257 with N[0,0] corresponding to z-z, N[1,1] to y-y and N[2,2] x-x
259 F: 3x3 array of floats
260 fabric tensor of the third kind (deviatoric part)
261 with F[0,0] corresponding to z-z, F[1,1] to y-y and F[2,2] x-x
263 a: float
264 scalar anisotropy factor based on the deviatoric part F
266 Note
267 ----
268 see [Kanatani, 1984] for more information on the fabric tensor
269 and [Gu et al, 2017] for the scalar anisotropy factor
271 Function contibuted by Max Wiebicke (Dresden University)
272 """
273 # from http://stackoverflow.com/questions/2850743/numpy-how-to-quickly-normalize-many-vectors
274 norms = numpy.apply_along_axis(numpy.linalg.norm, 1, orientations)
275 orientations = orientations / norms.reshape(-1, 1)
276 # create an empty array
277 N = numpy.zeros((3, 3))
278 size = len(orientations)
279 for i in range(size):
280 orientation = orientations[i]
281 tensProd = numpy.outer(orientation, orientation)
282 N[:, :] = N[:, :] + tensProd
283 # fabric tensor of the first kind
284 N = N / size
285 # fabric tensor of third kind
286 F = (N - (numpy.trace(N) * (1.0 / 3.0)) * numpy.eye(3, 3)) * (15.0 / 2.0)
287 # scalar anisotropy factor
288 a = numpy.sqrt(3.0 / 2.0 * numpy.tensordot(F, F, axes=2))
290 return N, F, a
293def projectOrientation(vector, coordSystem, projectionSystem):
294 """
295 This functions projects a 3D vector from a given coordinate system into a 2D plane given by a defined projection.
297 Parameters
298 ----------
299 vector: 1x3 array of floats
300 Vector to be projected
301 For cartesian system: ZYX
302 For spherical system: r, tetha (inclination), phi (azimuth) in Radians
304 coordSystem: string
305 Coordinate system of the vector
306 Either "cartesian" or "spherical"
308 projectionSystem : string
309 Projection to be used
310 Either "lambert", "stereo" or "equidistant"
312 Returns
313 -------
314 projection_xy: 1x2 array of floats
315 X and Y coordinates of the projected vector
317 projection_theta_r: 1x2 array of floats
318 Theta and R coordinates of the projected vector in radians
320 """
322 projection_xy_local = numpy.zeros(2)
323 projection_theta_r_local = numpy.zeros(2)
325 # Reshape the vector and check for errors in shape
326 try:
327 vector = numpy.reshape(vector, (3, 1))
328 except Exception:
329 print("\n spam.orientations.projectOrientation: The vector must be an array of 1x3")
330 return
332 if coordSystem == "spherical":
333 # unpack vector
335 r, theta, phi = vector
337 x = r * numpy.sin(theta) * numpy.cos(phi)
338 y = r * numpy.sin(theta) * numpy.sin(phi)
339 z = r * numpy.cos(theta)
341 elif coordSystem == "cartesian":
343 # unpack vector
344 z, y, x = vector
345 # we're in cartesian coordinates, (x-y-z mode) Calculate spherical coordinates
346 # passing to 3d spherical coordinates too...
347 # From: https://en.wikipedia.org/wiki/Spherical_coordinate_system
348 # Several different conventions exist for representing the three coordinates, and for the order in which they should be written.
349 # The use of (r, θ, φ) to denote radial distance, inclination (or elevation), and azimuth, respectively,
350 # is common practice in physics, and is specified by ISO standard 80000-2 :2009, and earlier in ISO 31-11 (1992).
351 r = numpy.sqrt(x**2 + y**2 + z**2)
352 theta = numpy.arccos(z / r) # inclination
353 phi = numpy.arctan2(y, x) # azimuth
355 else:
356 print("\n spam.orientations.projectOrientation: Wrong coordinate system")
357 return
359 if projectionSystem == "lambert": # dividing by sqrt(2) so that we're projecting onto a unit circle
360 projection_xy_local[0] = x * (numpy.sqrt(2 / (1 + z)))
361 projection_xy_local[1] = y * (numpy.sqrt(2 / (1 + z)))
363 # sperhical coordinates -- CAREFUL as per this wikipedia page: https://en.wikipedia.org/wiki/Lambert_azimuthal_equal-area_projection
364 # the symbols for inclination and azimuth ARE INVERTED WITH RESPEST TO THE SPHERICAL COORDS!!!
365 projection_theta_r_local[0] = phi
366 # HACK: doing numpy.pi - angle in order for the +z to be projected to 0,0
367 projection_theta_r_local[1] = 2 * numpy.cos((numpy.pi - theta) / 2)
369 # cylindrical coordinates
370 # projection_theta_r_local[0] = phi
371 # projection_theta_r_local[1] = numpy.sqrt( 2.0 * ( 1 + z ) )
373 elif projectionSystem == "stereo":
374 projection_xy_local[0] = x / (1 - z)
375 projection_xy_local[1] = y / (1 - z)
377 # https://en.wikipedia.org/wiki/Stereographic_projection uses a different standard from the page on spherical coord Spherical_coordinate_system
378 projection_theta_r_local[0] = phi
379 # HACK: doing numpy.pi - angle in order for the +z to be projected to 0,0
380 # HACK: doing numpy.pi - angle in order for the +z to be projected to 0,0
381 projection_theta_r_local[1] = numpy.sin(numpy.pi - theta) / (1 - numpy.cos(numpy.pi - theta))
383 elif projectionSystem == "equidistant":
384 # https://en.wikipedia.org/wiki/Azimuthal_equidistant_projection
385 # TODO: To be checked, but this looks like it should -- a straight down projection.
386 projection_xy_local[0] = numpy.sin(phi)
387 projection_xy_local[1] = numpy.cos(phi)
389 projection_theta_r_local[0] = phi
390 projection_theta_r_local[1] = numpy.cos(theta - numpy.pi / 2)
392 else:
393 print("\n spam.orientations.projectOrientation: Wrong projection system")
394 return
396 return projection_xy_local, projection_theta_r_local