From 39a0665f637191152605911cf209fc16a36e5ae9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Micha=C3=ABl=20Defferrard?= Date: Mon, 30 Nov 2020 18:04:35 +0100 Subject: [PATCH] SphereIcosahedral & SphereHealpix: allow subdivisions that aren't powers of 2, consistent with SphereCubed --- doc/background/spherical_samplings.rst | 22 ++++++++--------- pygsp/graphs/nngraphs/spherehealpix.py | 28 +++++++++++----------- pygsp/graphs/nngraphs/sphereicosahedral.py | 17 +++++++++---- pygsp/tests/test_graphs.py | 24 +++++++++++-------- 4 files changed, 51 insertions(+), 40 deletions(-) diff --git a/doc/background/spherical_samplings.rst b/doc/background/spherical_samplings.rst index ee5529c1..03f4a6fc 100644 --- a/doc/background/spherical_samplings.rst +++ b/doc/background/spherical_samplings.rst @@ -148,12 +148,12 @@ Subdivision of the icosahedron (:class:`~pygsp.graphs.SphereIcosahedral`) ------------------------------------------------------------------------- The sampling is made of an `icosahedron`_ (a `regular and convex polyhedron`_ -made of 12 vertices and 20 triangular faces) whose faces are recursively -subdivided (into four triangles) to the desired resolution. -The result is a :math:`\{3,5+\}_{2^m,0}` `geodesic polyhedron`_ (made of -:math:`20⋅4^m` triangles, :math:`n=10⋅4^m+2` vertices) or its dual, a -:math:`\{5+,3\}_{2^m,0}` `Goldberg polyhedron`_ (made of :math:`10⋅4^m-10` -hexagons, 12 pentagons, :math:`n=20⋅4^m` vertices). +made of 12 vertices and 20 triangular faces) whose faces are subdivided into +:math:`m^2` triangles projected to the sphere. +The result is a :math:`\{3,5+\}_{m,0}` `geodesic polyhedron`_ (made of +:math:`20⋅m^2` triangles, :math:`n=10⋅m^2+2` vertices) or its dual, a +:math:`\{5+,3\}_{m,0}` `Goldberg polyhedron`_ (made of :math:`10⋅m^2-10` +hexagons, 12 pentagons, :math:`n=20⋅m^2` vertices). The resulting `polyhedral graph`_ (i.e., a 3-vertex-connected planar graph) is the 1-`skeleton`_ of the polyhedron. All have `icosahedral symmetry`_. @@ -176,8 +176,8 @@ SHT. Subdivision of the cube (:class:`~pygsp.graphs.SphereCubed`) ------------------------------------------------------------ -The sampling is made of a cube whose faces are subdivided into finer -quadrilaterals to the desired resolution. +The sampling is made of a cube whose faces are subdivided into :math:`m^2` +finer quadrilaterals projected to the sphere. The result is a convex polyhedron made of :math:`n=6⋅m^2` faces. The graph vertices represent the quadrilateral faces. @@ -190,9 +190,9 @@ Subdivision of the rhombic dodecahedron (:class:`~pygsp.graphs.SphereHealpix`) ------------------------------------------------------------------------------ The sampling is made of a `rhombic dodecahedron`_ (a convex polyhedron made of -12 rhombic faces) whose faces are recursively subdivided (into four -quadrilaterals) to the desired resolution. -The result is a convex polyhedron made of :math:`n=12⋅4^m` faces. +12 rhombic faces) whose faces are subdivided into :math:`m^2` finer +quadrilaterals projected to the sphere. +The result is a convex polyhedron made of :math:`n=12⋅m^2` faces. The graph vertices represent the quadrilateral faces. The Hierarchical Equal Area isoLatitude Pixelisation (`HEALPix`_) [Go]_ was diff --git a/pygsp/graphs/nngraphs/spherehealpix.py b/pygsp/graphs/nngraphs/spherehealpix.py index 039d9d60..46b5fb7e 100644 --- a/pygsp/graphs/nngraphs/spherehealpix.py +++ b/pygsp/graphs/nngraphs/spherehealpix.py @@ -27,8 +27,9 @@ class SphereHealpix(NNGraph): Parameters ---------- subdivisions : int - Number of recursive subdivisions. Known as ``order`` in HEALPix - terminology, with ``nside=2**order`` and ``npix=12*4**order``. + Number of edges the dodecahedron's edges are divided into (``nside``), + resulting in ``12*subdivisions**2`` vertices (``npix``). + It must be a power of 2 if ``nest=True`` (to be hierarchical). indexes : array_like of int Indexes of the pixels from which to build a graph. Useful to build a graph from a subset of the pixels, e.g., for partial sky observations. @@ -66,7 +67,7 @@ class SphereHealpix(NNGraph): Examples -------- >>> import matplotlib.pyplot as plt - >>> G = graphs.SphereHealpix(1, k=8) + >>> G = graphs.SphereHealpix() >>> fig = plt.figure() >>> ax1 = fig.add_subplot(131) >>> ax2 = fig.add_subplot(132, projection='3d') @@ -80,27 +81,26 @@ class SphereHealpix(NNGraph): >>> import matplotlib.pyplot as plt >>> fig, axes = plt.subplots(1, 2) - >>> graph = graphs.SphereHealpix(1, nest=False, k=8) + >>> graph = graphs.SphereHealpix(2, nest=False, k=8) >>> graph.set_coordinates('sphere', dim=2) >>> _ = graph.plot(indices=True, ax=axes[0], title='RING ordering') - >>> graph = graphs.SphereHealpix(1, nest=True, k=8) + >>> graph = graphs.SphereHealpix(2, nest=True, k=8) >>> graph.set_coordinates('sphere', dim=2) >>> _ = graph.plot(indices=True, ax=axes[1], title='NESTED ordering') """ - def __init__(self, subdivisions=1, indexes=None, nest=False, **kwargs): + def __init__(self, subdivisions=2, indexes=None, nest=False, **kwargs): hp = _import_hp() - nside = hp.order2nside(subdivisions) self.subdivisions = subdivisions self.nest = nest if indexes is None: - npix = hp.nside2npix(nside) + npix = hp.nside2npix(subdivisions) indexes = np.arange(npix) - x, y, z = hp.pix2vec(nside, indexes, nest=nest) + x, y, z = hp.pix2vec(subdivisions, indexes, nest=nest) coords = np.stack([x, y, z], axis=1) k = kwargs.pop('k', 20) @@ -108,16 +108,16 @@ def __init__(self, subdivisions=1, indexes=None, nest=False, **kwargs): kernel_width = kwargs.pop('kernel_width') except KeyError: try: - kernel_width = _OPTIMAL_KERNEL_WIDTHS[k][nside] + kernel_width = _OPTIMAL_KERNEL_WIDTHS[k][subdivisions] except KeyError: - raise ValueError('No known optimal kernel width for {} ' - 'neighbors and nside={}.'.format(k, nside)) + raise ValueError('No optimal kernel width for {} neighbors and' + ' {} subdivisions.'.format(k, subdivisions)) super(SphereHealpix, self).__init__(coords, k=k, kernel_width=kernel_width, **kwargs) - lat, lon = hp.pix2ang(nside, indexes, nest=nest, lonlat=False) + lat, lon = hp.pix2ang(subdivisions, indexes, nest=nest, lonlat=False) self.signals['lat'] = np.pi/2 - lat # colatitude to latitude self.signals['lon'] = lon @@ -130,7 +130,7 @@ def _get_extra_repr(self): return attrs -# TODO: find an interpolation between nside and k (#neighbors). +# TODO: find an interpolation between subdivisions (nside) and k (#neighbors). _OPTIMAL_KERNEL_WIDTHS = { 8: { 1: 0.02500 * 32, # extrapolated diff --git a/pygsp/graphs/nngraphs/sphereicosahedral.py b/pygsp/graphs/nngraphs/sphereicosahedral.py index 64b3715a..17747aa8 100644 --- a/pygsp/graphs/nngraphs/sphereicosahedral.py +++ b/pygsp/graphs/nngraphs/sphereicosahedral.py @@ -26,7 +26,10 @@ class SphereIcosahedral(NNGraph): Parameters ---------- subdivisions : int - Number of recursive subdivisions. + Number of edges the icosahedron's edges are divided into, resulting in + ``10*subdivisions**2+2`` vertices (or ``20*subdivisions**2`` if + ``dual=True``). + It must be a power of 2 in the current implementation. dual : bool Whether the graph vertices correspond to the vertices (``dual=False``) or the triangular faces (``dual=True``) of the subdivided icosahedron. @@ -77,15 +80,15 @@ class SphereIcosahedral(NNGraph): >>> import matplotlib.pyplot as plt >>> fig, axes = plt.subplots(1, 2) - >>> graph = graphs.SphereIcosahedral(0, dual=False, k=5) + >>> graph = graphs.SphereIcosahedral(1, dual=False, k=5) >>> graph.set_coordinates('sphere', dim=2) >>> _ = graph.plot(indices=True, ax=axes[0], title='Icosahedron') - >>> graph = graphs.SphereIcosahedral(0, dual=True, k=3) + >>> graph = graphs.SphereIcosahedral(1, dual=True, k=3) >>> graph.set_coordinates('sphere', dim=2) >>> _ = graph.plot(indices=True, ax=axes[1], title='Dodecahedron') """ - def __init__(self, subdivisions=1, dual=False, **kwargs): + def __init__(self, subdivisions=2, dual=False, **kwargs): self.subdivisions = subdivisions self.dual = dual @@ -115,7 +118,11 @@ def normalize(vertices): """Project the vertices on the sphere.""" vertices /= np.linalg.norm(vertices, axis=1)[:, None] - for _ in range(subdivisions): + if np.log2(subdivisions) % 1 != 0: + raise NotImplementedError('Only recursive subdivisions by two are ' + 'implemented. Choose a power of two.') + + for _ in range(int(np.log2(subdivisions))): mesh = mesh.subdivide() # TODO: shall we project between subdivisions? Some do, some don't. # Projecting pushes points away from the 12 base vertices, which diff --git a/pygsp/tests/test_graphs.py b/pygsp/tests/test_graphs.py index e2f75a38..59bc2dfa 100644 --- a/pygsp/tests/test_graphs.py +++ b/pygsp/tests/test_graphs.py @@ -556,21 +556,22 @@ def test_sphere_gausslegendre(self, nrings=11): self.assertRaises(NotImplementedError, graphs.SphereGaussLegendre, reduced='glesp-equal-area') - def test_sphere_icosahedral(self, subdivisions=3): - n_faces = 20 * 4**subdivisions - n_edges = 30 * 4**subdivisions + def test_sphere_icosahedral(self, subdivisions=8): + n_faces = 20 * subdivisions**2 + n_edges = 30 * subdivisions**2 n_vertices = n_edges - n_faces + 2 graph = graphs.SphereIcosahedral(subdivisions, kind='radius', - radius=1.6/2**subdivisions) + radius=1.6/subdivisions) self.assertEqual(graph.n_vertices, n_vertices) self.assertEqual(graph.n_edges, n_edges) self._test_sphere(graph) - self.assertEqual(np.sum(graph.signals['lat'] == 0), 4*2**subdivisions) + self.assertEqual(np.sum(graph.signals['lat'] == 0), 4*subdivisions) graph = graphs.SphereIcosahedral(subdivisions, dual=True, k=3) self.assertEqual(graph.n_vertices, n_faces) self.assertEqual(graph.n_edges, n_edges) self._test_sphere(graph) - self.assertEqual(np.sum(graph.signals['lat'] == 0), 4*2**subdivisions) + self.assertEqual(np.sum(graph.signals['lat'] == 0), 4*subdivisions) + self.assertRaises(NotImplementedError, graphs.SphereIcosahedral, 3) def test_sphere_cubed(self, subdivisions=9): for spacing in ['equiangular', 'equidistant']: @@ -585,17 +586,20 @@ def test_sphere_cubed(self, subdivisions=9): [-1, 0, 0], [0, -1, 0], [0, 0, -1]] np.testing.assert_allclose(graph.coords, coords) - def test_sphere_healpix(self, subdivisions=2): - nside = 2**subdivisions + def test_sphere_healpix(self, subdivisions=4): graph = graphs.SphereHealpix(subdivisions) - self.assertEqual(graph.n_vertices, 12*nside**2) + self.assertEqual(graph.n_vertices, 12*subdivisions**2) self._test_sphere(graph) graphs.SphereHealpix(subdivisions, k=20) graphs.SphereHealpix(subdivisions, k=21, kernel_width=0.1) self.assertRaises(ValueError, graphs.SphereHealpix, subdivisions, k=21) graphs.SphereHealpix(subdivisions, indexes=range(10), k=8) graphs.SphereHealpix(subdivisions, nest=True) - self.assertEqual(np.sum(graph.signals['lat'] == 0), 4*nside) + self.assertEqual(np.sum(graph.signals['lat'] == 0), 4*subdivisions) + # Subdivisions must be powers of 2 for NESTED ordering. + graphs.SphereHealpix(3, nest=False, kernel_width=1) + self.assertRaises(ValueError, graphs.SphereHealpix, 3, nest=True, + kernel_width=1) def test_twomoons(self): graphs.TwoMoons(moontype='standard')