From aaa1c9259ee30fea3bea6c7cfa5941070724132e Mon Sep 17 00:00:00 2001 From: Justus Magin Date: Fri, 27 Sep 2024 13:09:34 +0200 Subject: [PATCH] use the dot-product identity to compute angles between vectors (#47) * use the dot-product identity for angles between vectors * use the new function to compute the angular distances * rename a test * clip to the range of arccos --- healpix_convolution/distances.py | 15 +++++++++++---- healpix_convolution/tests/test_distances.py | 12 +++++++++++- 2 files changed, 22 insertions(+), 5 deletions(-) diff --git a/healpix_convolution/distances.py b/healpix_convolution/distances.py index 4eb016d..e8f654e 100644 --- a/healpix_convolution/distances.py +++ b/healpix_convolution/distances.py @@ -16,6 +16,16 @@ def cell_ids2vectors(cell_ids, nside, nest): return np.reshape(vecs, cell_ids.shape + (3,)) +def angle_between_vectors(a, b, axis): + length_a = np.linalg.norm(a, axis=axis) + length_b = np.linalg.norm(b, axis=axis) + dot_product = np.sum(a * b, axis=axis) + + argument = np.clip(dot_product / (length_a * length_b), -1, 1) + + return np.arccos(argument) + + def _distances(a, b, axis, nside, nest): vec_a = cell_ids2vectors(a, nside, nest) @@ -24,10 +34,7 @@ def _distances(a, b, axis, nside, nest): vec_b_ = cell_ids2vectors(np.where(mask, b, 0), nside, nest) vec_b = np.where(mask[..., None], vec_b_, np.nan) - dot_product = np.abs(np.sum(vec_a * vec_b, axis=axis)) - cross_product = np.linalg.norm(np.cross(vec_a, vec_b, axis=axis), axis=axis) - - return np.arctan2(cross_product, dot_product) + return angle_between_vectors(vec_a, vec_b, axis=axis) def angular_distances(neighbours, *, resolution, indexing_scheme="nested", axis=None): diff --git a/healpix_convolution/tests/test_distances.py b/healpix_convolution/tests/test_distances.py index 1688226..cdfa5e4 100644 --- a/healpix_convolution/tests/test_distances.py +++ b/healpix_convolution/tests/test_distances.py @@ -4,6 +4,16 @@ import healpix_convolution.distances as hds +def test_angle_between_vectors(): + angles = np.linspace(0, 2 * np.pi, 10) + vectors = np.stack([np.cos(angles), np.sin(angles), np.zeros_like(angles)], axis=-1) + + actual = hds.angle_between_vectors(vectors[:1, :], vectors, axis=-1) + expected = np.pi - abs(angles - np.pi) + + np.testing.assert_allclose(actual, expected) + + @pytest.mark.parametrize( ["neighbours", "expected"], ( @@ -11,7 +21,7 @@ (np.array([[2, 71, 8]]), np.array([[0, 0.25637566, 0.25574741]])), ), ) -def test_distances_numpy(neighbours, expected): +def test_angular_distances_numpy(neighbours, expected): actual = hds.angular_distances(neighbours, resolution=2, indexing_scheme="nested") np.testing.assert_allclose(actual, expected)