From f4c4b13ec93c3f46e608743fa22f991d6abbe892 Mon Sep 17 00:00:00 2001 From: MACarlsen Date: Fri, 19 Sep 2025 09:24:00 +0200 Subject: [PATCH 01/12] Perform histogramming for density-plots in equi-angle cubed-sphere coordinates to avoid unpleaseant artifacts from to highly isotropic bins in current version. Signed-off-by: MACarlsen --- orix/sampling/S2_sampling.py | 90 ++++++++++++++++++++++++++++++++++++ 1 file changed, 90 insertions(+) diff --git a/orix/sampling/S2_sampling.py b/orix/sampling/S2_sampling.py index a0c24518..9baaeaa4 100644 --- a/orix/sampling/S2_sampling.py +++ b/orix/sampling/S2_sampling.py @@ -392,6 +392,96 @@ 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*. From a1ca645808b1ee4584fe5961cf04b0aff6472943 Mon Sep 17 00:00:00 2001 From: MACarlsen Date: Fri, 19 Sep 2025 11:02:40 +0200 Subject: [PATCH 02/12] subroutines to assign cube-coordinates to vectors and generate bincenters and vertexes --- orix/measure/pole_density_function.py | 107 ++++++++++++++++++++++++ orix/tests/test_measure/test_measure.py | 52 ++++++++++++ 2 files changed, 159 insertions(+) diff --git a/orix/measure/pole_density_function.py b/orix/measure/pole_density_function.py index d9613660..2dc3cb55 100644 --- a/orix/measure/pole_density_function.py +++ b/orix/measure/pole_density_function.py @@ -222,3 +222,110 @@ def pole_density_function( 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 assignes to Index 0 and Index 1 respectively. + + Then each point is given 2D coordinates on the respective + faces. Coordinates are in angles wrt. cube edges. + """ + + if np.any(vectors.norm == 0.0): + #TODO Ask the maintainers what they normally do here. + raise ZeroDivisionError + + # 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 are sure that the denominator is non-zero + # 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 diff --git a/orix/tests/test_measure/test_measure.py b/orix/tests/test_measure/test_measure.py index 0145de1e..f7cf83b3 100644 --- a/orix/tests/test_measure/test_measure.py +++ b/orix/tests/test_measure/test_measure.py @@ -25,6 +25,8 @@ from orix.measure import pole_density_function from orix.quaternion import symmetry from orix.vector import Vector3d +from orix.measure.pole_density_function import _cube_gnom_coordinates +from orix.sampling.S2_sampling import sample_S2_equiangle_cube_mesh_face_centers @pytest.fixture( @@ -158,3 +160,53 @@ def test_pole_density_function_empty_vector_raises(self): ValueError, match="`azimuth` and `polar` angles have 0 size" ): pole_density_function(v) + + +class TestGnomCubeRouties: + + def test_corner_edge_assignment(): + """ Make sire we get useable result 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(): + """ Make sure grids get assigned 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, :]) From b560f153b9bdb64aca285301e7696b9fbcbe1cdf Mon Sep 17 00:00:00 2001 From: MACarlsen Date: Fri, 19 Sep 2025 14:07:41 +0200 Subject: [PATCH 03/12] First working version --- orix/measure/pole_density_function.py | 155 +++++++++++++++++++++----- orix/plot/inverse_pole_figure_plot.py | 2 +- 2 files changed, 126 insertions(+), 31 deletions(-) diff --git a/orix/measure/pole_density_function.py b/orix/measure/pole_density_function.py index 2dc3cb55..eeb1e6d4 100644 --- a/orix/measure/pole_density_function.py +++ b/orix/measure/pole_density_function.py @@ -23,7 +23,7 @@ from orix.projections.stereographic import StereographicProjection from orix.quaternion.symmetry import Symmetry from orix.vector.vector3d import Vector3d - +from scipy.interpolate import RegularGridInterpolator def pole_density_function( *args: np.ndarray | Vector3d, @@ -88,7 +88,6 @@ def 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]) @@ -114,36 +113,89 @@ def pole_density_function( # initial sampling is performed at half the angular resolution resolution /= 2 - azimuth, polar, _ = v.to_polar() + # azimuth, polar, _ = v.to_polar() + face_index_array, face_coordinates = _cube_gnom_coordinates(v) # 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.") + face_index_array = face_index_array.ravel() + face_coordinate_1 = face_coordinates[0].ravel() + face_coordinate_2 = face_coordinates[1].ravel() + + bins = int(90 / resolution) + bin_edges = np.linspace(-np.pi/4, np.pi/4, bins+1) + hist = np.zeros((6, bins, bins)) + for face_index in range(6): + this_face = face_index_array == face_index + hist[face_index], _ = np.histogramdd( + (face_coordinate_1[this_face], face_coordinate_2[this_face]), + bins=(bin_edges, bin_edges), + ) - # Generate equal area mesh on S2 + # 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) + solid_angle_term *= v.size + hist = hist / solid_angle_term[np.newaxis, ...] + + # TODO: If the plot resolution is very high, and the smoothing kernel is very broad + # this blurring has bad performance. In those cases, overwrite the users choice + # for resolution, and then upsample in the end. + + if sigma/resolution > 20 and resolution<1.0: + print('Performance is bad when smoothing-kernel is much larger than\ + the resoltion. Consider using a lower resolution.') + + if sigma != 0.0: + # If smoothing is only a bit, du 60 small steps, otherwise do a max step-size + if (sigma / resolution)**2 <= 20: + N = 60 + t = 1 / N * (sigma / resolution)**2 + else: + t = 1 / 3 + N = int(1 / t * (sigma / resolution)**2) + + hist = _smooth_gnom_cube_histograms(hist, t, N) + + # TODO For now, avoid touching the plotting code, by returning the + # expected format. I will have a deep dive in the plotting later. azimuth_coords, polar_coords = _sample_S2_equal_area_coordinates( - resolution, - hemisphere=hemisphere, + 1, + hemisphere='upper', azimuth_endpoint=True, ) + azimuth_coords 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, + + 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 - # "wrap" along azimuthal axis, "reflect" along polar axis - mode = ("wrap", "reflect") - # apply broadening in angular space - hist = gaussian_filter(hist, sigma / resolution, mode=mode) + grid_face_index, grid_face_coords = _cube_gnom_coordinates(v_center_grid) + hist_grid = np.zeros(v_center_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) + this_face = grid_face_index == face_index + hist_grid[this_face] = interpolator(grid_face_coords[:, this_face].T) + + hist = hist_grid + + # TODO: I don't understand how symmetrization is done. Could just sum over + # rotated histograms instead. # 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( @@ -201,6 +253,10 @@ def pole_density_function( # calculate bin vertices x, y = sp.vector2xy(v_res2_grid) x, y = x.reshape(v_res2_grid.shape), y.reshape(v_res2_grid.shape) + + # This was missing before + hist = hist / symmetry.size + else: # all points valid in stereographic projection hist = np.ma.array(hist, mask=np.zeros_like(hist, dtype=bool)) @@ -209,14 +265,6 @@ def pole_density_function( 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() - if log: # +1 to avoid taking the log of 0 hist = np.log(hist + 1) @@ -235,10 +283,10 @@ def _cube_gnom_coordinates(vectors: Vector3d) -> tuple[np.ndarray[int], np.ndarr 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 assignes to Index 0 and Index 1 respectively. + 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 edges. + faces. Coordinates are in angles wrt. cube face center-lines. """ if np.any(vectors.norm == 0.0): @@ -329,3 +377,50 @@ def _cube_gnom_coordinates(vectors: Vector3d) -> tuple[np.ndarray[int], np.ndarr np.arctan(-unit_vectors.y[this_face] / unit_vectors.z[this_face]) return face_index, coordinates + + +def _smooth_gnom_cube_histograms(histograms, step_parameter, iterations=1): + """ Histograms shape is (6, n_nbins, n_bins) and edge connectivity + is as according to the rest of this file. + """ + sub_histogram_shape = histograms[0][0].shape + output_histogram = np.copy(histograms) + diffused_weight = np.zeros(histograms.shape) + + for n in range(iterations): + + diffused_weight[...] = 0 + + # Diffuse on faces + for face_index in range(6): + + diffused_weight[face_index, 1:, :] += output_histogram[face_index, :-1, :] + diffused_weight[face_index, :-1, :] += output_histogram[face_index, 1:, :] + diffused_weight[face_index, :, 1:] += output_histogram[face_index, :, :-1] + diffused_weight[face_index, :, :-1] += output_histogram[face_index, :, 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/plot/inverse_pole_figure_plot.py b/orix/plot/inverse_pole_figure_plot.py index be01a729..7ceb8cc4 100644 --- a/orix/plot/inverse_pole_figure_plot.py +++ b/orix/plot/inverse_pole_figure_plot.py @@ -109,7 +109,7 @@ def _edge_patch(self) -> mpatches.Patch: def pole_density_function( self, *args: np.ndarray | Vector3d, - resolution: float = 0.25, + resolution: float = 1.0, sigma: float = 5, log: bool = False, colorbar: bool = True, From dc826971724afb3e7db73184435f3bc185ed2c06 Mon Sep 17 00:00:00 2001 From: MACarlsen Date: Fri, 19 Sep 2025 16:16:37 +0200 Subject: [PATCH 04/12] better symmetrization --- orix/measure/pole_density_function.py | 184 +++++++++----------------- orix/plot/inverse_pole_figure_plot.py | 2 +- 2 files changed, 66 insertions(+), 120 deletions(-) diff --git a/orix/measure/pole_density_function.py b/orix/measure/pole_density_function.py index eeb1e6d4..eb4b2da8 100644 --- a/orix/measure/pole_density_function.py +++ b/orix/measure/pole_density_function.py @@ -24,6 +24,7 @@ from orix.quaternion.symmetry import Symmetry from orix.vector.vector3d import Vector3d from scipy.interpolate import RegularGridInterpolator +from orix.sampling.S2_sampling import _sample_S2_uv_mesh_coordinates def pole_density_function( *args: np.ndarray | Vector3d, @@ -85,7 +86,6 @@ 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} @@ -106,29 +106,49 @@ def pole_density_function( "Accepts only one (Vector3d) or two (azimuth, polar) input arguments." ) - 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() - face_index_array, face_coordinates = _cube_gnom_coordinates(v) - # np.histogram2d expects 1d arrays - face_index_array = face_index_array.ravel() - face_coordinate_1 = face_coordinates[0].ravel() - face_coordinate_2 = face_coordinates[1].ravel() - - bins = int(90 / resolution) + # 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)) - for face_index in range(6): - this_face = face_index_array == face_index - hist[face_index], _ = np.histogramdd( - (face_coordinate_1[this_face], face_coordinate_2[this_face]), - bins=(bin_edges, bin_edges), - ) + + if symmetry is None: + face_index_array, face_coordinates = _cube_gnom_coordinates(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 + hist[face_index] = np.histogramdd( + (face_coordinate_1[this_face], face_coordinate_2[this_face]), + bins=(bin_edges, bin_edges), + )[0] + + # Explicit symmetrization + elif symmetry is not None: + 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 + hist[face_index] += np.histogramdd( + (face_coordinate_1[this_face], face_coordinate_2[this_face]), + bins=(bin_edges, bin_edges), + )[0] / symmetry.size + + # Bins are not all same solid angle area, so we need to normalize. if mrd: @@ -140,130 +160,56 @@ def pole_density_function( solid_angle_term *= v.size hist = hist / solid_angle_term[np.newaxis, ...] - # TODO: If the plot resolution is very high, and the smoothing kernel is very broad - # this blurring has bad performance. In those cases, overwrite the users choice - # for resolution, and then upsample in the end. - - if sigma/resolution > 20 and resolution<1.0: - print('Performance is bad when smoothing-kernel is much larger than\ - the resoltion. Consider using a lower resolution.') - + # Smoothing if sigma != 0.0: # If smoothing is only a bit, du 60 small steps, otherwise do a max step-size - if (sigma / resolution)**2 <= 20: + if (sigma / histogram_resolution)**2 <= 20: N = 60 - t = 1 / N * (sigma / resolution)**2 + t = 1 / N * (sigma / histogram_resolution)**2 else: t = 1 / 3 - N = int(1 / t * (sigma / resolution)**2) + N = int(1 / t * (sigma / histogram_resolution)**2) hist = _smooth_gnom_cube_histograms(hist, t, N) - # TODO For now, avoid touching the plotting code, by returning the - # expected format. I will have a deep dive in the plotting later. - azimuth_coords, polar_coords = _sample_S2_equal_area_coordinates( - 1, + + + # Make plot grid + azimuth_coords, polar_coords = _sample_S2_uv_mesh_coordinates( + plot_resolution, hemisphere='upper', azimuth_endpoint=True, ) - azimuth_coords 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", ) - v_center_grid = Vector3d.from_polar( + + 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 - grid_face_index, grid_face_coords = _cube_gnom_coordinates(v_center_grid) - hist_grid = np.zeros(v_center_grid.shape) - bin_middles = np.linspace(-np.pi/4 + np.pi/4/bins, np.pi/4 - np.pi/4/bins, bins) + if symmetry is not None: + + v_grid = v_grid.in_fundamental_sector(symmetry) + v_grid_vertexes = v_grid_vertexes.in_fundamental_sector(symmetry) + # 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) this_face = grid_face_index == face_index hist_grid[this_face] = interpolator(grid_face_coords[:, this_face].T) - hist = hist_grid - # TODO: I don't understand how symmetrization is done. Could just sum over - # rotated histograms instead. - - # 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 - ) - - # 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 - - # compute histogram data array as masked array - hist = np.ma.array( - temp, mask=~(v_center_res2_grid <= symmetry.fundamental_sector) - ) - # calculate bin vertices - x, y = sp.vector2xy(v_res2_grid) - x, y = x.reshape(v_res2_grid.shape), y.reshape(v_res2_grid.shape) - - # This was missing before - hist = hist / symmetry.size - - 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) + # 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 diff --git a/orix/plot/inverse_pole_figure_plot.py b/orix/plot/inverse_pole_figure_plot.py index 7ceb8cc4..be01a729 100644 --- a/orix/plot/inverse_pole_figure_plot.py +++ b/orix/plot/inverse_pole_figure_plot.py @@ -109,7 +109,7 @@ def _edge_patch(self) -> mpatches.Patch: def pole_density_function( self, *args: np.ndarray | Vector3d, - resolution: float = 1.0, + resolution: float = 0.25, sigma: float = 5, log: bool = False, colorbar: bool = True, From b3093ec7e22377b6f9004be71a846c087587cd3d Mon Sep 17 00:00:00 2001 From: MACarlsen Date: Sat, 20 Sep 2025 08:21:17 +0200 Subject: [PATCH 05/12] fixed most tests --- orix/measure/pole_density_function.py | 41 +++++++---- orix/tests/test_measure/test_measure.py | 91 ++++++++++++------------- 2 files changed, 71 insertions(+), 61 deletions(-) diff --git a/orix/measure/pole_density_function.py b/orix/measure/pole_density_function.py index eb4b2da8..2ba2d8af 100644 --- a/orix/measure/pole_density_function.py +++ b/orix/measure/pole_density_function.py @@ -36,11 +36,12 @@ def pole_density_function( 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 ---------- @@ -106,6 +107,12 @@ def pole_density_function( "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: @@ -128,9 +135,11 @@ def pole_density_function( 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] # Explicit symmetrization @@ -143,13 +152,13 @@ def pole_density_function( 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) @@ -157,7 +166,10 @@ def pole_density_function( 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) - solid_angle_term *= v.size + if weights is not None: + solid_angle_term *= np.sum(weights) + else: + solid_angle_term *= v.size hist = hist / solid_angle_term[np.newaxis, ...] # Smoothing @@ -172,8 +184,6 @@ def pole_density_function( hist = _smooth_gnom_cube_histograms(hist, t, N) - - # Make plot grid azimuth_coords, polar_coords = _sample_S2_uv_mesh_coordinates( plot_resolution, @@ -193,7 +203,7 @@ def pole_density_function( v_grid_vertexes = Vector3d.from_polar(azimuth=azimuth_grid, polar=polar_grid).unit if symmetry is not None: - + mask=~(v_grid <= symmetry.fundamental_sector) v_grid = v_grid.in_fundamental_sector(symmetry) v_grid_vertexes = v_grid_vertexes.in_fundamental_sector(symmetry) @@ -207,6 +217,12 @@ def pole_density_function( hist_grid[this_face] = interpolator(grid_face_coords[:, this_face].T) hist = hist_grid + # Mask out points outside funamental region if symmetry is given. + if symmetry is not None: + hist = np.ma.array(hist, mask=mask) + else: + hist = np.ma.array(hist, mask=np.zeros_like(hist, dtype=bool)) + # 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) @@ -233,11 +249,12 @@ def _cube_gnom_coordinates(vectors: Vector3d) -> tuple[np.ndarray[int], np.ndarr 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): - #TODO Ask the maintainers what they normally do here. - raise ZeroDivisionError + raise ValueError # Assign face index to each vector face_index = np.zeros(vectors.shape, dtype=int) diff --git a/orix/tests/test_measure/test_measure.py b/orix/tests/test_measure/test_measure.py index f7cf83b3..f989eac7 100644 --- a/orix/tests/test_measure/test_measure.py +++ b/orix/tests/test_measure/test_measure.py @@ -27,7 +27,7 @@ from orix.vector import Vector3d from orix.measure.pole_density_function import _cube_gnom_coordinates from orix.sampling.S2_sampling import sample_S2_equiangle_cube_mesh_face_centers - +from orix.sampling.S2_sampling import _sample_S2_uv_mesh_coordinates, sample_S2_random_mesh @pytest.fixture( params=[ @@ -47,13 +47,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 @@ -66,43 +65,25 @@ def test_pole_density_function(self): assert isinstance(hist2, np.ma.MaskedArray) assert hist2.mask.sum() > 0 - assert np.allclose(hist1.mean(), hist2.mean()) - - 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) @@ -132,10 +113,11 @@ 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 @@ -157,15 +139,15 @@ def test_pole_density_function_empty_vector_raises(self): assert not v.size with pytest.raises( - ValueError, match="`azimuth` and `polar` angles have 0 size" + ValueError ): pole_density_function(v) -class TestGnomCubeRouties: +class TestGnomCubeRoutines: - def test_corner_edge_assignment(): - """ Make sire we get useable result for corner-cases. + def test_corner_edge_assignment(self): + """ Make sure we get useable results for corner-cases. """ corners = Vector3d( [[1, 1, 1], @@ -199,8 +181,8 @@ def test_corner_edge_assignment(): 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(): - """ Make sure grids get assigned correct faces and coordinates. + 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) @@ -210,3 +192,14 @@ def test_grid_correct_mapping(): 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 From 4eaa9ed39dcbb6e404653d5d4576d2f26bbf6cfb Mon Sep 17 00:00:00 2001 From: MACarlsen Date: Sat, 20 Sep 2025 08:21:57 +0200 Subject: [PATCH 06/12] fixed seed for tests --- orix/tests/test_measure/test_measure.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/orix/tests/test_measure/test_measure.py b/orix/tests/test_measure/test_measure.py index f989eac7..197831f3 100644 --- a/orix/tests/test_measure/test_measure.py +++ b/orix/tests/test_measure/test_measure.py @@ -86,21 +86,24 @@ def test_pole_density_function_mrd_norm(self, point_groups, resolution): 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) From 199ff64ab51898cf1dad14c5a2044f31b8d038f4 Mon Sep 17 00:00:00 2001 From: MACarlsen Date: Sat, 20 Sep 2025 08:49:58 +0200 Subject: [PATCH 07/12] flake8 warnings and type hints --- orix/measure/pole_density_function.py | 86 +++++++++++++++---------- orix/tests/test_measure/test_measure.py | 81 ++++++++++++++--------- 2 files changed, 102 insertions(+), 65 deletions(-) diff --git a/orix/measure/pole_density_function.py b/orix/measure/pole_density_function.py index 2ba2d8af..54d4e608 100644 --- a/orix/measure/pole_density_function.py +++ b/orix/measure/pole_density_function.py @@ -18,7 +18,6 @@ # import numpy as np -from scipy.ndimage import gaussian_filter from orix.projections.stereographic import StereographicProjection from orix.quaternion.symmetry import Symmetry @@ -26,6 +25,7 @@ from scipy.interpolate import RegularGridInterpolator from orix.sampling.S2_sampling import _sample_S2_uv_mesh_coordinates + def pole_density_function( *args: np.ndarray | Vector3d, resolution: float = 1, @@ -104,7 +104,8 @@ def pole_density_function( 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: @@ -137,7 +138,8 @@ def pole_density_function( 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]), + (face_coordinate_1[this_face], + face_coordinate_2[this_face]), bins=(bin_edges, bin_edges), weights=w, )[0] @@ -146,7 +148,8 @@ def pole_density_function( elif symmetry is not None: for rotation in symmetry: - face_index_array, face_coordinates = _cube_gnom_coordinates(rotation*v) + 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() @@ -154,17 +157,20 @@ def pole_density_function( 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]), + (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) + 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 / (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) @@ -174,7 +180,8 @@ def pole_density_function( # Smoothing if sigma != 0.0: - # If smoothing is only a bit, du 60 small steps, otherwise do a max step-size + # 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 @@ -190,7 +197,11 @@ def pole_density_function( hemisphere='upper', azimuth_endpoint=True, ) - azimuth_grid, polar_grid = np.meshgrid(azimuth_coords, polar_coords, indexing="ij") + 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, @@ -200,19 +211,28 @@ def pole_density_function( 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 + v_grid_vertexes = Vector3d.from_polar( + azimuth=azimuth_grid, + polar=polar_grid + ).unit if symmetry is not None: - mask=~(v_grid <= symmetry.fundamental_sector) + mask = ~(v_grid <= symmetry.fundamental_sector) v_grid = v_grid.in_fundamental_sector(symmetry) v_grid_vertexes = v_grid_vertexes.in_fundamental_sector(symmetry) # 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) + 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) + interpolator = RegularGridInterpolator( + (bin_middles, bin_middles), + hist[face_index], + bounds_error=False, + fill_value=None + ) this_face = grid_face_index == face_index hist_grid[this_face] = interpolator(grid_face_coords[:, this_face].T) hist = hist_grid @@ -234,7 +254,8 @@ def pole_density_function( return hist, (x, y) -def _cube_gnom_coordinates(vectors: Vector3d) -> tuple[np.ndarray[int], np.ndarray[float]]: +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: @@ -274,35 +295,33 @@ def _cube_gnom_coordinates(vectors: Vector3d) -> tuple[np.ndarray[int], np.ndarr indx = np.all([vectors.y > vectors.x, vectors.y >= -vectors.x, vectors.y >= vectors.z, - vectors.y > -vectors.z,], axis = 0) + 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) + 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) + 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) + 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 are sure that the denominator is non-zero + # 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]) @@ -342,11 +361,14 @@ def _cube_gnom_coordinates(vectors: Vector3d) -> tuple[np.ndarray[int], np.ndarr return face_index, coordinates -def _smooth_gnom_cube_histograms(histograms, step_parameter, iterations=1): +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. """ - sub_histogram_shape = histograms[0][0].shape output_histogram = np.copy(histograms) diffused_weight = np.zeros(histograms.shape) @@ -355,13 +377,11 @@ def _smooth_gnom_cube_histograms(histograms, step_parameter, iterations=1): diffused_weight[...] = 0 # Diffuse on faces - for face_index in range(6): - - diffused_weight[face_index, 1:, :] += output_histogram[face_index, :-1, :] - diffused_weight[face_index, :-1, :] += output_histogram[face_index, 1:, :] - diffused_weight[face_index, :, 1:] += output_histogram[face_index, :, :-1] - diffused_weight[face_index, :, :-1] += output_histogram[face_index, :, 1:] - + 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 @@ -379,8 +399,8 @@ def _smooth_gnom_cube_histograms(histograms, step_parameter, iterations=1): ) 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] + 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\ diff --git a/orix/tests/test_measure/test_measure.py b/orix/tests/test_measure/test_measure.py index 197831f3..f87eea05 100644 --- a/orix/tests/test_measure/test_measure.py +++ b/orix/tests/test_measure/test_measure.py @@ -17,8 +17,6 @@ # along with orix. If not, see . # -from copy import deepcopy - import numpy as np import pytest @@ -26,8 +24,12 @@ from orix.quaternion import symmetry from orix.vector import Vector3d from orix.measure.pole_density_function import _cube_gnom_coordinates -from orix.sampling.S2_sampling import sample_S2_equiangle_cube_mesh_face_centers -from orix.sampling.S2_sampling import _sample_S2_uv_mesh_coordinates, sample_S2_random_mesh +from orix.sampling.S2_sampling import ( + sample_S2_equiangle_cube_mesh_face_centers, + _sample_S2_uv_mesh_coordinates, + sample_S2_random_mesh, + ) + @pytest.fixture( params=[ @@ -80,8 +82,14 @@ def test_pole_density_function_mrd_norm(self, point_groups, resolution): 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) + 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) @@ -116,9 +124,10 @@ def test_pole_density_function_weights(self): # the same because MRD normalizes by average assert np.allclose(hist0, hist2) - #TDOD: Ask about expected behaviour for 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) + # 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) @@ -131,7 +140,8 @@ def test_PDF_IPDF_equivalence(self): v = Vector3d.random(100_000) hist_pdf, _ = pole_density_function(v, weights=None) - hist_ipdf, _ = pole_density_function(v, weights=None, symmetry=symmetry.C1) + hist_ipdf, _ = pole_density_function(v, weights=None, + symmetry=symmetry.C1) # in testing this test passes at tolerance of 1% for 100_000 # vectors, but raise tolerance to 2% to ensure pass @@ -154,28 +164,28 @@ def test_corner_edge_assignment(self): """ 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],] + [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],] + [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) @@ -190,19 +200,26 @@ def test_grid_correct_mapping(self): 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]) + 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, :]) + 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])) + 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) + 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 From 2913226a08f8d07c839a43557221159d0b289b8c Mon Sep 17 00:00:00 2001 From: MACarlsen Date: Sat, 20 Sep 2025 08:57:41 +0200 Subject: [PATCH 08/12] simplyfy code by different handling of no symmetry case --- orix/measure/pole_density_function.py | 61 +++++++++---------------- orix/tests/test_measure/test_measure.py | 1 + 2 files changed, 23 insertions(+), 39 deletions(-) diff --git a/orix/measure/pole_density_function.py b/orix/measure/pole_density_function.py index 54d4e608..679e8ef1 100644 --- a/orix/measure/pole_density_function.py +++ b/orix/measure/pole_density_function.py @@ -20,7 +20,7 @@ import numpy as np from orix.projections.stereographic import StereographicProjection -from orix.quaternion.symmetry import Symmetry +from orix.quaternion.symmetry import Symmetry, C1 from orix.vector.vector3d import Vector3d from scipy.interpolate import RegularGridInterpolator from orix.sampling.S2_sampling import _sample_S2_uv_mesh_coordinates @@ -32,7 +32,7 @@ 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]]: @@ -92,6 +92,10 @@ def pole_density_function( 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): @@ -128,40 +132,23 @@ def pole_density_function( bin_edges = np.linspace(-np.pi/4, np.pi/4, bins+1) hist = np.zeros((6, bins, bins)) - if symmetry is None: - face_index_array, face_coordinates = _cube_gnom_coordinates(v) + # 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] - - # Explicit symmetrization - elif symmetry is not None: - 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 + 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: @@ -216,10 +203,9 @@ def pole_density_function( polar=polar_grid ).unit - if symmetry is not None: - mask = ~(v_grid <= symmetry.fundamental_sector) - v_grid = v_grid.in_fundamental_sector(symmetry) - v_grid_vertexes = v_grid_vertexes.in_fundamental_sector(symmetry) + mask = ~(v_grid <= symmetry.fundamental_sector) + v_grid = v_grid.in_fundamental_sector(symmetry) + v_grid_vertexes = v_grid_vertexes.in_fundamental_sector(symmetry) # Interpolation from histograms to plot grid grid_face_index, grid_face_coords = _cube_gnom_coordinates(v_grid) @@ -237,11 +223,8 @@ def pole_density_function( hist_grid[this_face] = interpolator(grid_face_coords[:, this_face].T) hist = hist_grid - # Mask out points outside funamental region if symmetry is given. - if symmetry is not None: - hist = np.ma.array(hist, mask=mask) - else: - hist = np.ma.array(hist, mask=np.zeros_like(hist, dtype=bool)) + # Mask out points outside funamental region. + hist = np.ma.array(hist, mask=mask) # Transform grdi to mystery coordinates used by plotting routine x, y = sp.vector2xy(v_grid_vertexes) diff --git a/orix/tests/test_measure/test_measure.py b/orix/tests/test_measure/test_measure.py index f87eea05..7ceab8bb 100644 --- a/orix/tests/test_measure/test_measure.py +++ b/orix/tests/test_measure/test_measure.py @@ -33,6 +33,7 @@ @pytest.fixture( params=[ + None, symmetry.D2h, symmetry.S6, symmetry.D3d, From c90e5fe336eeb706ebdb36c0a3c145b2e1bc02b4 Mon Sep 17 00:00:00 2001 From: MACarlsen Date: Sat, 20 Sep 2025 09:01:34 +0200 Subject: [PATCH 09/12] black --- orix/measure/pole_density_function.py | 209 +++++++++++++++--------- orix/sampling/S2_sampling.py | 14 +- orix/tests/test_measure/test_measure.py | 124 ++++++++------ 3 files changed, 212 insertions(+), 135 deletions(-) diff --git a/orix/measure/pole_density_function.py b/orix/measure/pole_density_function.py index 679e8ef1..80b72ee1 100644 --- a/orix/measure/pole_density_function.py +++ b/orix/measure/pole_density_function.py @@ -129,35 +129,40 @@ def pole_density_function( # Do actual histogramming bins = int(90 / histogram_resolution) - bin_edges = np.linspace(-np.pi/4, np.pi/4, bins+1) + 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_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 + 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) + 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 + / (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) @@ -169,26 +174,26 @@ def pole_density_function( 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: + if (sigma / histogram_resolution) ** 2 <= 20: N = 60 - t = 1 / N * (sigma / histogram_resolution)**2 + t = 1 / N * (sigma / histogram_resolution) ** 2 else: t = 1 / 3 - N = int(1 / t * (sigma / histogram_resolution)**2) + 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='upper', + hemisphere="upper", azimuth_endpoint=True, ) 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, @@ -198,10 +203,7 @@ def pole_density_function( 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 + v_grid_vertexes = Vector3d.from_polar(azimuth=azimuth_grid, polar=polar_grid).unit mask = ~(v_grid <= symmetry.fundamental_sector) v_grid = v_grid.in_fundamental_sector(symmetry) @@ -210,15 +212,16 @@ def pole_density_function( # 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) + 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 - ) + fill_value=None, + ) this_face = grid_face_index == face_index hist_grid[this_face] = interpolator(grid_face_coords[:, this_face].T) hist = hist_grid @@ -237,9 +240,10 @@ def pole_density_function( 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``` +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.) @@ -263,40 +267,70 @@ def _cube_gnom_coordinates(vectors: Vector3d)\ # 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) + 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) + 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) + 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) + 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) + 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) + 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 @@ -306,40 +340,52 @@ def _cube_gnom_coordinates(vectors: Vector3d)\ # 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]) + 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]) + 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]) + 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]) + 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]) + 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]) + 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 @@ -349,7 +395,7 @@ def _smooth_gnom_cube_histograms( step_parameter: float, iterations: int = 1, ) -> np.ndarray[float]: - """ Histograms shape is (6, n_nbins, n_bins) and edge connectivity + """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) @@ -386,7 +432,8 @@ def _smooth_gnom_cube_histograms( diffused_weight[edge_2] += output_histogram[edge_1] # Add to output - output_histogram = (1-step_parameter)*output_histogram\ - + diffused_weight/4*step_parameter + 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 9baaeaa4..7f5c4700 100644 --- a/orix/sampling/S2_sampling.py +++ b/orix/sampling/S2_sampling.py @@ -417,15 +417,15 @@ def sample_S2_equiangle_cube_mesh_vertexes( """ bins_along_edge = int(90 / resolution) - bin_edges = np.linspace(-np.pi/4, np.pi/4, bins_along_edge+1) + 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)) + 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 + 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)) + 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) @@ -462,8 +462,10 @@ def sample_S2_equiangle_cube_mesh_face_centers( """ 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)) + 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 diff --git a/orix/tests/test_measure/test_measure.py b/orix/tests/test_measure/test_measure.py index 7ceab8bb..d0ad9e78 100644 --- a/orix/tests/test_measure/test_measure.py +++ b/orix/tests/test_measure/test_measure.py @@ -28,7 +28,7 @@ sample_S2_equiangle_cube_mesh_face_centers, _sample_S2_uv_mesh_coordinates, sample_S2_random_mesh, - ) +) @pytest.fixture( @@ -76,21 +76,20 @@ def test_pole_density_function_mrd_norm(self, point_groups, resolution): # Make plot grid _, polar_coords = _sample_S2_uv_mesh_coordinates( resolution, - hemisphere='upper', + 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]: + 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) + ) + 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) @@ -141,8 +140,7 @@ def test_PDF_IPDF_equivalence(self): v = Vector3d.random(100_000) hist_pdf, _ = pole_density_function(v, weights=None) - hist_ipdf, _ = pole_density_function(v, weights=None, - symmetry=symmetry.C1) + hist_ipdf, _ = pole_density_function(v, weights=None, symmetry=symmetry.C1) # in testing this test passes at tolerance of 1% for 100_000 # vectors, but raise tolerance to 2% to ensure pass @@ -152,67 +150,97 @@ def test_pole_density_function_empty_vector_raises(self): v = Vector3d.empty() assert not v.size - with pytest.raises( - ValueError - ): + 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. - """ + """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],] - ) + [ + [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],] - ) + [ + [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,]) + 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. - """ + """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]) + 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, :]) + 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. - """ + """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]: @@ -221,6 +249,6 @@ def test_blurring_kernel(self): vectors, sigma=sigma, resolution=resolution, - ) + ) assert hist[0, 0] / hist[0, s] > 2.3 assert hist[0, 0] / hist[0, s] < 3.0 From 9813725bb687900485aa56cdb6f4c3282b486eb0 Mon Sep 17 00:00:00 2001 From: MACarlsen Date: Sat, 20 Sep 2025 09:05:46 +0200 Subject: [PATCH 10/12] isort --- orix/measure/pole_density_function.py | 6 +++--- orix/tests/test_measure/test_measure.py | 6 +++--- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/orix/measure/pole_density_function.py b/orix/measure/pole_density_function.py index 80b72ee1..0e55003e 100644 --- a/orix/measure/pole_density_function.py +++ b/orix/measure/pole_density_function.py @@ -18,12 +18,12 @@ # import numpy as np +from scipy.interpolate import RegularGridInterpolator from orix.projections.stereographic import StereographicProjection -from orix.quaternion.symmetry import Symmetry, C1 -from orix.vector.vector3d import Vector3d -from scipy.interpolate import RegularGridInterpolator +from orix.quaternion.symmetry import C1, Symmetry from orix.sampling.S2_sampling import _sample_S2_uv_mesh_coordinates +from orix.vector.vector3d import Vector3d def pole_density_function( diff --git a/orix/tests/test_measure/test_measure.py b/orix/tests/test_measure/test_measure.py index d0ad9e78..0813fb29 100644 --- a/orix/tests/test_measure/test_measure.py +++ b/orix/tests/test_measure/test_measure.py @@ -21,14 +21,14 @@ import pytest from orix.measure import pole_density_function -from orix.quaternion import symmetry -from orix.vector import Vector3d from orix.measure.pole_density_function import _cube_gnom_coordinates +from orix.quaternion import symmetry from orix.sampling.S2_sampling import ( - sample_S2_equiangle_cube_mesh_face_centers, _sample_S2_uv_mesh_coordinates, + sample_S2_equiangle_cube_mesh_face_centers, sample_S2_random_mesh, ) +from orix.vector import Vector3d @pytest.fixture( From c9afd0860579ce5910e09068918f3e4462a6005a Mon Sep 17 00:00:00 2001 From: MACarlsen Date: Sun, 21 Sep 2025 07:14:47 +0200 Subject: [PATCH 11/12] fix for ``hemisphere='lower'`` --- orix/measure/pole_density_function.py | 3 ++- orix/tests/test_measure/test_measure.py | 6 ++++++ 2 files changed, 8 insertions(+), 1 deletion(-) diff --git a/orix/measure/pole_density_function.py b/orix/measure/pole_density_function.py index 0e55003e..0810d5da 100644 --- a/orix/measure/pole_density_function.py +++ b/orix/measure/pole_density_function.py @@ -106,6 +106,7 @@ 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)\ @@ -186,7 +187,7 @@ def pole_density_function( # Make plot grid azimuth_coords, polar_coords = _sample_S2_uv_mesh_coordinates( plot_resolution, - hemisphere="upper", + hemisphere=hemisphere, azimuth_endpoint=True, ) azimuth_grid, polar_grid = np.meshgrid( diff --git a/orix/tests/test_measure/test_measure.py b/orix/tests/test_measure/test_measure.py index 0813fb29..97c9a37d 100644 --- a/orix/tests/test_measure/test_measure.py +++ b/orix/tests/test_measure/test_measure.py @@ -68,6 +68,12 @@ def test_output_format(self): assert isinstance(hist2, np.ma.MaskedArray) assert hist2.mask.sum() > 0 + 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 + @pytest.mark.parametrize("resolution", [0.5, 1.0]) def test_pole_density_function_mrd_norm(self, point_groups, resolution): pg = point_groups From 9dff947d7d642f2396663150c9eee56323fc09dd Mon Sep 17 00:00:00 2001 From: MACarlsen <44873424+MACarlsen@users.noreply.github.com> Date: Thu, 25 Sep 2025 08:49:34 +0200 Subject: [PATCH 12/12] replace restriction to fundam zone with masking --- orix/measure/pole_density_function.py | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/orix/measure/pole_density_function.py b/orix/measure/pole_density_function.py index 0810d5da..af2aae89 100644 --- a/orix/measure/pole_density_function.py +++ b/orix/measure/pole_density_function.py @@ -206,10 +206,6 @@ def pole_density_function( ).unit v_grid_vertexes = Vector3d.from_polar(azimuth=azimuth_grid, polar=polar_grid).unit - mask = ~(v_grid <= symmetry.fundamental_sector) - v_grid = v_grid.in_fundamental_sector(symmetry) - v_grid_vertexes = v_grid_vertexes.in_fundamental_sector(symmetry) - # Interpolation from histograms to plot grid grid_face_index, grid_face_coords = _cube_gnom_coordinates(v_grid) hist_grid = np.zeros(v_grid.shape) @@ -228,6 +224,7 @@ def pole_density_function( 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