diff --git a/orix/measure/pole_density_function.py b/orix/measure/pole_density_function.py index d9613660..af2aae89 100644 --- a/orix/measure/pole_density_function.py +++ b/orix/measure/pole_density_function.py @@ -18,10 +18,11 @@ # import numpy as np -from scipy.ndimage import gaussian_filter +from scipy.interpolate import RegularGridInterpolator from orix.projections.stereographic import StereographicProjection -from orix.quaternion.symmetry import Symmetry +from orix.quaternion.symmetry import C1, Symmetry +from orix.sampling.S2_sampling import _sample_S2_uv_mesh_coordinates from orix.vector.vector3d import Vector3d @@ -31,15 +32,16 @@ def pole_density_function( sigma: float = 5, weights: np.ndarray | None = None, hemisphere: str = "upper", - symmetry: Symmetry | None = None, + symmetry: Symmetry | None = C1, log: bool = False, mrd: bool = True, ) -> tuple[np.ma.MaskedArray, tuple[np.ndarray, np.ndarray]]: - """Compute the Pole Density Function (PDF) of vectors in the - stereographic projection. See :cite:`rohrer2004distribution`. + """Compute the Pole Density Function (PDF) of vectors on a + cubed-sphere grid and return a map in polar coordinates for + plotting. - If ``symmetry`` is defined then the PDF is folded back into the - point group fundamental sector and accumulated. + If ``symmetry`` is defined then the PDF is symmetrizes and the + density map is returned only on the fundamental sector. Parameters ---------- @@ -85,13 +87,15 @@ def pole_density_function( orix.plot.StereographicPlot.pole_density_function orix.vector.Vector3d.pole_density_function """ - from orix.sampling.S2_sampling import _sample_S2_equal_area_coordinates hemisphere = hemisphere.lower() - poles = {"upper": -1, "lower": 1} sp = StereographicProjection(poles[hemisphere]) + # If user explicitly passes symmetry=None + if symmetry is None: + symmetry = C1 + if len(args) == 1: v = args[0] if not isinstance(v, Vector3d): @@ -102,123 +106,332 @@ def pole_density_function( elif len(args) == 2: # azimuth and polar angles v = Vector3d.from_polar(*args) + else: raise ValueError( - "Accepts only one (Vector3d) or two (azimuth, polar) input arguments." + "Accepts only one (Vector3d) or two (azimuth, polar)\ + input arguments." + ) + + if v.size == 0: + raise ValueError + + if np.any(v.norm) == 0.0: + raise ValueError + + # IF we blur by a lot, save some compute time by doing the histograms + # at lower resolution. + if resolution < 0.2 * sigma and resolution < 1.0: + histogram_resolution = 0.2 * sigma + plot_resolution = resolution + else: + histogram_resolution = resolution + plot_resolution = resolution + + # Do actual histogramming + bins = int(90 / histogram_resolution) + bin_edges = np.linspace(-np.pi / 4, np.pi / 4, bins + 1) + hist = np.zeros((6, bins, bins)) + + # Explicit symmetrization + for rotation in symmetry: + + face_index_array, face_coordinates = _cube_gnom_coordinates(rotation * v) + face_index_array = face_index_array.ravel() + face_coordinate_1 = face_coordinates[0].ravel() + face_coordinate_2 = face_coordinates[1].ravel() + for face_index in range(6): + this_face = face_index_array == face_index + w = weights[this_face] if weights is not None else None + hist[face_index] += ( + np.histogramdd( + (face_coordinate_1[this_face], face_coordinate_2[this_face]), + bins=(bin_edges, bin_edges), + weights=w, + )[0] + / symmetry.size + ) + + # Bins are not all same solid angle area, so we need to normalize. + if mrd: + bin_middles = np.linspace( + -np.pi / 4 + np.pi / 4 / bins, np.pi / 4 - np.pi / 4 / bins, bins + ) + y_ang, z_ang = np.meshgrid(bin_middles, bin_middles) + solid_angle_term = ( + 1 + / (np.tan(y_ang) ** 2 + np.tan(z_ang) ** 2 + 1) + / (np.cos(y_ang) * np.cos(z_ang)) + / (1 - 0.5 * (np.sin(z_ang) * np.sin(y_ang)) ** 2) ) + solid_angle_term *= 1 / 6 / np.sum(solid_angle_term) + if weights is not None: + solid_angle_term *= np.sum(weights) + else: + solid_angle_term *= v.size + hist = hist / solid_angle_term[np.newaxis, ...] - if symmetry is not None: - v = v.in_fundamental_sector(symmetry) - # To help with aliasing after reprojection into point group - # fundamental sector in the inverse pole figure case, the - # initial sampling is performed at half the angular resolution - resolution /= 2 - - azimuth, polar, _ = v.to_polar() - # np.histogram2d expects 1d arrays - azimuth, polar = np.ravel(azimuth), np.ravel(polar) - if not azimuth.size: - raise ValueError("`azimuth` and `polar` angles have 0 size.") - - # Generate equal area mesh on S2 - azimuth_coords, polar_coords = _sample_S2_equal_area_coordinates( - resolution, + # Smoothing + if sigma != 0.0: + # If smoothing is only a bit, du 60 small steps, + # otherwise do a max step-size. + if (sigma / histogram_resolution) ** 2 <= 20: + N = 60 + t = 1 / N * (sigma / histogram_resolution) ** 2 + else: + t = 1 / 3 + N = int(1 / t * (sigma / histogram_resolution) ** 2) + + hist = _smooth_gnom_cube_histograms(hist, t, N) + + # Make plot grid + azimuth_coords, polar_coords = _sample_S2_uv_mesh_coordinates( + plot_resolution, hemisphere=hemisphere, azimuth_endpoint=True, ) - azimuth_grid, polar_grid = np.meshgrid(azimuth_coords, polar_coords, indexing="ij") - # generate histogram in angular space - hist, *_ = np.histogram2d( - azimuth, - polar, - bins=(azimuth_coords, polar_coords), - density=False, - weights=weights, - ) - - # "wrap" along azimuthal axis, "reflect" along polar axis - mode = ("wrap", "reflect") - # apply broadening in angular space - hist = gaussian_filter(hist, sigma / resolution, mode=mode) - - # In the case of inverse pole figure, accumulate all values outside - # of the point group fundamental sector back into correct bin within - # fundamental sector - if symmetry is not None: - # compute histogram bin centers in azimuth and polar coords - azimuth_center_grid, polar_center_grid = np.meshgrid( - azimuth_coords[:-1] + np.diff(azimuth_coords) / 2, - polar_coords[:-1] + np.diff(polar_coords) / 2, - indexing="ij", - ) - v_center_grid = Vector3d.from_polar( - azimuth=azimuth_center_grid, polar=polar_center_grid - ).unit - # fold back in into fundamental sector - v_center_grid_fs = v_center_grid.in_fundamental_sector(symmetry) - azimuth_center_fs, polar_center_fs, _ = v_center_grid_fs.to_polar() - azimuth_center_fs = azimuth_center_fs.ravel() - polar_center_fs = polar_center_fs.ravel() - - # Generate coorinates with user-defined resolution. - # When `symmetry` is defined, the initial grid was calculated - # with `resolution = resolution / 2` - azimuth_coords_res2, polar_coords_res2 = _sample_S2_equal_area_coordinates( - 2 * resolution, - hemisphere=hemisphere, - azimuth_endpoint=True, - ) - azimuth_res2_grid, polar_res2_grid = np.meshgrid( - azimuth_coords_res2, polar_coords_res2, indexing="ij" - ) - v_res2_grid = Vector3d.from_polar( - azimuth=azimuth_res2_grid, polar=polar_res2_grid - ) + azimuth_grid, polar_grid = np.meshgrid( + azimuth_coords, + polar_coords, + indexing="ij", + ) + azimuth_center_grid, polar_center_grid = np.meshgrid( + azimuth_coords[:-1] + np.diff(azimuth_coords) / 2, + polar_coords[:-1] + np.diff(polar_coords) / 2, + indexing="ij", + ) - # calculate histogram values for vectors folded back into - # fundamental sector - i = np.digitize(azimuth_center_fs, azimuth_coords_res2[1:-1]) - j = np.digitize(polar_center_fs, polar_coords_res2[1:-1]) - # recompute histogram - temp = np.zeros((azimuth_coords_res2.size - 1, polar_coords_res2.size - 1)) - # add hist data to new histogram without buffering - np.add.at(temp, (i, j), hist.ravel()) - - # get new histogram bins centers to compute histogram mask - azimuth_center_res2_grid, polar_center_res2_grid = np.meshgrid( - azimuth_coords_res2[:-1] + np.ediff1d(azimuth_coords_res2) / 2, - polar_coords_res2[:-1] + np.ediff1d(polar_coords_res2) / 2, - indexing="ij", - ) - v_center_res2_grid = Vector3d.from_polar( - azimuth=azimuth_center_res2_grid, polar=polar_center_res2_grid - ).unit + v_grid = Vector3d.from_polar( + azimuth=azimuth_center_grid, polar=polar_center_grid + ).unit + v_grid_vertexes = Vector3d.from_polar(azimuth=azimuth_grid, polar=polar_grid).unit - # compute histogram data array as masked array - hist = np.ma.array( - temp, mask=~(v_center_res2_grid <= symmetry.fundamental_sector) + # Interpolation from histograms to plot grid + grid_face_index, grid_face_coords = _cube_gnom_coordinates(v_grid) + hist_grid = np.zeros(v_grid.shape) + bin_middles = np.linspace( + -np.pi / 4 + np.pi / 4 / bins, np.pi / 4 - np.pi / 4 / bins, bins + ) + for face_index in range(6): + interpolator = RegularGridInterpolator( + (bin_middles, bin_middles), + hist[face_index], + bounds_error=False, + fill_value=None, ) - # calculate bin vertices - x, y = sp.vector2xy(v_res2_grid) - x, y = x.reshape(v_res2_grid.shape), y.reshape(v_res2_grid.shape) - else: - # all points valid in stereographic projection - hist = np.ma.array(hist, mask=np.zeros_like(hist, dtype=bool)) - # calculate bin vertices - v_grid = Vector3d.from_polar(azimuth=azimuth_grid, polar=polar_grid).unit - x, y = sp.vector2xy(v_grid) - x, y = x.reshape(v_grid.shape), y.reshape(v_grid.shape) - - # Normalize by the average number of counts per cell on the - # unit sphere to calculate in terms of Multiples of Random - # Distribution (MRD). See :cite:`rohrer2004distribution`. - # as `hist` is a masked array, only valid (unmasked) values are - # used in this computation - if mrd: - hist = hist / hist.mean() + this_face = grid_face_index == face_index + hist_grid[this_face] = interpolator(grid_face_coords[:, this_face].T) + hist = hist_grid + + # Mask out points outside funamental region. + mask = ~(v_grid <= symmetry.fundamental_sector) + hist = np.ma.array(hist, mask=mask) + + # Transform grdi to mystery coordinates used by plotting routine + x, y = sp.vector2xy(v_grid_vertexes) + x, y = x.reshape(v_grid_vertexes.shape), y.reshape(v_grid_vertexes.shape) if log: # +1 to avoid taking the log of 0 hist = np.log(hist + 1) return hist, (x, y) + + +def _cube_gnom_coordinates( + vectors: Vector3d, +) -> tuple[np.ndarray[int], np.ndarray[float]]: + """Assigns an index (0 to 5) to an array of ```Vector3d``` + assigning them each to a face of a cube in the followiung way: + + Index 0 is positive x face. (Includes +x+y and +x+z edges.) + Index 1 is negative x face. (Includes -x-y and -x-z edges.) + Index 2 is positive y face. (Includes +y+z and +y-x edges.) + Index 3 is negative y face. (Includes -y-z and -y+x edges.) + Index 4 is positive z face. (Includes -y+z and -x+z edges.) + Index 5 is negative z face. (Includes +y-z and +x-z edges.) + + The two "extra" corners are assigned to Index 0 and Index 1 respectively. + + Then each point is given 2D coordinates on the respective + faces. Coordinates are in angles wrt. cube face center-lines. + Always ordered with increasing coordinates. + First coordinate comes first x -> y -> z. + """ + + if np.any(vectors.norm == 0.0): + raise ValueError + + # Assign face index to each vector + face_index = np.zeros(vectors.shape, dtype=int) + + indx = np.all( + [ + vectors.x >= vectors.y, + vectors.x >= -vectors.y, + vectors.x >= vectors.z, + vectors.x >= -vectors.z, + ], + axis=0, + ) + face_index[indx] = 0 + + indx = np.all( + [ + vectors.x <= vectors.y, + vectors.x <= -vectors.y, + vectors.x <= vectors.z, + vectors.x <= -vectors.z, + ], + axis=0, + ) + face_index[indx] = 1 + + indx = np.all( + [ + vectors.y > vectors.x, + vectors.y >= -vectors.x, + vectors.y >= vectors.z, + vectors.y > -vectors.z, + ], + axis=0, + ) + face_index[indx] = 2 + + indx = np.all( + [ + vectors.y < vectors.x, + vectors.y <= -vectors.x, + vectors.y <= vectors.z, + vectors.y < -vectors.z, + ], + axis=0, + ) + face_index[indx] = 3 + + indx = np.all( + [ + vectors.z > vectors.x, + vectors.z >= -vectors.x, + vectors.z > vectors.y, + vectors.z >= -vectors.y, + ], + axis=0, + ) + face_index[indx] = 4 + + indx = np.all( + [ + vectors.z < vectors.x, + vectors.z <= -vectors.x, + vectors.z < vectors.y, + vectors.z <= -vectors.y, + ], + axis=0, + ) + face_index[indx] = 5 + + # Assign coordinates + coordinates = np.zeros((2,) + vectors.shape) + unit_vectors = vectors.unit + + # Comment: no need for np.arctan2. We know denom is pos + # so np.arctan should be faster. + this_face = face_index == 0 + coordinates[0, this_face] = np.arctan( + unit_vectors.y[this_face] / unit_vectors.x[this_face] + ) + coordinates[1, this_face] = np.arctan( + unit_vectors.z[this_face] / unit_vectors.x[this_face] + ) + + this_face = face_index == 1 + coordinates[0, this_face] = np.arctan( + -unit_vectors.y[this_face] / unit_vectors.x[this_face] + ) + coordinates[1, this_face] = np.arctan( + -unit_vectors.z[this_face] / unit_vectors.x[this_face] + ) + + this_face = face_index == 2 + coordinates[0, this_face] = np.arctan( + unit_vectors.x[this_face] / unit_vectors.y[this_face] + ) + coordinates[1, this_face] = np.arctan( + unit_vectors.z[this_face] / unit_vectors.y[this_face] + ) + + this_face = face_index == 3 + coordinates[0, this_face] = np.arctan( + -unit_vectors.x[this_face] / unit_vectors.y[this_face] + ) + coordinates[1, this_face] = np.arctan( + -unit_vectors.z[this_face] / unit_vectors.y[this_face] + ) + + this_face = face_index == 4 + coordinates[0, this_face] = np.arctan( + unit_vectors.x[this_face] / unit_vectors.z[this_face] + ) + coordinates[1, this_face] = np.arctan( + unit_vectors.y[this_face] / unit_vectors.z[this_face] + ) + + this_face = face_index == 5 + coordinates[0, this_face] = np.arctan( + -unit_vectors.x[this_face] / unit_vectors.z[this_face] + ) + coordinates[1, this_face] = np.arctan( + -unit_vectors.y[this_face] / unit_vectors.z[this_face] + ) + + return face_index, coordinates + + +def _smooth_gnom_cube_histograms( + histograms: np.ndarray[float], + step_parameter: float, + iterations: int = 1, +) -> np.ndarray[float]: + """Histograms shape is (6, n_nbins, n_bins) and edge connectivity + is as according to the rest of this file. + """ + output_histogram = np.copy(histograms) + diffused_weight = np.zeros(histograms.shape) + + for n in range(iterations): + + diffused_weight[...] = 0 + + # Diffuse on faces + for fi in range(6): + diffused_weight[fi, 1:, :] += output_histogram[fi, :-1, :] + diffused_weight[fi, :-1, :] += output_histogram[fi, 1:, :] + diffused_weight[fi, :, 1:] += output_histogram[fi, :, :-1] + diffused_weight[fi, :, :-1] += output_histogram[fi, :, 1:] + + connected_edge_pairs = ( + ((2, slice(None), -1), (4, slice(None), -1)), # +y+z + ((3, slice(None), -1), (4, slice(None), 0)), # -y+z + ((2, slice(None), 0), (5, slice(None), -1)), # +y-z + ((3, slice(None), 0), (5, slice(None), 0)), # -y-z + ((0, slice(None), -1), (4, -1, slice(None))), # +x+z + ((1, slice(None), -1), (4, 0, slice(None))), # -x+z + ((0, slice(None), 0), (5, -1, slice(None))), # +x-z + ((1, slice(None), 0), (5, 0, slice(None))), # -x-z + ((0, -1, slice(None)), (2, -1, slice(None))), # +x+y + ((1, -1, slice(None)), (2, 0, slice(None))), # -x+y + ((0, 0, slice(None)), (3, -1, slice(None))), # +x-y + ((1, 0, slice(None)), (3, 0, slice(None))), # -x-y + ) + + for edge_1, edge_2 in connected_edge_pairs: + diffused_weight[edge_1] += output_histogram[edge_2] + diffused_weight[edge_2] += output_histogram[edge_1] + + # Add to output + output_histogram = ( + 1 - step_parameter + ) * output_histogram + diffused_weight / 4 * step_parameter + + return output_histogram diff --git a/orix/sampling/S2_sampling.py b/orix/sampling/S2_sampling.py index a0c24518..7f5c4700 100644 --- a/orix/sampling/S2_sampling.py +++ b/orix/sampling/S2_sampling.py @@ -392,6 +392,98 @@ def sample_S2_cube_mesh( return Vector3d(np.vstack((bottom, top, east, west, south, north, m_c))).unit +def sample_S2_equiangle_cube_mesh_vertexes( + resolution: float, grid_type: str = "spherified_corner" +) -> Vector3d: + """Return vectors of the vertices of an equi-angle gnomonic mesh on + a unit sphere *S2*. Same grid as ``sample_S2_cube_mesh`` with + ``grid_type="spherified_edge"`` but includes duplicate entries at edges and + corners. + + Parameters + ---------- + resolution + Angle between neighbour grid points, in degrees. + grid_type + Type of cube grid: ``"normalized"``, ``"spherified_edge"`` or + ``"spherified_corner"``. + + Returns + ------- + vec + 3D array of Vectors that sample the unit sphere. First dimension is + the index of the cube face. Second two form a rectangular grid on + the respective face. + """ + + bins_along_edge = int(90 / resolution) + bin_edges = np.linspace(-np.pi / 4, np.pi / 4, bins_along_edge + 1) + + dim_1, dim_2 = np.meshgrid(np.tan(bin_edges), np.tan(bin_edges)) + norm_term = np.sqrt(dim_1**2 + dim_2**2 + 1) + dim_1 /= norm_term + dim_2 /= norm_term + dim_3 = np.ones((bins_along_edge + 1, bins_along_edge + 1)) / norm_term + + data = np.zeros((6, bins_along_edge + 1, bins_along_edge + 1, 3)) + + # Assign points on each edge + data[0, ...] = np.stack([dim_3, dim_2, dim_1], axis=-1) + data[1, ...] = np.stack([-dim_3, dim_2, dim_1], axis=-1) + data[2, ...] = np.stack([dim_2, dim_3, dim_1], axis=-1) + data[3, ...] = np.stack([dim_2, -dim_3, dim_1], axis=-1) + data[4, ...] = np.stack([dim_2, dim_1, dim_3], axis=-1) + data[5, ...] = np.stack([dim_2, dim_1, -dim_3], axis=-1) + + return Vector3d(data=data) + + +def sample_S2_equiangle_cube_mesh_face_centers( + resolution: float, grid_type: str = "spherified_corner" +) -> Vector3d: + """Return vectors of the face_centers of an equi-angle gnomonic mesh on + a unit sphere *S2*. Similar to ``sample_S2_cube_mesh`` but does not include + edges and corners. + + Parameters + ---------- + resolution + Angle between neighbour grid points, in degrees. + grid_type + Type of cube grid: ``"normalized"``, ``"spherified_edge"`` or + ``"spherified_corner"``. + + Returns + ------- + vec + 3D array of Vectors that sample the unit sphere. First dimension is + the index of the cube face. Second two form a rectangular grid on + the respective face. + """ + + bins = int(90 / resolution) + bin_middles = np.linspace( + -np.pi / 4 + np.pi / 4 / bins, np.pi / 4 - np.pi / 4 / bins, bins + ) + dim_1, dim_2 = np.meshgrid(np.tan(bin_middles), np.tan(bin_middles)) + norm_term = np.sqrt(dim_1**2 + dim_2**2 + 1) + dim_1 /= norm_term + dim_2 /= norm_term + dim_3 = np.ones((bins, bins)) / norm_term + + data = np.zeros((6, bins, bins, 3)) + + # Assign points on each edge + data[0, ...] = np.stack([dim_3, dim_2, dim_1], axis=-1) + data[1, ...] = np.stack([-dim_3, dim_2, dim_1], axis=-1) + data[2, ...] = np.stack([dim_2, dim_3, dim_1], axis=-1) + data[3, ...] = np.stack([dim_2, -dim_3, dim_1], axis=-1) + data[4, ...] = np.stack([dim_2, dim_1, dim_3], axis=-1) + data[5, ...] = np.stack([dim_2, dim_1, -dim_3], axis=-1) + + return Vector3d(data=data) + + def sample_S2_hexagonal_mesh(resolution: float) -> Vector3d: """Return vectors of a hexagonal bipyramid mesh projected on a unit sphere *S2*. diff --git a/orix/tests/test_measure/test_measure.py b/orix/tests/test_measure/test_measure.py index 0145de1e..97c9a37d 100644 --- a/orix/tests/test_measure/test_measure.py +++ b/orix/tests/test_measure/test_measure.py @@ -17,18 +17,23 @@ # along with orix. If not, see . # -from copy import deepcopy - import numpy as np import pytest from orix.measure import pole_density_function +from orix.measure.pole_density_function import _cube_gnom_coordinates from orix.quaternion import symmetry +from orix.sampling.S2_sampling import ( + _sample_S2_uv_mesh_coordinates, + sample_S2_equiangle_cube_mesh_face_centers, + sample_S2_random_mesh, +) from orix.vector import Vector3d @pytest.fixture( params=[ + None, symmetry.D2h, symmetry.S6, symmetry.D3d, @@ -45,13 +50,12 @@ def point_groups(request): class TestMeasurePoleDensityFunction: - def test_pole_density_function(self): - v = Vector3d.random(10_000) + def test_output_format(self): + v = sample_S2_random_mesh(1, seed=954) hist1, (x1, y1) = pole_density_function(v) assert hist1.shape[0] + 1 == x1.shape[0] == y1.shape[0] assert hist1.shape[1] + 1 == x1.shape[1] == y1.shape[1] - assert np.allclose(hist1.mean(), 1, rtol=1e-3) assert isinstance(hist1, np.ma.MaskedArray) assert hist1.mask.sum() == 0 @@ -64,60 +68,56 @@ def test_pole_density_function(self): assert isinstance(hist2, np.ma.MaskedArray) assert hist2.mask.sum() > 0 - assert np.allclose(hist1.mean(), hist2.mean()) + hist3, (x3, y3) = pole_density_function(v, hemisphere="lower") + assert hist3.shape[0] + 1 == x3.shape[0] == y3.shape[0] + assert hist3.shape[1] + 1 == x3.shape[1] == y3.shape[1] + assert isinstance(hist3, np.ma.MaskedArray) + assert hist3.mask.sum() == 0 - def test_pole_density_function_symmetry(self, point_groups): + @pytest.mark.parametrize("resolution", [0.5, 1.0]) + def test_pole_density_function_mrd_norm(self, point_groups, resolution): pg = point_groups - v = Vector3d.random(10_000) - - hist, _ = pole_density_function(v, symmetry=pg, mrd=False) - assert np.allclose(hist.sum(), v.size, rtol=0.01) - - def test_pole_density_function_hemisphere(self): - v = Vector3d.random(11_234) - - hist1_upper, _ = pole_density_function(v, hemisphere="upper") - assert np.allclose(hist1_upper.mean(), 1) - - hist1_lower, _ = pole_density_function(v, hemisphere="lower") - assert np.allclose(hist1_lower.mean(), 1) - - hist2_upper, _ = pole_density_function(v, hemisphere="upper", mrd=False) - hist2_lower, _ = pole_density_function(v, hemisphere="lower", mrd=False) - assert np.allclose(hist2_upper.sum() + hist2_lower.sum(), v.size) - - @pytest.mark.parametrize("n", [10, 1000, 10_000, 12_546]) - def test_pole_density_function_values(self, n): - # vectors only on upper hemisphere - v = Vector3d.random(n) - v2 = deepcopy(v) - v2.z[v2.z < 0] *= -1 - - hist1, _ = pole_density_function(v2, mrd=False) - assert np.allclose(hist1.sum(), n, atol=0.1) - - hist2, _ = pole_density_function(v, symmetry=symmetry.Th, mrd=False) - assert np.allclose(hist2.sum(), n, rtol=0.1) - - hist3, _ = pole_density_function(v2, symmetry=symmetry.Th, mrd=False) - assert np.allclose(hist3.sum(), n, rtol=0.1) + v = sample_S2_random_mesh(1.0, seed=230) + + # Make plot grid + _, polar_coords = _sample_S2_uv_mesh_coordinates( + resolution, + hemisphere="upper", + azimuth_endpoint=True, + ) + polar_coords = polar_coords[:-1] + np.diff(polar_coords) / 2 + solid_angle = np.abs(np.sin(polar_coords))[np.newaxis, :] + + for sigma in [2 * resolution, 5 * resolution]: + hist, _ = pole_density_function( + v, + symmetry=pg, + resolution=resolution, + sigma=sigma, + ) + mean_value = np.sum(solid_angle * hist) / np.sum(solid_angle * ~hist.mask) + print(mean_value) + assert np.allclose(mean_value, 1.0, rtol=0.01) def test_pole_density_function_log(self): - v = Vector3d.random(11_234) + # v = Vector3d.random(11_234) + v = sample_S2_random_mesh(1.0, seed=230) hist1, _ = pole_density_function(v, log=False) hist2, _ = pole_density_function(v, log=True) assert not np.allclose(hist1, hist2) def test_pole_density_function_sigma(self): - v = Vector3d.random(11_234) + # v = Vector3d.random(11_234) + v = sample_S2_random_mesh(1.0, seed=230) hist1, _ = pole_density_function(v, sigma=2.5) hist2, _ = pole_density_function(v, sigma=5) assert not np.allclose(hist1, hist2) def test_pole_density_function_weights(self): - v = Vector3d.random(11_234) + # v = Vector3d.random(11_234) + v = sample_S2_random_mesh(1.0, seed=230) v.z[v.z < 0] *= -1 hist0, _ = pole_density_function(v, weights=None) @@ -130,10 +130,12 @@ def test_pole_density_function_weights(self): # the same because MRD normalizes by average assert np.allclose(hist0, hist2) - hist0_counts, _ = pole_density_function(v, weights=None, mrd=False) - hist2_counts, _ = pole_density_function(v, weights=weights2, mrd=False) + # TDOD: Ask about expected behaviour for mrd=False + # hist0_counts, _ = pole_density_function(v, weights=None, mrd=False) + # hist2_counts, _ = pole_density_function(v, weights=weights2, + # mrd=False) # not the same because hist values are not normalized - assert not np.allclose(hist0_counts, hist2_counts) + # assert not np.allclose(hist0_counts, hist2_counts) # non-uniform weights weights2[54] *= 1.01 @@ -154,7 +156,105 @@ def test_pole_density_function_empty_vector_raises(self): v = Vector3d.empty() assert not v.size - with pytest.raises( - ValueError, match="`azimuth` and `polar` angles have 0 size" - ): + with pytest.raises(ValueError): pole_density_function(v) + + +class TestGnomCubeRoutines: + + def test_corner_edge_assignment(self): + """Make sure we get useable results for corner-cases.""" + corners = Vector3d( + [ + [1, 1, 1], + [1, 1, -1], + [1, -1, 1], + [1, -1, -1], + [-1, 1, 1], + [-1, 1, -1], + [-1, -1, 1], + [-1, -1, -1], + ] + ) + + edges = Vector3d( + [ + [0, 1, 1], + [0, 1, -1], + [0, -1, 1], + [0, -1, -1], + [1, 0, 1], + [1, 0, -1], + [-1, 0, 1], + [-1, 0, -1], + [1, 1, 0], + [1, -1, 0], + [-1, 1, 0], + [-1, -1, 0], + ] + ) + + c_index, c_coordinates = _cube_gnom_coordinates(corners) + e_index, e_coordinates = _cube_gnom_coordinates(edges) + + assert np.all( + c_index + == [ + 0, + 5, + 0, + 3, + 2, + 1, + 4, + 1, + ] + ) + assert np.all( + e_index + == [ + 2, + 5, + 4, + 3, + 0, + 5, + 4, + 1, + 0, + 3, + 2, + 1, + ] + ) + + def test_grid_correct_mapping(self): + """Make sure grids get assigned to the correct faces and coordinates.""" + faces_grid = sample_S2_equiangle_cube_mesh_face_centers(15) + index_faces, corrdinates_faces = _cube_gnom_coordinates(faces_grid) + + exp_coords = np.array( + [-0.65449847, -0.39269908, -0.13089969, 0.13089969, 0.39269908, 0.65449847] + ) + for face_index in range(6): + assert np.all(index_faces[face_index] == face_index) + assert np.allclose( + corrdinates_faces[0, face_index], exp_coords[:, np.newaxis] + ) + assert np.allclose( + corrdinates_faces[1, face_index], exp_coords[np.newaxis, :] + ) + + def test_blurring_kernel(self): + """Check that the smoothing gives us roughly the correct width.""" + vectors = Vector3d(np.array([0.0, 0, 1])) + for resolution in [0.25, 0.5, 1.0]: + for s in [5, 10, 20]: + sigma = resolution * s + hist, _ = pole_density_function( + vectors, + sigma=sigma, + resolution=resolution, + ) + assert hist[0, 0] / hist[0, s] > 2.3 + assert hist[0, 0] / hist[0, s] < 3.0