Source code for pycsphere.mesh

# #############################################################################
# mesh.py
# =======
# Author : Matthieu Simeoni [matthieu.simeoni@gmail.com]
# #############################################################################

r"""
Routines for deterministic and random spherical point sets.
"""

import numpy as np
from abc import ABC, abstractmethod
import scipy.spatial as sp
import scipy.sparse as sparse
import matplotlib.pyplot as plt
import mpl_toolkits.mplot3d.axes3d as plt3
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
import healpy as hp
from matplotlib.cm import get_cmap
import matplotlib.colors as mplcols
from typing import Optional, List, Union, Tuple
from astropy.coordinates import Angle
import astropy.units as u
from scipy.special import lambertw


[docs]class SphericalPointSet(ABC): r""" Base class for spherical point sets. Any subclass/instance of this class must implement the abstract method ``generate_point_set`` which returns the point set cartesian/spherical coordinates. A spherical point set is a collection :math:`\Theta_N=\{\mathbf{r}_1, \ldots, \mathbf{r}_N\}\subset \mathbb{S}^2` of unit vectors (directions). Notes ----- The class ``SphericalPointSet`` possesses a few class/object attributes which allow to describe the properties of the point set. The terminology is the same as [PointSets]_: * A sequence of point sets :math:`\{\Theta_N\}_{n\in\mathbb{N}}\subset \mathcal{P}(\mathbb{S}^2)` is called ``equidistributed`` if the sequence of *normalised atomic measures* :math:`\nu_N(A):=|A\cap \Theta_N|/N` converges w.r.t. the weak*-topology towards the Lebesgue measure as :math:`N` grows to infinity. * The *separation* of a point set is defined as :math:`\delta(\Theta_N):=\min_{\mathbf{r}\neq \mathbf{s}\in\Theta_N} \|\mathbf{r}- \mathbf{s}\|_2`. A set is said ``well_separated`` if :math:`\delta(\Theta_N)\geq c N^{-1/2}` for some ``separation_cst`` constant :math:`c>0` and :math:`N\geq 2`. * The ``nodal_width`` of a point set is defined as :math:`\eta(\Theta_N):=\max_{\mathbf{s}\in\mathbb{S}^2}\min_{\mathbf{r}\in\Theta_N} \|\mathbf{r}- \mathbf{s}\|_2`. A set is said to have ``good_covering`` if :math:`\eta(\Theta_N)\leq C N^{-1/2}` for some ``nodal_width_cst`` :math:`C>0` and :math:`N\geq 2`. * The ``mesh_ratio`` of a point set is the ratio between its nodal width and separation. If the mesh ratio is asymptotically bounded, then the point set is said ``quasi_uniform``. A lower bound on the mesh ratio is ``minimal_mesh_ratio = (np.sqrt(5) - 1) / 2``. """ minimal_mesh_ratio = (np.sqrt(5) - 1) / 2 separation_cst = None nodal_width_cst = None mesh_ratio = None well_separated = None good_covering = None quasi_uniform = None equidistributed = None equal_area = None isolatitude = None hierarchical = None
[docs] @classmethod def arc2angle(cls, arc_length: Union[np.ndarray, float], deg: bool = False) -> Union[np.ndarray, float]: r""" Convert an arc length to an angular section. Parameters ---------- arc_length: Union[np.ndarray, float] Input arc length. deg: bool Whether the output angle is in degrees (``True``) or radians (``False``, default). Returns ------- Union[np.ndarray, float] Angle corresponding to the arc. """ if deg is False: return 2 * np.arcsin(arc_length / 2) else: return 2 * np.arcsin(arc_length / 2) * 180 / np.pi
[docs] @classmethod def angle2arc(cls, angle: Union[np.ndarray, float], deg: bool = False) -> Union[np.ndarray, float]: r""" Compute the arc length of a given angular section. Parameters ---------- angle: Union[np.ndarray, float] Angle defining the arc. deg: bool Whether the input angle is in degrees (``True``) or radians (``False``, default). Returns ------- Union[np.ndarray, float] Arc length of the angular section. """ if deg is True: return 2 * np.sin((np.pi / 180) * angle / 2) else: return 2 * np.sin(angle / 2)
[docs] @classmethod def angle2res(cls, angular_resolution: float, deg: bool = True) -> int: r""" Compute the minimal point set size required to achieve a desired (angular) nodal width. Parameters ---------- angular_resolution: float Desired angular nodal width. deg: bool Whether the angular nodal width is in degrees or not. Returns ------- int Minimal point set resolution. """ nodal_width = cls.angle2arc(angular_resolution, deg=deg) return np.ceil(cls.nodal_width_cst / nodal_width) ** 2
[docs] def __init__(self, resolution: int, separation: Optional[float] = None, nodal_width: Optional[float] = None, lonlat: bool = False): r""" Parameters ---------- resolution: int Resolution (size) of the point set. separation: Optional[float] Separation of the point set (``None`` if unknown). nodal_width: Optional[float] Nodal width of the point set (``None`` if unknown). lonlat: bool If ``True``, longitude/latitude spherical coordinates (in degrees) are used. Otherwise, longitude/co-latitude (in radians) are used. """ self.resolution = resolution self.separation = separation self.nodal_width = nodal_width self._lonlat = lonlat self.vec, self.dir = self.generate_point_set()
[docs] @abstractmethod def generate_point_set(self) -> Tuple[np.ndarray, np.ndarray]: r""" Generate the cartesian/spherical coordinates of the point set. Returns ------- Tuple[np.ndarray, np.ndarray] Cartesian/spherical coordinates of the point set. The latter are also stored in the attributes ``self.vec`` and ``self.dir`` respectively. The attribute ``self.vec`` has shape ``(3, N)`` and contains the cartesian coordinates of the points in the set. The attribute ``self.dir`` has shape ``(2, N)`` and contains the spherical coordinates of the points in the set: * If ``self.latlon==True``, longitudes and latitudes in degrees, * If ``self.latlon==False``, longitudes and co-latitudes in radians. """ pass
[docs] def plot_tessellation(self, facecmap: Optional[str] = 'Blues_r', ncols: int = 20, edgecolor: str = '#32519D', cell_centers: bool = True, markercolor: str = 'k', markersize: float = 2, linewidths=1, alpha: float = 0.9, seed: int = 0): r""" Plots the tessellation defined by the point set. Parameters ---------- facecmap: Optional[str] Color map to be used for random coloring the tessellation faces. ncols: int Number of colors in the color map. edgecolor: str Color of the cell edges. cell_centers: bool Whether or not the cell centers are displayed. markercolor: str Color of the cell centers. markersize: float Size of the cell centers. linewidths: float Width of the cell edges. alpha: float Transparency of the cell faces. seed: int Seed for random coloring reproducibility. Warnings -------- This method can be **very slow** for large point sets (>10'000 points)! """ self.plot_voronoi_cells(facecmap=facecmap, ncols=ncols, edgecolor=edgecolor, cell_centers=cell_centers, markercolor=markercolor, markersize=markersize, linewidths=linewidths, alpha=alpha, seed=seed)
[docs] def plot_delaunay_cells(self, facecmap: Optional[str] = 'Blues_r', ncols: int = 20, edgecolor: str = '#32519D', cell_centers: bool = True, markercolor: str = 'k', markersize: float = 2, linewidths=1, alpha: float = 0.9, seed: int = 0): r""" Plots the tessellation defined by the set's Delaunay triangulation. Parameters ---------- facecmap: Optional[str] Color map to be used for random coloring the tessellation faces. ncols: int Number of colors in the color map. edgecolor: str Color of the cell edges. cell_centers: bool Whether or not the cell centers are displayed. markercolor: str Color of the cell centers. markersize: float Size of the cell centers. linewidths: float Width of the cell edges. alpha: float Transparency of the cell faces. seed: int Seed for random coloring reproducibility. Warnings -------- This method can be **very slow** for large point sets (>10'000 points)! """ vec = self.vec.transpose() delaunay = sp.ConvexHull(vec) vertices = [[vec[delaunay.simplices[i]]] for i in range(delaunay.simplices.shape[0])] self._plot_spherical_polygons(vertices, facecmap=facecmap, ncols=ncols, edgecolor=edgecolor, cell_centers=cell_centers, markercolor=markercolor, markersize=markersize, linewidths=linewidths, alpha=alpha, seed=seed)
[docs] def plot_voronoi_cells(self, facecmap: Optional[str] = 'Blues_r', ncols: int = 20, edgecolor: str = '#32519D', cell_centers: bool = True, markercolor: str = 'k', markersize: float = 2, linewidths=1, alpha: float = 0.9, seed: int = 0): r""" Plots the tessellation defined by the set's Voronoi decomposition. Parameters ---------- facecmap: Optional[str] Color map to be used for random coloring the tessellation faces. ncols: int Number of colors in the color map. edgecolor: str Color of the cell edges. cell_centers: bool Whether or not the cell centers are displayed. markercolor: str Color of the cell centers. markersize: float Size of the cell centers. linewidths: float Width of the cell edges. alpha: float Transparency of the cell faces. seed: int Seed for random coloring reproducibility. Warnings -------- This method can be **very slow** for large point sets (>10'000 points)! """ voronoi = sp.SphericalVoronoi(self.vec.transpose(), radius=1) voronoi.sort_vertices_of_regions() vertices = [[voronoi.vertices[region]] for region in voronoi.regions] self._plot_spherical_polygons(vertices, facecmap=facecmap, ncols=ncols, edgecolor=edgecolor, cell_centers=cell_centers, markercolor=markercolor, markersize=markersize, linewidths=linewidths, alpha=alpha, seed=seed)
def _plot_spherical_polygons(self, vertices: List[List[np.ndarray]], facecmap: Optional[str] = 'Blues_r', ncols: int = 20, edgecolor: str = '#32519D', cell_centers: bool = True, markercolor: str = 'k', markersize: float = 3, linewidths=1, alpha: float = 0.9, seed: int = 0): rng = np.random.default_rng(seed=seed) cmap_obj = get_cmap(facecmap, ncols) normalise_data = mplcols.Normalize(vmin=0, vmax=ncols - 1) fig = plt.figure() ax3 = plt3.Axes3D(fig) ax3.scatter3D(1, 1, 1, s=0) ax3.scatter3D(-1, -1, -1, s=0) ax3.view_init(elev=80, azim=0) if cell_centers is True: ax3.scatter3D(self.vec[0], self.vec[1], self.vec[2], '.', color=markercolor, s=markersize) for verts in vertices: color_index = rng.integers(0, ncols - 1, size=1) polygon = Poly3DCollection(verts, linewidths=linewidths, alpha=alpha) polygon.set_facecolor(cmap_obj(normalise_data(color_index))) polygon.set_edgecolor(edgecolor) ax3.add_collection3d(polygon) plt.axis('off') plt.show() @property def angular_nodal_width(self) -> Optional[float]: r""" Angular nodal width. """ if self.nodal_width is None: return None else: return self.arc2angle(self.nodal_width, deg=self.lonlat) @property def lonlat(self) -> bool: return self._lonlat @lonlat.setter def lonlat(self, value): self._lonlat = value self.dir = hp.vec2dir(self.vec, lonlat=value)
[docs] def compute_empirical_nodal_width(self, mode='mean') -> float: r""" Compute the set's nodal width. """ cvx_hull = sp.ConvexHull(self.vec.transpose()) cols = np.roll(cvx_hull.simplices, shift=1, axis=-1).reshape(-1) rows = cvx_hull.simplices.reshape(-1) # Form sparse coo_matrix from extracted pairs affinity_matrix = sparse.coo_matrix((cols * 0 + 1, (rows, cols)), shape=(cvx_hull.points.shape[0], cvx_hull.points.shape[0])) # Symmetrize the matrix to obtain an undirected graph. extended_row = np.concatenate([affinity_matrix.row, affinity_matrix.col]) extended_col = np.concatenate([affinity_matrix.col, affinity_matrix.row]) affinity_matrix.row = extended_row affinity_matrix.col = extended_col affinity_matrix.data = np.concatenate([affinity_matrix.data, affinity_matrix.data]) affinity_matrix = affinity_matrix.tocsr().tocoo() # Delete potential duplicate pairs distance = np.linalg.norm(cvx_hull.points[affinity_matrix.row, :] - cvx_hull.points[affinity_matrix.col, :], axis=-1) if mode is 'mean': nodal_distance = np.mean(distance) # average distance to neighbors elif mode is 'max': nodal_distance = np.max(distance) elif mode is 'median': nodal_distance = np.median(distance) else: raise TypeError("Parameter mode must be one of ['mean', 'max', 'median']") return nodal_distance
[docs]class FibonacciPointSet(SphericalPointSet): r""" Fibonacci point set. Examples -------- .. plot:: from pycsphere.mesh import FibonacciPointSet N = FibonacciPointSet.angle2N(angular_resolution=10) fib = FibonacciPointSet(N, lonlat=True) fib.plot_delaunay_cells() fib.plot_tessellation() Notes ----- Points in the Fibonacci point set are arranged uniformly along a spiral pattern on the sphere linking the two poles (see Section 2 of [PointSets]_ for a definition of the point set). The Fibonacci point set is defined for odd resolutions :math:`2N+1` only. It is well-separated, has good covering and is hence quasi-uniform. It is also equidistributed. Its Voronoi cells are however irregular and do not have the same shape/area. The points are not arranged on isolatitude tracks. Finally it is non hierarchical (Fibonacci point sets at two different resolutions are not contained in one another). .. todo:: Add support for local Fibonacci meshes. """ separation_cst = 3.09206862 nodal_width_cst = 2.72812463 mesh_ratio = 0.882298 well_separated = True good_covering = True quasi_uniform = True equidistributed = True equal_area = False isolatitude = False hierarchical = False
[docs] @classmethod def angle2N(cls, angular_resolution: float, deg: bool = True) -> int: r""" Minimal parameter :math:`N` to achieve a prescribed angular nodal width. Parameters ---------- angular_resolution: float Desired angular nodal width. deg: bool Whether or not the nodal width is provided in degrees. Returns ------- int Minimal value of the parameter :math:`N` defining the point set resolution :math:`2N+ 1`. """ resolution = cls.angle2res(angular_resolution=angular_resolution, deg=deg) return np.ceil((resolution - 1) / 2)
[docs] def __init__(self, N: int, lonlat: bool = False): r""" Parameters ---------- N: int Defines the point set resolution :math:`2N+1`. lonlat: bool Convention for the spherical coordinates (see help of :py:func:`~pycsphere.mesh.SphericalPointSet.__init__`). """ separation = self.separation_cst / np.sqrt(2 * N + 1) nodal_width = self.nodal_width_cst / np.sqrt(2 * N + 1) self.golden_ratio = (1 + np.sqrt(5)) / 2 self.golden_angle = 2 * np.pi * (1 - 1 / self.golden_ratio) super(FibonacciPointSet, self).__init__(resolution=2 * N + 1, separation=separation, nodal_width=nodal_width, lonlat=lonlat)
[docs] def generate_point_set(self) -> Tuple[np.ndarray, np.ndarray]: step = np.arange(self.resolution) phi = step * self.golden_angle phi = Angle(phi * u.rad).wrap_at(2 * np.pi * u.rad) theta = np.arccos(1 - (2 * step / self.resolution)) vec = hp.dir2vec(theta=theta, phi=phi, lonlat=False) dir = hp.vec2dir(vec, lonlat=self.lonlat) return vec, dir
[docs]class HEALPixPointSet(SphericalPointSet): r""" HEALPix point set. Examples -------- .. plot:: from pycsphere.mesh import HEALPixPointSet nside = HEALPixPointSet.angle2nside(angular_resolution=10) healpix = HEALPixPointSet(nside, lonlat=True) healpix.plot_delaunay_cells() healpix.plot_voronoi_cells() healpix.plot_tessellation() Notes ----- Developed by NASA for fast data analysis of the cosmic microwave background (CMB), the *Hierarchical Equal Area iso-Latitude Pixelization (HEALPix)* was designed to have three properties essential for computational efficiency in discretising functions on :math:`\mathbb{S}^2`: 1. The sphere is hierarchically tessellated into curvilinear quadrilaterals. 2. The pixelisation is an equal area partition of :math:`\mathbb{S}^2`. 3. The point sets are distributed along iso-latitude lines. A formal definition of the point set is provided Section 2 of [PointSets]_. We use here the :py:mod:`healpy` package. The HEALPix point set is defined for resolutions :math:`N=12k^2` only (parameter :math:`k` is called ``nside``). It is well-separated, has good covering and is hence quasi-uniform. It is also equidistributed. Its tessellation cells are moreover regular and have the same shape/area. Finally it is hierarchical (point sets at two different resolutions are nested in one another). .. todo:: Add spatial filtering routines for local HEALPix meshes. """ separation_cst = 2.8345 nodal_width_cst = 2.8345 mesh_ratio = 1 well_separated = True good_covering = True quasi_uniform = True equidistributed = True equal_area = True isolatitude = True hierarchical = True
[docs] @classmethod def angle2nside(cls, angular_resolution: float, deg: bool = True) -> int: r""" Minimal parameter :math:`k` (``nside``) to achieve a prescribed angular nodal width. Parameters ---------- angular_resolution: float Desired angular nodal width. deg: bool Whether or not the nodal width is provided in degrees. Returns ------- int Minimal value of the parameter :math:`k` defining the point set resolution :math:`N=12k^2`. """ resolution = cls.angle2res(angular_resolution=angular_resolution, deg=deg) return (2 ** np.ceil(np.log2(np.sqrt(resolution / 12)))).astype(int)
[docs] def __init__(self, nside: int, lonlat: bool = False, nest=False): r""" Parameters ---------- nside: int Defines the point set resolution ``12 * nside ** 2``. ``nside`` must be a power of two. lonlat: bool Convention for the spherical coordinates (see help of :py:func:`~pycsphere.mesh.SphericalPointSet.__init__`). nest: bool If ``True``, the HEALPix mesh is stored in NESTED ordering (efficient for pixel querying routines). If ``False``, the HEALPix mesh is stored in RING ordering (efficient for spherical harmonics transforms). """ resolution = int(hp.nside2npix(nside)) separation = self.separation_cst / resolution nodal_width = self.nodal_width_cst / resolution self.nside = nside self.nrings = 4 * self.nside - 1 self.nest = nest super(HEALPixPointSet, self).__init__(resolution=resolution, separation=separation, nodal_width=nodal_width, lonlat=lonlat)
[docs] def generate_point_set(self) -> Tuple[np.ndarray, np.ndarray]: vec = np.stack(hp.pix2vec(self.nside, np.arange(self.resolution)), axis=0) dir = hp.vec2dir(vec, lonlat=self.lonlat) return vec, dir
[docs] def plot_tessellation(self, facecmap: Optional[str] = 'Blues_r', ncols: int = 20, edgecolor: str = '#32519D', cell_centers: bool = True, markercolor: str = 'k', markersize: float = 2, linewidths=1, alpha: float = 0.9, seed: int = 0): vertices = [[hp.boundaries(self.nside, i).transpose()] for i in range(self.resolution)] self._plot_spherical_polygons(vertices, facecmap=facecmap, ncols=ncols, edgecolor=edgecolor, cell_centers=cell_centers, markercolor=markercolor, markersize=markersize, linewidths=linewidths, alpha=alpha, seed=seed)
[docs]class RandomPointSet(SphericalPointSet): r""" Random uniform point set. Examples -------- .. plot:: from pycsphere.mesh import RandomPointSet res = RandomPointSet.angle2res(angular_resolution=10) rnd = RandomPointSet(res, lonlat=True) rnd.plot_delaunay_cells() rnd.plot_tessellation() Notes ----- Properties of spherical random point sets are discussed in Section 2 of [PointSets]_. They are equidistributed but not quasi-uniform, well-separated or with good-coverage. They are also of course not hierarchical and yield highly irregular tessellation cells with different shapes and areas. """ separation_cst = np.sqrt(2 * np.pi) nodal_width_cst = 2 mesh_ratio = None well_separated = False good_covering = False quasi_uniform = False equidistributed = True equal_area = False isolatitude = False hierarchical = False
[docs] @classmethod def angle2res(cls, angular_resolution: float, deg: bool = True) -> int: nodal_width = cls.angle2arc(angular_resolution, deg=deg) return np.ceil(np.exp(-np.real(lambertw(-1 / (cls.nodal_width_cst / nodal_width) ** 2, k=-1)))).astype(int)
[docs] def __init__(self, N: int, seed: int = 0, lonlat: bool = False): r""" Parameters ---------- N: int Resolution of the random point set. seed: int Seed for reproducibility of the random point set. lonlat: bool Convention for the spherical coordinates (see help of :py:func:`~pycsphere.mesh.SphericalPointSet.__init__`). """ resolution = N separation = self.separation_cst / resolution nodal_width = self.nodal_width_cst / np.sqrt(resolution / np.log(resolution)) self.seed = seed super(RandomPointSet, self).__init__(resolution=resolution, separation=separation, nodal_width=nodal_width, lonlat=lonlat)
[docs] def generate_point_set(self) -> Tuple[np.ndarray, np.ndarray]: rng = np.random.default_rng(self.seed) lon = 360 * rng.random(size=self.resolution) - 180 uniform = rng.random(size=self.resolution) lat = np.arcsin(2 * uniform - 1) * 180 / np.pi vec = hp.dir2vec(theta=lon, phi=lat, lonlat=True) dir = hp.vec2dir(vec, lonlat=self.lonlat) return vec, dir