From 4c96903ab426efa1114ec326dcc83edb5861d3ed Mon Sep 17 00:00:00 2001 From: Patol75 Date: Thu, 29 Feb 2024 04:29:05 +0000 Subject: [PATCH] deploy: 590ac5215f2741f4e18f13aa8f9dcf5a77bda0da --- pydrex/stats.html | 1064 +++++++++++++++++++++++---------------------- 1 file changed, 534 insertions(+), 530 deletions(-) diff --git a/pydrex/stats.html b/pydrex/stats.html index 10aac6b8..e1be26ff 100644 --- a/pydrex/stats.html +++ b/pydrex/stats.html @@ -159,312 +159,314 @@

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

53 # Cumulative volume fractions. 54 frac_ascending = frac[sort_ascending] 55 cumfrac = frac_ascending.cumsum() -56 # Number of new samples with volume less than each cumulative fraction. -57 count_less = np.searchsorted(cumfrac, rng.random(n_samples)) -58 out_orientations[i, ...] = orient[sort_ascending][count_less] -59 out_fractions[i, ...] = frac_ascending[count_less] -60 return out_orientations, out_fractions +56 # Force cumfrac[-1] to be equal to sum(frac_ascending) i.e. 1. +57 cumfrac[-1] = 1.0 +58 # Number of new samples with volume less than each cumulative fraction. +59 count_less = np.searchsorted(cumfrac, rng.random(n_samples)) +60 out_orientations[i, ...] = orient[sort_ascending][count_less] +61 out_fractions[i, ...] = frac_ascending[count_less] +62 return out_orientations, out_fractions @@ -557,50 +561,50 @@

- 
 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)
 

 
 
@@ -637,77 +641,77 @@
- 
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
 

 
 
@@ -734,74 +738,74 @@
- 
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    )
 

 
 
@@ -840,18 +844,18 @@
- 
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
 

 
 
@@ -871,14 +875,14 @@
- 
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)
 

 
 
@@ -898,14 +902,14 @@
- 
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)
 

 
 
@@ -925,12 +929,12 @@
- 
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)
 

 
 
@@ -950,13 +954,13 @@
- 
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)