.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/features/plot_point_moments.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_examples_features_plot_point_moments.py: .. _point_moments_example: How to compute local geometric moments ======================================= We use the :meth:`~skshapes.PolyData.point_moments` method to compute local averages and covariance matrices. This is useful for estimating point normals and curvatures. .. GENERATED FROM PYTHON SOURCE LINES 11-13 First, we create a simple point cloud in 2D: .. GENERATED FROM PYTHON SOURCE LINES 13-23 .. code-block:: Python import torch from matplotlib import pyplot as plt import skshapes as sks n_points, dim = 3, 2 shape = sks.PolyData(torch.rand(n_points, dim)) print(shape) .. rst-class:: sphx-glr-script-out .. code-block:: none skshapes.PolyData (0x7f6d1790f6d0 on cpu, float32), a 2D point cloud with: - 3 points - center (0.345, 0.390), radius 0.384 - bounds x=(0.201, 0.608), y=(0.111, 0.668), .. GENERATED FROM PYTHON SOURCE LINES 24-36 The method :meth:`~skshapes.PolyData.point_neighborhoods` associates to each point :math:`x_i` in the shape a weight distribution :math:`\nu_i(x)` over the ambient space. Typically, this corresponds to a normalized Gaussian window at scale :math:`\sigma`. If :math:`\mu(x)` denotes the uniform distribution on the support of the shape, then: .. math:: \nu_i(x) \propto \exp\left(-\frac{||x-x_i||^2}{2\sigma^2}\right) \mu(x) The method :meth:`~skshapes.PolyData.point_moments` takes as input the same arguments, and returns a :class:`skshapes.features.Moments` object that computes the moments of order 0, 1 and 2 of these distributions. .. GENERATED FROM PYTHON SOURCE LINES 36-39 .. code-block:: Python moments = shape.point_moments(scale=0.5) .. GENERATED FROM PYTHON SOURCE LINES 40-41 The total mass :math:`m_i` of :math:`\nu_i` is an estimate of the local point density: .. GENERATED FROM PYTHON SOURCE LINES 41-43 .. code-block:: Python print(moments.masses) .. rst-class:: sphx-glr-script-out .. code-block:: none tensor([1.8844, 2.2160, 2.6714]) .. GENERATED FROM PYTHON SOURCE LINES 44-46 The local mean :math:`\overline{x}_i = \mathbb{E}_{x\sim \nu_i}[x]` is a point in the vicinity of :math:`x_i`: .. GENERATED FROM PYTHON SOURCE LINES 46-48 .. code-block:: Python print(moments.means) .. rst-class:: sphx-glr-script-out .. code-block:: none tensor([[0.4216, 0.2971], [0.2916, 0.4616], [0.3206, 0.4127]]) .. GENERATED FROM PYTHON SOURCE LINES 49-52 The local covariance matrix :math:`\Sigma_i = \mathbb{E}_{x\sim \nu_i}[(x-\overline{x}_i)(x-\overline{x}_i)^{\intercal}]` is symmetric and positive semi-definite: .. GENERATED FROM PYTHON SOURCE LINES 52-54 .. code-block:: Python print(moments.covariances) .. rst-class:: sphx-glr-script-out .. code-block:: none tensor([[[ 0.0393, -0.0384], [-0.0384, 0.0478]], [[ 0.0245, -0.0257], [-0.0257, 0.0450]], [[ 0.0311, -0.0313], [-0.0313, 0.0480]]]) .. GENERATED FROM PYTHON SOURCE LINES 55-56 Its eigenvalues are non-negative and typically of order :math:`\sigma^2`: .. GENERATED FROM PYTHON SOURCE LINES 56-60 .. code-block:: Python L = moments.covariance_eigenvalues print(L) .. rst-class:: sphx-glr-script-out .. code-block:: none tensor([[0.0049, 0.0822], [0.0071, 0.0624], [0.0071, 0.0720]]) .. GENERATED FROM PYTHON SOURCE LINES 61-63 Its eigenvectors are orthogonal and represent the principal directions of the local distribution of points: .. GENERATED FROM PYTHON SOURCE LINES 63-67 .. code-block:: Python Q = moments.covariance_eigenvectors print(Q) .. rst-class:: sphx-glr-script-out .. code-block:: none tensor([[[-0.7452, -0.6668], [-0.6668, 0.7452]], [[-0.8279, -0.5608], [-0.5608, 0.8279]], [[-0.7938, -0.6081], [-0.6081, 0.7938]]]) .. GENERATED FROM PYTHON SOURCE LINES 68-69 The eigenvectors are stored column-wise, and sorted by increasing eigenvalue: .. GENERATED FROM PYTHON SOURCE LINES 69-75 .. code-block:: Python LQt = L.view(n_points, dim, 1) * Q.transpose(1, 2) QLQt = Q @ LQt print(f"Reconstruction error: {(QLQt - moments.covariances).abs().max():.2e}") .. rst-class:: sphx-glr-script-out .. code-block:: none Reconstruction error: 3.73e-09 .. GENERATED FROM PYTHON SOURCE LINES 76-83 .. note:: These attributes are computed when required, and cached in memory. On a 2D curve -------------- Going further, let's consider a sampled curve in 2D: .. GENERATED FROM PYTHON SOURCE LINES 83-87 .. code-block:: Python shape = sks.doc.wiggly_curve(n_points=32, dim=2) print(shape) .. rst-class:: sphx-glr-script-out .. code-block:: none skshapes.PolyData (0x7f6d26ea60d0 on cpu, float32), a 2D point cloud with: - 32 points - center (-0.078, 0.224), radius 1.283 - bounds x=(-1.076, 1.196), y=(-0.923, 1.098), .. GENERATED FROM PYTHON SOURCE LINES 88-104 Intuitively, computing local moments :math:`m_i` of order 0, :math:`\overline{x}_i` of order 1 and :math:`\Sigma_i` of order 2 is equivalent to fitting a Gaussian distribution to each local neighborhood :math:`\nu_i`. In dimension :math:`D=2` or :math:`D=3`, if :math:`\lambda_{\mathbb{R}^D}` denotes the Lebesgue measure on the ambient space (i.e. the area or the volume), then: .. math:: \text{d}\nu_i(x) ~&\simeq~ m_i\,\text{d}\mathcal{N}(\overline{x}_i, \Sigma_i)(x) \\ &= ~\frac{m_i}{(2\pi\sigma^2)^{D/2} } \, \exp\left[-\tfrac{1}{2}(x-\overline{x}_i)^{\intercal}\Sigma_i^{-1}(x-\overline{x}_i)\right]\, \text{d}\lambda_{\mathbb{R}^D}(x)~ We visualize these descriptors by drawing ellipses centered at :math:`\overline{x}_i`, colored by the local densities :math:`m_i` and oriented in :math:`\mathbb{R}^2` or :math:`\mathbb{R}^3` with axes that are aligned with the eigenvectors of :math:`\Sigma_i` and whose lengths are proportional to the square root of its eigenvalues: .. GENERATED FROM PYTHON SOURCE LINES 104-117 .. code-block:: Python plt.figure(figsize=(12, 10)) for i, sigma in enumerate([0.1, 0.2, 0.5, 1.0]): moments = shape.point_moments(scale=sigma) # Plot the ellipses ax = plt.subplot(2, 2, i + 1) ax.set_title(f"Moments at scale = {sigma:.1f}") sks.doc.display_covariances(shape.points, moments) plt.tight_layout() .. image-sg:: /auto_examples/features/images/sphx_glr_plot_point_moments_001.png :alt: Moments at scale = 0.1, Moments at scale = 0.2, Moments at scale = 0.5, Moments at scale = 1.0 :srcset: /auto_examples/features/images/sphx_glr_plot_point_moments_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 118-126 As evidenced here, :meth:`~skshapes.PolyData.point_moments` describes the local shape context at scale :math:`\sigma`. We use it to compute point normals and curvatures that are robust to noise and sampling artifacts. On a 3D surface --------------- We can perform the same analysis on a 3D surface: .. GENERATED FROM PYTHON SOURCE LINES 126-136 .. code-block:: Python import pyvista as pv shape = ( sks.PolyData(pv.examples.download_bunny()) .resample(n_points=5000) .normalize() ) print(shape) .. rst-class:: sphx-glr-script-out .. code-block:: none skshapes.PolyData (0x7f6d16e79050 on cpu, float32), a 3D triangle mesh with: - 5,000 points, 14,952 edges, 9,949 triangles - center (-0.000, 0.000, 0.000), radius 1.000 - bounds x=(-0.581, 0.767), y=(-0.534, 0.807), z=(-0.597, 0.453) .. GENERATED FROM PYTHON SOURCE LINES 137-138 As expected, the point moments of order 1 and 2 now refer to 3D vectors and matrices: .. GENERATED FROM PYTHON SOURCE LINES 138-144 .. code-block:: Python moments = shape.point_moments(scale=0.05) print("moments.masses: ", moments.masses.shape) print("moments.means: ", moments.means.shape) print("moments.covariances:", moments.covariances.shape) .. rst-class:: sphx-glr-script-out .. code-block:: none moments.masses: torch.Size([5000]) moments.means: torch.Size([5000, 3]) moments.covariances: torch.Size([5000, 3, 3]) .. GENERATED FROM PYTHON SOURCE LINES 145-146 We visualize the ellipsoids on a subset of our full point cloud: .. GENERATED FROM PYTHON SOURCE LINES 146-157 .. code-block:: Python landmark_indices = list(range(0, shape.n_points, 50)) pl = pv.Plotter() sks.doc.display(plotter=pl, shape=shape, opacity=0.5) sks.doc.display_covariances( shape.points, moments, landmark_indices=landmark_indices, plotter=pl ) pl.show() .. image-sg:: /auto_examples/features/images/sphx_glr_plot_point_moments_002.png :alt: plot point moments :srcset: /auto_examples/features/images/sphx_glr_plot_point_moments_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 158-164 Using the covariance eigenvalues :math:`\lambda_1 \leqslant \lambda_2 \leqslant \lambda_3`, we can compute simple shape descriptors such as the plateness which is equal to 0 if the local ellipsoid is a sphere and 1 if it is a 2D disk: .. math:: \text{plateness} = 1 - \frac{\lambda_1}{\sqrt{\lambda_2\lambda_3}}~. .. GENERATED FROM PYTHON SOURCE LINES 164-188 .. code-block:: Python scales = [0.05, 0.1, 0.15, 0.2] pl = pv.Plotter(shape=(2, 2)) for i, scale in enumerate(scales): moments = shape.point_moments(scale=scale) eigs = moments.covariance_eigenvalues shape.point_data[f"plateness_{i}"] = ( 1 - eigs[:, 0] / (eigs[:, 1] * eigs[:, 2]).sqrt() ) pl.subplot(i // 2, i % 2) sks.doc.display( plotter=pl, shape=shape, scalars=f"plateness_{i}", scalar_bar=True, smooth=1, clim=(0, 1), title=f"Plateness at scale {scale:.2f}", ) pl.show() .. image-sg:: /auto_examples/features/images/sphx_glr_plot_point_moments_003.png :alt: plot point moments :srcset: /auto_examples/features/images/sphx_glr_plot_point_moments_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 189-194 Likewise, we can compute the tubeness, which is equal to 0 if the local ellipsoid is a sphere or a disk, and 1 if it is a 1D segment: .. math:: \text{tubeness} = 1 - \frac{\lambda_2}{\lambda_3}~. .. GENERATED FROM PYTHON SOURCE LINES 194-212 .. code-block:: Python pl = pv.Plotter(shape=(2, 2)) for i, scale in enumerate(scales): moments = shape.point_moments(scale=scale) eigs = moments.covariance_eigenvalues shape.point_data[f"tubeness_{i}"] = 1 - (eigs[:, 1] / eigs[:, 2]) pl.subplot(i // 2, i % 2) sks.doc.display( plotter=pl, shape=shape, scalars=f"tubeness_{i}", scalar_bar=True, smooth=1, clim=(0, 1), title=f"Tubeness at scale {scale:.2f}", ) pl.show() .. image-sg:: /auto_examples/features/images/sphx_glr_plot_point_moments_004.png :alt: plot point moments :srcset: /auto_examples/features/images/sphx_glr_plot_point_moments_004.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 23.067 seconds) .. _sphx_glr_download_auto_examples_features_plot_point_moments.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_point_moments.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_point_moments.py ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_