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