- 79def misorientation_hist(orientations, system: _geo.LatticeSystem, bins=None):
- 80 r"""Calculate misorientation histogram for polycrystal orientations.
- 81
- 82 The `bins` argument is passed to `numpy.histogram`.
- 83 If left as `None`, 1° bins will be used as recommended by the reference paper.
- 84 The `symmetry` argument specifies the lattice system which determines intrinsic
- 85 symmetry degeneracies and the maximum allowable misorientation angle.
- 86 See `_geo.LatticeSystem` for supported systems.
- 87
- 88 .. warning::
- 89 This method must be able to allocate an array of shape
- 90 $ \frac{N!}{2(N-2)!}× M^{2} $
- 91 for N the length of `orientations` and M the number of symmetry operations for
- 92 the given `system`.
- 93
- 94 See [Skemer et al. (2005)](https://doi.org/10.1016/j.tecto.2005.08.023).
+ 81def misorientation_hist(orientations, system: _geo.LatticeSystem, bins=None):
+ 82 r"""Calculate misorientation histogram for polycrystal orientations.
+ 83
+ 84 The `bins` argument is passed to `numpy.histogram`.
+ 85 If left as `None`, 1° bins will be used as recommended by the reference paper.
+ 86 The `symmetry` argument specifies the lattice system which determines intrinsic
+ 87 symmetry degeneracies and the maximum allowable misorientation angle.
+ 88 See `_geo.LatticeSystem` for supported systems.
+ 89
+ 90 .. warning::
+ 91 This method must be able to allocate an array of shape
+ 92 $ \frac{N!}{2(N-2)!}× M^{2} $
+ 93 for N the length of `orientations` and M the number of symmetry operations for
+ 94 the given `system`.
95
- 96 """
- 97 symmetry_ops = _geo.symmetry_operations(system)
- 98 # Compute and bin misorientation angles from orientation data.
- 99 q1_array = np.empty(
-100 (sp.comb(len(orientations), 2, exact=True), len(symmetry_ops), 4)
-101 )
-102 q2_array = np.empty(
-103 (sp.comb(len(orientations), 2, exact=True), len(symmetry_ops), 4)
-104 )
-105 for i, e in enumerate(
-106 it.combinations(Rotation.from_matrix(orientations).as_quat(), 2)
-107 ):
-108 q1, q2 = list(e)
-109 for j, qs in enumerate(symmetry_ops):
-110 if qs.shape == (4, 4): # Reflection, not a proper rotation.
-111 q1_array[i, j] = qs @ q1
-112 q2_array[i, j] = qs @ q2
-113 else:
-114 q1_array[i, j] = _utils.quat_product(qs, q1)
-115 q2_array[i, j] = _utils.quat_product(qs, q2)
-116
-117 _log.debug("calculating misorientations...")
-118 _log.debug("largest array size: %s GB", q1_array.nbytes / 1e9)
-119
-120 misorientations_data = _geo.misorientation_angles(q1_array, q2_array)
-121 θmax = _stats._max_misorientation(system)
-122 return np.histogram(misorientations_data, bins=θmax, range=(0, θmax), density=True)
+ 96 See [Skemer et al. (2005)](https://doi.org/10.1016/j.tecto.2005.08.023).
+ 97
+ 98 """
+ 99 symmetry_ops = _geo.symmetry_operations(system)
+100 # Compute and bin misorientation angles from orientation data.
+101 q1_array = np.empty(
+102 (sp.comb(len(orientations), 2, exact=True), len(symmetry_ops), 4)
+103 )
+104 q2_array = np.empty(
+105 (sp.comb(len(orientations), 2, exact=True), len(symmetry_ops), 4)
+106 )
+107 for i, e in enumerate(
+108 it.combinations(Rotation.from_matrix(orientations).as_quat(), 2)
+109 ):
+110 q1, q2 = list(e)
+111 for j, qs in enumerate(symmetry_ops):
+112 if qs.shape == (4, 4): # Reflection, not a proper rotation.
+113 q1_array[i, j] = qs @ q1
+114 q2_array[i, j] = qs @ q2
+115 else:
+116 q1_array[i, j] = _utils.quat_product(qs, q1)
+117 q2_array[i, j] = _utils.quat_product(qs, q2)
+118
+119 _log.debug("calculating misorientations...")
+120 _log.debug("largest array size: %s GB", q1_array.nbytes / 1e9)
+121
+122 misorientations_data = _geo.misorientation_angles(q1_array, q2_array)
+123 θmax = _stats._max_misorientation(system)
+124 return np.histogram(misorientations_data, bins=θmax, range=(0, θmax), density=True)
125def misorientations_random(low, high, system: _geo.LatticeSystem):
-126 """Get expected count of misorientation angles for an isotropic aggregate.
-127
-128 Estimate the expected number of misorientation angles between grains
-129 that would fall within $($`low`, `high`$)$ in degrees for an aggregate
-130 with randomly oriented grains, where `low` $∈ [0, $`high`$)$,
-131 and `high` is bounded by the maximum theoretical misorientation angle
-132 for the given lattice symmetry system.
-133 See `_geo.LatticeSystem` for supported systems.
-134
-135 """
-136 # TODO: Add cubic system: [Handscomb 1958](https://doi.org/10.4153/CJM-1958-010-0)
-137 max_θ = _max_misorientation(system)
-138 M, N = system.value
-139 if not 0 <= low <= high <= max_θ:
-140 raise ValueError(
-141 f"bounds must obey `low` ∈ [0, `high`) and `high` < {max_θ}.\n"
-142 + f"You've supplied (`low`, `high`) = ({low}, {high})."
-143 )
-144
-145 counts_low = 0 # Number of counts at the lower bin edge.
-146 counts_high = 0 # ... at the higher bin edge.
-147 counts_both = [counts_low, counts_high]
-148
-149 # Some constant factors.
-150 a = np.tan(np.deg2rad(90 / M))
-151 b = 2 * np.rad2deg(np.arctan(np.sqrt(1 + a**2)))
-152 c = round(2 * np.rad2deg(np.arctan(np.sqrt(1 + 2 * a**2))))
-153
-154 for i, edgeval in enumerate([low, high]):
-155 d = np.deg2rad(edgeval)
-156
-157 if 0 <= edgeval <= (180 / M):
-158 counts_both[i] += (N / 180) * (1 - np.cos(d))
-159
-160 elif (180 / M) <= edgeval <= (180 * M / N):
-161 counts_both[i] += (N / 180) * a * np.sin(d)
-162
-163 elif 90 <= edgeval <= b:
-164 counts_both[i] += (M / 90) * ((M + a) * np.sin(d) - M * (1 - np.cos(d)))
-165
-166 elif b <= edgeval <= c:
-167 ν = np.tan(np.deg2rad(edgeval / 2)) ** 2
-168
-169 counts_both[i] = (M / 90) * (
-170 (M + a) * np.sin(d)
-171 - M * (1 - np.cos(d))
-172 + (M / 180)
-173 * (
-174 (1 - np.cos(d))
-175 * (
-176 np.rad2deg(
-177 np.arccos((1 - ν * np.cos(np.deg2rad(180 / M))) / (ν - 1))
-178 )
-179 + 2
-180 * np.rad2deg(
-181 np.arccos(a / (np.sqrt(ν - a**2) * np.sqrt(ν - 1)))
-182 )
-183 )
-184 - 2
-185 * np.sin(d)
-186 * (
-187 2 * np.rad2deg(np.arccos(a / np.sqrt(ν - 1)))
-188 + a * np.rad2deg(np.arccos(1 / np.sqrt(ν - a**2)))
-189 )
-190 )
-191 )
-192 else:
-193 assert False # Should never happen.
-194
-195 return np.sum(counts_both) / 2
+ 127def misorientations_random(low, high, system: _geo.LatticeSystem):
+128 """Get expected count of misorientation angles for an isotropic aggregate.
+129
+130 Estimate the expected number of misorientation angles between grains
+131 that would fall within $($`low`, `high`$)$ in degrees for an aggregate
+132 with randomly oriented grains, where `low` $∈ [0, $`high`$)$,
+133 and `high` is bounded by the maximum theoretical misorientation angle
+134 for the given lattice symmetry system.
+135 See `_geo.LatticeSystem` for supported systems.
+136
+137 """
+138 # TODO: Add cubic system: [Handscomb 1958](https://doi.org/10.4153/CJM-1958-010-0)
+139 max_θ = _max_misorientation(system)
+140 M, N = system.value
+141 if not 0 <= low <= high <= max_θ:
+142 raise ValueError(
+143 f"bounds must obey `low` ∈ [0, `high`) and `high` < {max_θ}.\n"
+144 + f"You've supplied (`low`, `high`) = ({low}, {high})."
+145 )
+146
+147 counts_low = 0 # Number of counts at the lower bin edge.
+148 counts_high = 0 # ... at the higher bin edge.
+149 counts_both = [counts_low, counts_high]
+150
+151 # Some constant factors.
+152 a = np.tan(np.deg2rad(90 / M))
+153 b = 2 * np.rad2deg(np.arctan(np.sqrt(1 + a**2)))
+154 c = round(2 * np.rad2deg(np.arctan(np.sqrt(1 + 2 * a**2))))
+155
+156 for i, edgeval in enumerate([low, high]):
+157 d = np.deg2rad(edgeval)
+158
+159 if 0 <= edgeval <= (180 / M):
+160 counts_both[i] += (N / 180) * (1 - np.cos(d))
+161
+162 elif (180 / M) <= edgeval <= (180 * M / N):
+163 counts_both[i] += (N / 180) * a * np.sin(d)
+164
+165 elif 90 <= edgeval <= b:
+166 counts_both[i] += (M / 90) * ((M + a) * np.sin(d) - M * (1 - np.cos(d)))
+167
+168 elif b <= edgeval <= c:
+169 ν = np.tan(np.deg2rad(edgeval / 2)) ** 2
+170
+171 counts_both[i] = (M / 90) * (
+172 (M + a) * np.sin(d)
+173 - M * (1 - np.cos(d))
+174 + (M / 180)
+175 * (
+176 (1 - np.cos(d))
+177 * (
+178 np.rad2deg(
+179 np.arccos((1 - ν * np.cos(np.deg2rad(180 / M))) / (ν - 1))
+180 )
+181 + 2
+182 * np.rad2deg(
+183 np.arccos(a / (np.sqrt(ν - a**2) * np.sqrt(ν - 1)))
+184 )
+185 )
+186 - 2
+187 * np.sin(d)
+188 * (
+189 2 * np.rad2deg(np.arccos(a / np.sqrt(ν - 1)))
+190 + a * np.rad2deg(np.arccos(1 / np.sqrt(ν - a**2)))
+191 )
+192 )
+193 )
+194 else:
+195 assert False # Should never happen.
+196
+197 return np.sum(counts_both) / 2
213def point_density(
-214 x_data,
-215 y_data,
-216 z_data,
-217 gridsteps=101,
-218 weights=1,
-219 kernel="linear_inverse_kamb",
-220 axial=True,
-221 **kwargs,
-222):
-223 """Estimate point density of orientation data on the unit sphere.
-224
-225 Estimates the density of orientations on the unit sphere by counting the input data
-226 that falls within small areas around a uniform grid of spherical counting locations.
-227 The input data is expected in cartesian coordinates, and the contouring is performed
-228 using kernel functions defined in [Vollmer 1995](https://doi.org/10.1016/0098-3004(94)00058-3).
-229 The following optional parameters control the contouring method:
-230 - `gridsteps` (int) — the number of steps, i.e. number of points along a diameter of
-231 the spherical counting grid
-232 - `weights` (array) — auxiliary weights for each data point
-233 - `kernel` (string) — the name of the kernel function to use, see
-234 `SPHERICAL_COUNTING_KERNELS`
-235 - `axial` (bool) — toggle axial versions of the kernel functions
-236 (for crystallographic data this should normally be kept as `True`)
-237
-238 Any other keyword arguments are passed to the kernel function calls.
-239 Most kernels accept a parameter `σ` to control the degree of smoothing.
-240
-241 """
-242 if kernel not in SPHERICAL_COUNTING_KERNELS:
-243 raise ValueError(f"kernel '{kernel}' is not supported")
-244 weights = np.asarray(weights, dtype=np.float64)
-245
-246 # Create a grid of counters on a cylinder.
-247 ρ_grid, h_grid = np.mgrid[-np.pi : np.pi : gridsteps * 1j, -1 : 1 : gridsteps * 1j]
-248 # Project onto the sphere using the equal-area projection with centre at (0, 0).
-249 λ_grid = ρ_grid
-250 ϕ_grid = np.arcsin(h_grid)
-251 x_counters, y_counters, z_counters = _geo.to_cartesian(
-252 np.pi / 2 - λ_grid.ravel(), np.pi / 2 - ϕ_grid.ravel()
-253 )
-254
-255 # Basically, we can't model this as a convolution as we're not in Euclidean space,
-256 # so we have to iterate through and call the kernel function at each "counter".
-257 data = np.column_stack([x_data, y_data, z_data])
-258 counters = np.column_stack([x_counters, y_counters, z_counters])
-259 totals = np.empty(counters.shape[0])
-260 for i, counter in enumerate(counters):
-261 products = np.dot(data, counter)
-262 if axial:
-263 products = np.abs(products)
-264 density, scale = SPHERICAL_COUNTING_KERNELS[kernel](
-265 products, axial=axial, **kwargs
-266 )
-267 density *= weights
-268 totals[i] = (density.sum() - 0.5) / scale
-269
-270 X_counters, Y_counters = _geo.lambert_equal_area(x_counters, y_counters, z_counters)
+ 215def point_density(
+216 x_data,
+217 y_data,
+218 z_data,
+219 gridsteps=101,
+220 weights=1,
+221 kernel="linear_inverse_kamb",
+222 axial=True,
+223 **kwargs,
+224):
+225 """Estimate point density of orientation data on the unit sphere.
+226
+227 Estimates the density of orientations on the unit sphere by counting the input data
+228 that falls within small areas around a uniform grid of spherical counting locations.
+229 The input data is expected in cartesian coordinates, and the contouring is performed
+230 using kernel functions defined in [Vollmer 1995](https://doi.org/10.1016/0098-3004(94)00058-3).
+231 The following optional parameters control the contouring method:
+232 - `gridsteps` (int) — the number of steps, i.e. number of points along a diameter of
+233 the spherical counting grid
+234 - `weights` (array) — auxiliary weights for each data point
+235 - `kernel` (string) — the name of the kernel function to use, see
+236 `SPHERICAL_COUNTING_KERNELS`
+237 - `axial` (bool) — toggle axial versions of the kernel functions
+238 (for crystallographic data this should normally be kept as `True`)
+239
+240 Any other keyword arguments are passed to the kernel function calls.
+241 Most kernels accept a parameter `σ` to control the degree of smoothing.
+242
+243 """
+244 if kernel not in SPHERICAL_COUNTING_KERNELS:
+245 raise ValueError(f"kernel '{kernel}' is not supported")
+246 weights = np.asarray(weights, dtype=np.float64)
+247
+248 # Create a grid of counters on a cylinder.
+249 ρ_grid, h_grid = np.mgrid[-np.pi : np.pi : gridsteps * 1j, -1 : 1 : gridsteps * 1j]
+250 # Project onto the sphere using the equal-area projection with centre at (0, 0).
+251 λ_grid = ρ_grid
+252 ϕ_grid = np.arcsin(h_grid)
+253 x_counters, y_counters, z_counters = _geo.to_cartesian(
+254 np.pi / 2 - λ_grid.ravel(), np.pi / 2 - ϕ_grid.ravel()
+255 )
+256
+257 # Basically, we can't model this as a convolution as we're not in Euclidean space,
+258 # so we have to iterate through and call the kernel function at each "counter".
+259 data = np.column_stack([x_data, y_data, z_data])
+260 counters = np.column_stack([x_counters, y_counters, z_counters])
+261 totals = np.empty(counters.shape[0])
+262 for i, counter in enumerate(counters):
+263 products = np.dot(data, counter)
+264 if axial:
+265 products = np.abs(products)
+266 density, scale = SPHERICAL_COUNTING_KERNELS[kernel](
+267 products, axial=axial, **kwargs
+268 )
+269 density *= weights
+270 totals[i] = (density.sum() - 0.5) / scale
271
-272 # Normalise to mean, which estimates the density for a "uniform" distribution.
-273 totals /= totals.mean()
-274 totals[totals < 0] = 0
-275 # print(totals.min(), totals.mean(), totals.max())
-276 return (
-277 np.reshape(X_counters, ρ_grid.shape),
-278 np.reshape(Y_counters, ρ_grid.shape),
-279 np.reshape(totals, ρ_grid.shape),
-280 )
+272 X_counters, Y_counters = _geo.lambert_equal_area(x_counters, y_counters, z_counters)
+273
+274 # Normalise to mean, which estimates the density for a "uniform" distribution.
+275 totals /= totals.mean()
+276 totals[totals < 0] = 0
+277 # print(totals.min(), totals.mean(), totals.max())
+278 return (
+279 np.reshape(X_counters, ρ_grid.shape),
+280 np.reshape(Y_counters, ρ_grid.shape),
+281 np.reshape(totals, ρ_grid.shape),
+282 )
296def exponential_kamb(cos_dist, σ=10, axial=True):
-297 """Kernel function from Vollmer 1995 for exponential smoothing."""
-298 n = float(cos_dist.size)
-299 if axial:
-300 f = 2 * (1.0 + n / σ**2)
-301 units = np.sqrt(n * (f / 2.0 - 1) / f**2)
-302 else:
-303 f = 1 + n / σ**2
-304 units = np.sqrt(n * (f - 1) / (4 * f**2))
-305
-306 count = np.exp(f * (cos_dist - 1))
-307 return count, units
+ 298def exponential_kamb(cos_dist, σ=10, axial=True):
+299 """Kernel function from Vollmer 1995 for exponential smoothing."""
+300 n = float(cos_dist.size)
+301 if axial:
+302 f = 2 * (1.0 + n / σ**2)
+303 units = np.sqrt(n * (f / 2.0 - 1) / f**2)
+304 else:
+305 f = 1 + n / σ**2
+306 units = np.sqrt(n * (f - 1) / (4 * f**2))
+307
+308 count = np.exp(f * (cos_dist - 1))
+309 return count, units
310def linear_inverse_kamb(cos_dist, σ=10, axial=True):
-311 """Kernel function from Vollmer 1995 for linear smoothing."""
-312 n = float(cos_dist.size)
-313 radius = _kamb_radius(n, σ, axial=axial)
-314 f = 2 / (1 - radius)
-315 cos_dist = cos_dist[cos_dist >= radius]
-316 count = f * (cos_dist - radius)
-317 return count, _kamb_units(n, radius)
+ 312def linear_inverse_kamb(cos_dist, σ=10, axial=True):
+313 """Kernel function from Vollmer 1995 for linear smoothing."""
+314 n = float(cos_dist.size)
+315 radius = _kamb_radius(n, σ, axial=axial)
+316 f = 2 / (1 - radius)
+317 cos_dist = cos_dist[cos_dist >= radius]
+318 count = f * (cos_dist - radius)
+319 return count, _kamb_units(n, radius)
320def square_inverse_kamb(cos_dist, σ=10, axial=True):
-321 """Kernel function from Vollmer 1995 for inverse square smoothing."""
-322 n = float(cos_dist.size)
-323 radius = _kamb_radius(n, σ, axial=axial)
-324 f = 3 / (1 - radius) ** 2
-325 cos_dist = cos_dist[cos_dist >= radius]
-326 count = f * (cos_dist - radius) ** 2
-327 return count, _kamb_units(n, radius)
+ 322def square_inverse_kamb(cos_dist, σ=10, axial=True):
+323 """Kernel function from Vollmer 1995 for inverse square smoothing."""
+324 n = float(cos_dist.size)
+325 radius = _kamb_radius(n, σ, axial=axial)
+326 f = 3 / (1 - radius) ** 2
+327 cos_dist = cos_dist[cos_dist >= radius]
+328 count = f * (cos_dist - radius) ** 2
+329 return count, _kamb_units(n, radius)
330def kamb_count(cos_dist, σ=10, axial=True):
-331 """Original Kamb 1959 kernel function (raw count within radius)."""
-332 n = float(cos_dist.size)
-333 dist = _kamb_radius(n, σ, axial=axial)
-334 count = (cos_dist >= dist).astype(float)
-335 return count, _kamb_units(n, dist)
+ 332def kamb_count(cos_dist, σ=10, axial=True):
+333 """Original Kamb 1959 kernel function (raw count within radius)."""
+334 n = float(cos_dist.size)
+335 dist = _kamb_radius(n, σ, axial=axial)
+336 count = (cos_dist >= dist).astype(float)
+337 return count, _kamb_units(n, dist)
338def schmidt_count(cos_dist, axial=None):
-339 """Schmidt (a.k.a. 1%) counting kernel function."""
-340 radius = 0.01
-341 count = ((1 - cos_dist) <= radius).astype(float)
-342 # To offset the count.sum() - 0.5 required for the kamb methods...
-343 count = 0.5 / count.size + count
-344 return count, (cos_dist.size * radius)
+ 340def schmidt_count(cos_dist, axial=None):
+341 """Schmidt (a.k.a. 1%) counting kernel function."""
+342 radius = 0.01
+343 count = ((1 - cos_dist) <= radius).astype(float)
+344 # To offset the count.sum() - 0.5 required for the kamb methods...
+345 count = 0.5 / count.size + count
+346 return count, (cos_dist.size * radius)