diff --git a/doc/background/spherical_samplings.rst b/doc/background/spherical_samplings.rst index 0b445b49..ee5529c1 100644 --- a/doc/background/spherical_samplings.rst +++ b/doc/background/spherical_samplings.rst @@ -173,6 +173,19 @@ SHT. .. _geodesic grid: https://en.wikipedia.org/wiki/Geodesic_grid .. _icosahedral symmetry: https://en.wikipedia.org/wiki/Icosahedral_symmetry +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 result is a convex polyhedron made of :math:`n=6⋅m^2` faces. +The graph vertices represent the quadrilateral faces. + +The sampling is used in weather and climate modeling because edges and faces +are of approximately equal length and area (better when +``spacing='equiangular'``). It is hierarchical (for subdivisions arranged in +powers of 2) but doesn't have a sampling theorem nor a fast SHT. + Subdivision of the rhombic dodecahedron (:class:`~pygsp.graphs.SphereHealpix`) ------------------------------------------------------------------------------ diff --git a/pygsp/graphs/__init__.py b/pygsp/graphs/__init__.py index 1697595e..730e4890 100644 --- a/pygsp/graphs/__init__.py +++ b/pygsp/graphs/__init__.py @@ -179,6 +179,7 @@ SphereEquiangular SphereGaussLegendre SphereIcosahedral + SphereCubed SphereHealpix SphereRandom @@ -219,6 +220,7 @@ 'SphereGaussLegendre', 'SphereHealpix', 'SphereIcosahedral', + 'SphereCubed', 'TwoMoons' ] diff --git a/pygsp/graphs/nngraphs/spherecubed.py b/pygsp/graphs/nngraphs/spherecubed.py new file mode 100644 index 00000000..ed5e0cca --- /dev/null +++ b/pygsp/graphs/nngraphs/spherecubed.py @@ -0,0 +1,131 @@ +# -*- coding: utf-8 -*- + +import numpy as np + +from pygsp.graphs import NNGraph # prevent circular import in Python < 3.5 +from pygsp import utils + + +class SphereCubed(NNGraph): + r"""Sphere sampled as a subdivided cube. + + Background information is found at :doc:`/background/spherical_samplings`. + + Parameters + ---------- + subdivisions : int + Number of edges the cube's edges are divided into, resulting in + ``6*subdivisions**2`` vertices. + spacing : {'equiangular', 'equidistant'} + Whether the vertices (on a cube's face) are spaced by equal angles (as + in [4]_) or equal distances (as in [1]_). + kwargs : dict + Additional keyword parameters are passed to :class:`NNGraph`. + + Attributes + ---------- + signals : dict + Vertex position as latitude ``'lat'`` in [-π/2,π/2] and longitude + ``'lon'`` in [0,2π[. + + See Also + -------- + SphereEquiangular, SphereGaussLegendre : based on quadrature theorems + SphereHealpix, SphereIcosahedral : based on subdivided polyhedra + SphereRandom : random uniform sampling + + Notes + ----- + Edge weights are computed by :class:`NNGraph`. Gaussian kernel widths have + however not been optimized for convolutions on the resulting graph to be + maximally equivariant to rotation [7]_. + + References + ---------- + .. [1] R. Sadourny, Conservative finite-difference approximations of the + primitive equations on quasi-uniform spherical grids, 1972. + .. [2] EM O'Neill, RE Laubscher, Extended studies of a quadrilateralized + spherical cube Earth data base, 1976. + .. [3] R. A. White, S. W. Stemwedel, The quadrilateralized spherical cube + and quad-tree for all sky data, 1992. + .. [4] C. Ronchi, R. Iacono, P. Paolucci, The “Cubed-Sphere:” a new method + for the solution of partial differential equations in spherical + geometry, 1995. + .. [5] W. M. Putmana, S.-J. Lin, Finite-volume transport on various + cubed-sphere grids, 2007. + .. [6] https://en.wikipedia.org/wiki/Quadrilateralized_spherical_cube + .. [7] M. Defferrard et al., DeepSphere: a graph-based spherical CNN, 2019. + + Examples + -------- + >>> import matplotlib.pyplot as plt + >>> G = graphs.SphereCubed() + >>> fig = plt.figure() + >>> ax1 = fig.add_subplot(131) + >>> ax2 = fig.add_subplot(132, projection='3d') + >>> ax3 = fig.add_subplot(133) + >>> _ = ax1.spy(G.W, markersize=1.5) + >>> _ = G.plot(ax=ax2) + >>> G.set_coordinates('sphere', dim=2) + >>> _ = G.plot(ax=ax3, indices=True) + + Equiangular and equidistant spacings: + + >>> import matplotlib.pyplot as plt + >>> fig, ax = plt.subplots() + >>> graph = graphs.SphereCubed(4, 'equiangular') + >>> graph.set_coordinates('sphere', dim=2) + >>> _ = graph.plot('C0', edges=False, ax=ax) + >>> graph = graphs.SphereCubed(4, 'equidistant') + >>> graph.set_coordinates('sphere', dim=2) + >>> _ = graph.plot('C1', edges=False, ax=ax) + >>> _ = ax.set_title('Equiangular and equidistant spacings') + + """ + + def __init__(self, subdivisions=3, spacing='equiangular', **kwargs): + + self.subdivisions = subdivisions + self.spacing = spacing + + def linspace(interval): + x = np.linspace(-interval, interval, subdivisions, endpoint=False) + x += interval / subdivisions + return x + + a = np.sqrt(3) / 3 + if spacing == 'equidistant': + x = linspace(a) + y = linspace(a) + elif spacing == 'equiangular': + x = a * np.tan(linspace(np.pi/4)) + y = a * np.tan(linspace(np.pi/4)) + + x, y = np.meshgrid(x, y) + x, y = x.flatten(), y.flatten() + z = a * np.ones_like(x) + + coords = np.concatenate([ + [-y, x, z], # North face. + [z, x, y], # Equatorial face centered at longitude 0. + [-x, z, y], # Equatorial face centered at longitude 90°. + [-z, -x, y], # Equatorial face centered at longitude 180°. + [x, -z, y], # Equatorial face centered at longitude 270°. + [y, x, -z], # South face. + ], axis=1).T + + coords /= np.linalg.norm(coords, axis=1)[:, np.newaxis] + + super().__init__(coords, **kwargs) + + lat, lon = utils.xyz2latlon(*coords.T) + self.signals['lat'] = lat + self.signals['lon'] = lon + + def _get_extra_repr(self): + attrs = { + 'subdivisions': self.subdivisions, + 'spacing': self.spacing, + } + attrs.update(super(SphereCubed, self)._get_extra_repr()) + return attrs diff --git a/pygsp/graphs/nngraphs/spheregausslegendre.py b/pygsp/graphs/nngraphs/spheregausslegendre.py index 1c9ed505..20e8e405 100644 --- a/pygsp/graphs/nngraphs/spheregausslegendre.py +++ b/pygsp/graphs/nngraphs/spheregausslegendre.py @@ -36,7 +36,8 @@ class SphereGaussLegendre(NNGraph): See Also -------- SphereEquiangular : based on quadrature theorems - SphereIcosahedral, SphereHealpix : based on subdivided polyhedra + SphereIcosahedral, SphereCubed, SphereHealpix : + based on subdivided polyhedra SphereRandom : random uniform sampling Notes diff --git a/pygsp/graphs/nngraphs/spherehealpix.py b/pygsp/graphs/nngraphs/spherehealpix.py index 9e16d377..039d9d60 100644 --- a/pygsp/graphs/nngraphs/spherehealpix.py +++ b/pygsp/graphs/nngraphs/spherehealpix.py @@ -46,7 +46,7 @@ class SphereHealpix(NNGraph): See Also -------- SphereEquiangular, SphereGaussLegendre : based on quadrature theorems - SphereIcosahedral : based on subdivided polyhedra + SphereIcosahedral, SphereCubed : based on subdivided polyhedra SphereRandom : random uniform sampling Notes diff --git a/pygsp/graphs/nngraphs/sphereicosahedral.py b/pygsp/graphs/nngraphs/sphereicosahedral.py index 2b196bc0..64b3715a 100644 --- a/pygsp/graphs/nngraphs/sphereicosahedral.py +++ b/pygsp/graphs/nngraphs/sphereicosahedral.py @@ -42,7 +42,7 @@ class SphereIcosahedral(NNGraph): See Also -------- SphereEquiangular, SphereGaussLegendre : based on quadrature theorems - SphereHealpix : based on subdivided polyhedra + SphereCubed, SphereHealpix : based on subdivided polyhedra SphereRandom : random uniform sampling Notes @@ -121,6 +121,7 @@ def normalize(vertices): # Projecting pushes points away from the 12 base vertices, which # may make the point density more uniform. # See "A Comparison of Popular Point Configurations on S^2". + # As the equiangular vs equidistant spacing on the subdivided cube. normalize(mesh.vertices) if not dual: diff --git a/pygsp/graphs/nngraphs/sphererandom.py b/pygsp/graphs/nngraphs/sphererandom.py index 731cb311..0b72d274 100644 --- a/pygsp/graphs/nngraphs/sphererandom.py +++ b/pygsp/graphs/nngraphs/sphererandom.py @@ -29,7 +29,8 @@ class SphereRandom(NNGraph): See Also -------- SphereEquiangular, SphereGaussLegendre : based on quadrature theorems - SphereIcosahedral, SphereHealpix : based on subdivided polyhedra + SphereIcosahedral, SphereCubed, SphereHealpix : + based on subdivided polyhedra CubeRandom : randomly sampled cube References diff --git a/pygsp/graphs/sphereequiangular.py b/pygsp/graphs/sphereequiangular.py index 7492fd05..507bf455 100644 --- a/pygsp/graphs/sphereequiangular.py +++ b/pygsp/graphs/sphereequiangular.py @@ -38,7 +38,8 @@ class SphereEquiangular(Graph): See Also -------- SphereGaussLegendre : based on quadrature theorems - SphereIcosahedral, SphereHealpix : based on subdivided polyhedra + SphereIcosahedral, SphereCubed, SphereHealpix : + based on subdivided polyhedra SphereRandom : random uniform sampling Notes diff --git a/pygsp/tests/test_graphs.py b/pygsp/tests/test_graphs.py index 098c9681..e2f75a38 100644 --- a/pygsp/tests/test_graphs.py +++ b/pygsp/tests/test_graphs.py @@ -572,6 +572,19 @@ def test_sphere_icosahedral(self, subdivisions=3): self._test_sphere(graph) self.assertEqual(np.sum(graph.signals['lat'] == 0), 4*2**subdivisions) + def test_sphere_cubed(self, subdivisions=9): + for spacing in ['equiangular', 'equidistant']: + graph = graphs.SphereCubed(subdivisions, spacing) + self.assertEqual(graph.n_vertices, 6*subdivisions**2) + self._test_sphere(graph) + n_pixels_at_equator = sum(abs(graph.signals['lat']) < 1e-10) + self.assertEqual(n_pixels_at_equator, 4*subdivisions) + # Ordering of the faces. + graph = graphs.SphereCubed(1, k=2) + coords = [[0, 0, 1], [1, 0, 0], [0, 1, 0], + [-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 graph = graphs.SphereHealpix(subdivisions)