Source code for geomstats.geometry.discrete_surfaces
"""Discrete Surfaces with Elastic metrics.
Lead authors: Emmanuel Hartman, Adele Myers.
+
+References
+----------
+.. [HSKCB2022] Emmanuel Hartman, Yashil Sukurdeep, Eric Klassen,
+ Nicolas Charon, and Martin Bauer.
+ "Elastic shape analysis of surfaces with second-order Sobolev metrics:
+ a comprehensive numerical framework". arXiv:2204.04238 [cs.CV], 25 Sep 2022
+.. [HPBDC2023] Emmanuel Hartman, Emery Pierson, Martin Bauer,
+ Mohamed Daoudi, and Nicolas Charon.
+ “Basis Restricted Elastic Shape Analysis on the Space of Unregistered Surfaces.”
+ arXiv, November 7, 2023. https://doi.org/10.48550/arXiv.2311.04382."""importmath
@@ -780,8 +791,10 @@
Source code for geomstats.geometry.discrete_surfaces
return gs.sqrt(gs.linalg.det(surface_metrics_bp))
+
+[docs]@staticmethod
- def_surface_metric_matrices_from_one_forms(one_forms):
+ defsurface_metric_matrices_from_one_forms(one_forms):"""Compute the surface metric matrices directly from the one_forms. This function is useful for efficiency purposes.
@@ -797,7 +810,8 @@
Source code for geomstats.geometry.discrete_surfaces
Surface metric matrices evaluated at each face of
the triangulated surface. """
- returngs.matmul(one_forms,Matrices.transpose(one_forms))
+ returngs.matmul(one_forms,Matrices.transpose(one_forms))
Source code for geomstats.geometry.discrete_surfaces
:math:`\Delta_q = - Tr(g_q^{-1} \nabla^2)`
where :math:`g_q` is the surface metric matrix of :math:`q`.
+ The area of the triangles is computed using Heron's formula.
+
Parameters ---------- point : array-like, shape=[..., n_vertices, 3]
@@ -880,12 +896,18 @@
Source code for geomstats.geometry.discrete_surfaces
cot_flatten =gs.expand_dims(gs.reshape(cot,point.shape[:-2]+(-1,)),axis=-1)def_laplacian(tangent_vec):
-"""Evaluate the mesh Laplacian operator.
+r"""Evaluate the mesh Laplacian operator. The operator is evaluated at a tangent vector at point to the manifold of DiscreteSurfaces. In other words, the operator is evaluated at a vector field defined on the surface point.
+ .. math::
+
+ \left(\Delta_q h\right)_{v_i}=\frac{1}{2} \sum_{\substack{j \mid(i, j)
+ \in E \\ \text { or }(j, i) \in E}}\left(\cot \left(\alpha_{i j}
+ \right)+\cot \left(\beta_{i j}\right)\right)\left(h_i-h_j\right)
+
Parameters ---------- tangent_vec : array-like, shape=[..., n_vertices, 3]
@@ -941,7 +963,7 @@
Source code for geomstats.geometry.discrete_surfaces
"""Elastic metric defined by a family of second order Sobolev metrics. Each individual discrete surface is represented by a 2D-array of shape
- `[n_vertices, 3]`. See [HSKCB2022]_ for details.
+ `[n_vertices, 3]`. See [HSKCB2022]_ and [HPBDC2023]_ (appendix) for details. The parameters a0, a1, b1, c1, d1, a2 (detailed below) are non-negative weighting coefficients for the different terms in the metric.
@@ -971,9 +993,14 @@
Source code for geomstats.geometry.discrete_surfaces
References
----------
- .. [HSKCB2022] "Elastic shape analysis of surfaces with second-order
- Sobolev metrics: a comprehensive numerical framework".
- arXiv:2204.04238 [cs.CV], 25 Sep 2022
+ .. [HSKCB2022] Emmanuel Hartman, Yashil Sukurdeep, Eric Klassen,
+ Nicolas Charon, and Martin Bauer.
+ "Elastic shape analysis of surfaces with second-order Sobolev metrics:
+ a comprehensive numerical framework". arXiv:2204.04238 [cs.CV], 25 Sep 2022
+ .. [HPBDC2023] Emmanuel Hartman, Emery Pierson, Martin Bauer,
+ Mohamed Daoudi, and Nicolas Charon.
+ “Basis Restricted Elastic Shape Analysis on the Space of Unregistered Surfaces.”
+ arXiv, November 7, 2023. https://doi.org/10.48550/arXiv.2311.04382. """def__init__(self,space,a0=1.0,a1=1.0,b1=1.0,c1=1.0,d1=1.0,a2=1.0):
@@ -985,26 +1012,27 @@
Source code for geomstats.geometry.discrete_surfaces
self.d1=d1self.a2=a2
- self.exp_solver=DiscreteSurfacesExpSolver(space,n_steps=10)
+ ifnotgs.__name__.endswith("numpy"):
+ self.exp_solver=DiscreteSurfacesExpSolver(space,n_steps=10)
- optimizer=ScipyMinimize(
- method="L-BFGS-B",
- jac="autodiff",
- options={"disp":False,"ftol":0.001},
- )
- self.log_solver=PathStraightening(space,n_nodes=10,optimizer=optimizer)
+ optimizer=ScipyMinimize(
+ method="L-BFGS-B",
+ jac="autodiff",
+ options={"disp":False,"ftol":0.001},
+ )
+ self.log_solver=PathStraightening(space,n_nodes=10,optimizer=optimizer)def_inner_product_a0(self,tangent_vec_a,tangent_vec_b,vertex_areas_bp):r"""Compute term of order 0 within the inner-product. Denote h and k the tangent vectors a and b respectively.
- Denote q the base point, i.e. the surface.
+ Denote q the base point.
+
+ This method computes :math:`G_{a_0} = a_0 <h, k>`, i.e.
- The equation of the inner-product is:
- :math:`\int_M (G_{a_0} + G_{a_1} + G_{b_1} + G_{c_1} + G_{d_1} + G_{a_2})vol_q`.
+ .. math::
- This method computes :math:`G_{a_0} = a_0 <h, k>`,
- with notations taken from [HSKCB2022]_.
+ G_{a_0}(h, h) = \sum_{i=1}^N\left\|h_i\right\|^2 \operatorname{vol}_{x_i} Parameters ----------
@@ -1017,18 +1045,11 @@
Source code for geomstats.geometry.discrete_surfaces
Returns
-------
- _ : array-like, shape=[...,]
+ inner_prod_a0 : array-like, shape=[...,] Term of order 0, and coefficient a0, of the inner-product.
-
- References
- ----------
- .. [HSKCB2022] "Elastic shape analysis of surfaces with second-order
- Sobolev metrics: a comprehensive numerical framework".
- arXiv:2204.04238 [cs.CV], 25 Sep 2022. """returnself.a0*gs.sum(
- vertex_areas_bp
- *gs.einsum("...bi,...bi->...b",tangent_vec_a,tangent_vec_b),
+ gs.dot(tangent_vec_a,tangent_vec_b)*vertex_areas_bp,axis=-1,)
@@ -1036,38 +1057,33 @@
Source code for geomstats.geometry.discrete_surfaces
r"""Compute a1 term of order 1 within the inner-product. Denote h and k the tangent vectors a and b respectively.
- Denote q the base point, i.e. the surface.
+ Denote q the base point.
+
+ This method computes :math:`G_{a_1} = a_1.g_q^{-1} <dh_m, dk_m>`, i.e.
- The equation of the inner-product is:
- :math:`\int_M (G_{a_0} + G_{a_1} + G_{b_1} + G_{c_1} + G_{d_1} + G_{a_2})vol_q`.
+ .. math::
- This method computes :math:`G_{a_1} = a_1.g_q^{-1} <dh_m, dk_m>`,
- with notations taken from [HSKCB2022]_.
+ G_{a_1}(h, h) = \sum_{f \in F} \operatorname{tr}\left(g_f^{-1}
+ \delta g_f g_f^{-1} \delta g_f\right) \operatorname{vol}_f Parameters ----------
- ginvdga : array-like, shape=[n_faces, 2, 2]
+ ginvdga : array-like, shape=[..., n_faces, 2, 2] Product of the inverse of the surface metric matrices with their differential at a.
- ginvdgb : array-like, shape=[n_faces, 2, 2]
+ ginvdgb : array-like, shape=[..., n_faces, 2, 2] Product of the inverse of the surface metric matrices with their differential at b.
- areas_bp : array-like, shape=[n_faces,]
+ areas_bp : array-like, shape=[..., n_faces,] Areas of the faces of the surface given by the base point. Returns -------
- _ : array-like, shape=[...,]
- Term of order 0, and coefficient a1, of the inner-product.
-
- References
- ----------
- .. [HSKCB2022] "Elastic shape analysis of surfaces with second-order
- Sobolev metrics: a comprehensive numerical framework".
- arXiv:2204.04238 [cs.CV], 25 Sep 2022.
+ inner_prod_a1 : array-like, shape=[...,]
+ Term of order 1, and coefficient a1, of the inner-product. """returnself.a1*gs.sum(
- gs.einsum("...bii->...b",gs.matmul(ginvdga,ginvdgb))*areas_bp,
+ gs.trace(gs.matmul(ginvdga,ginvdgb))*areas_bp,axis=-1,)
@@ -1075,40 +1091,33 @@
Source code for geomstats.geometry.discrete_surfaces
r"""Compute b1 term of order 1 within the inner-product. Denote h and k the tangent vectors a and b respectively.
- Denote q the base point, i.e. the surface.
+ Denote q the base point.
- The equation of the inner-product is:
- :math:`\int_M (G_{a_0} + G_{a_1} + G_{b_1} + G_{c_1} + G_{d_1} + G_{a_2})vol_q`.
+ This method computes :math:`G_{b_1} = b_1.g_q^{-1} <dh_+, dk_+>`, i.e.
- This method computes :math:`G_{b_1} = b_1.g_q^{-1} <dh_+, dk_+>`,
- with notations taken from [HSKCB2022]_.
+ .. math::
+
+ G_{b_1}(h, h) = \sum_{f \in F} \operatorname{tr}
+ \left(g_f^{-1} \delta g_f\right)^2 \operatorname{vol}_f Parameters ----------
- ginvdga : array-like, shape=[n_faces, 2, 2]
+ ginvdga : array-like, shape=[..., n_faces, 2, 2] Product of the inverse of the surface metric matrices with their differential at a.
- ginvdgb : array-like, shape=[n_faces, 2, 2]
+ ginvdgb : array-like, shape=[..., n_faces, 2, 2] Product of the inverse of the surface metric matrices with their differential at b.
- areas_bp : array-like, shape=[n_faces,]
+ areas_bp : array-like, shape=[..., n_faces,] Areas of the faces of the surface given by the base point. Returns -------
- _ : array-like, shape=[...,]
- Term of order 0, and coefficient b1, of the inner-product.
-
- References
- ----------
- .. [HSKCB2022] "Elastic shape analysis of surfaces with second-order
- Sobolev metrics: a comprehensive numerical framework".
- arXiv:2204.04238 [cs.CV], 25 Sep 2022.
+ inner_prod_b1 : array-like, shape=[...,]
+ Term of order 1, and coefficient b1, of the inner-product. """returnself.b1*gs.sum(
- gs.einsum("...bii->...b",ginvdga)
- *gs.einsum("...bii->...b",ginvdgb)
- *areas_bp,
+ gs.trace(ginvdga)*gs.trace(ginvdgb)*areas_bp,axis=-1,)
@@ -1116,13 +1125,15 @@
Source code for geomstats.geometry.discrete_surfaces
r"""Compute c1 term of order 1 within the inner-product. Denote h and k the tangent vectors a and b respectively.
- Denote q the base point, i.e. the surface.
-
- The equation of the inner-product is:
- :math:`\int_M (G_{a_0} + G_{a_1} + G_{b_1} + G_{c_1} + G_{d_1} + G_{a_2})vol_q`.
+ Denote q the base point. This method computes :math:`G_{c_1} = c_1.g_q^{-1} <dh_\perp, dk_\perp>`,
- with notations taken from [HSKCB2022]_.
+ i.e.
+
+ .. math::
+
+ G_{c_1}(h, h) = \sum_{f \in F}\left\langle\delta n_f,
+ \delta n_f\right\rangle \mathrm{vol}_f Parameters ----------
@@ -1130,94 +1141,80 @@
Source code for geomstats.geometry.discrete_surfaces
Point a corresponding to tangent vec a.
point_b : array-like, shape=[..., n_vertices, 3] Point b corresponding to tangent vec b.
- normals_bp : array-like, shape=[n_faces, 3]
+ normals_bp : array-like, shape=[..., n_faces, 3] Normals of each face of the surface given by the base point.
- areas_bp : array-like, shape=[n_faces,]
+ areas_bp : array-like, shape=[..., n_faces,] Areas of the faces of the surface given by the base point. Returns -------
- _ : array-like, shape=[...,]
- Term of order 0, and coefficient c1, of the inner-product.
-
- References
- ----------
- .. [HSKCB2022] "Elastic shape analysis of surfaces with second-order
- Sobolev metrics: a comprehensive numerical framework".
- arXiv:2204.04238 [cs.CV], 25 Sep 2022.
+ inner_prod_c1 : array-like, shape=[...,]
+ Term of order 1, and coefficient c1, of the inner-product. """dna=self._space.normals(point_a)-normals_bpdnb=self._space.normals(point_b)-normals_bp
- returnself.c1*gs.sum(
- gs.einsum("...bi,...bi->...b",dna,dnb)*areas_bp,axis=-1
- )
+
+ returnself.c1*gs.sum(gs.dot(dna,dnb)*areas_bp,axis=-1)def_inner_product_d1(
- self,one_forms_a,one_forms_b,one_forms_bp,areas_bp,inv_surface_metrics_bp
+ self,one_forms_a,one_forms_b,one_forms_bp,areas_bp,ginv_bp):r"""Compute d1 term of order 1 within the inner-product. Denote h and k the tangent vectors a and b respectively.
- Denote q the base point, i.e. the surface.
+ Denote q the base point.
- The equation of the inner-product is:
- :math:`\int_M (G_{a_0} + G_{a_1} + G_{b_1} + G_{c_1} + G_{d_1} + G_{a_2})vol_q`.
+ This method computes :math:`G_{d_1} = d_1.g_q^{-1} <dh_0, dk_0>`, i.e.
- This method computes :math:`G_{d_1} = d_1.g_q^{-1} <dh_0, dk_0>`,
- with notations taken from [HSKCB2022]_.
+ .. math::
+
+ G_{d_1}(h, h) = \sum_{f \in F} \operatorname{tr}\left(g_f^{-1}
+ \xi_f g_f^{-1} \xi_f^T\right) \operatorname{vol}_f
+
+ where :math:`\xi_f=d q_f^T d h_f-d h_f^T d q_f`. Parameters ----------
- one_forms_a : array-like, shape=[n_points, n_faces, 2, 3]
+ one_forms_a : array-like, shape=[..., n_faces, 2, 3] One forms at point a corresponding to tangent vec a.
- one_forms_b : array-like, shape=[n_points, n_faces, 2, 3]
+ one_forms_b : array-like, shape=[..., n_faces, 2, 3] One forms at point b corresponding to tangent vec b.
- one_forms_bp : array-like, shape=[n_faces, 2, 3]
+ one_forms_bp : array-like, shape=[..., 2, 3] One forms at base point.
- areas_bp : array-like, shape=[n_faces,]
+ areas_bp : array-like, shape=[..., n_faces,] Areas of the faces of the surface given by the base point.
- inv_surface_metrics_bp : array-like, shape=[n_faces, 2, 2]
+ ginv_bp : array-like, shape=[..., n_faces, 2, 2] Inverses of the surface metric matrices at each face. Returns -------
- _ : array-like, shape=[...,]
- Term of order 0, and coefficient d1, of the inner-product.
-
- References
- ----------
- .. [HSKCB2022] "Elastic shape analysis of surfaces with second-order
- Sobolev metrics: a comprehensive numerical framework".
- arXiv:2204.04238 [cs.CV], 25 Sep 2022.
+ inner_prod_d1 : array-like, shape=[...,]
+ Term of order 1, and coefficient d1, of the inner-product. """one_forms_bp_t=Matrices.transpose(one_forms_bp)
- one_forms_a_t=Matrices.transpose(one_forms_a)
- xa=one_forms_a_t-one_forms_bp_t
+ aux=gs.matmul(one_forms_bp_t,ginv_bp)
+ xa=one_forms_a-one_forms_bpxa_0=gs.matmul(
- gs.matmul(one_forms_bp_t,inv_surface_metrics_bp),
- gs.matmul(Matrices.transpose(xa),one_forms_bp_t)
- -gs.matmul(one_forms_bp,xa),
+ aux,
+ gs.matmul(xa,one_forms_bp_t)
+ -gs.matmul(one_forms_bp,Matrices.transpose(xa)),)
- one_forms_b_t=Matrices.transpose(one_forms_b)
- xb=one_forms_b_t-one_forms_bp_t
+ xb=one_forms_b-one_forms_bpxb_0=gs.matmul(
- gs.matmul(one_forms_bp_t,inv_surface_metrics_bp),
- gs.matmul(Matrices.transpose(xb),one_forms_bp_t)
- -gs.matmul(one_forms_bp,xb),
+ aux,
+ gs.matmul(xb,one_forms_bp_t)
+ -gs.matmul(one_forms_bp,Matrices.transpose(xb)),)returnself.d1*gs.sum(
- gs.einsum(
- "...bii->...b",
- gs.matmul(
+ gs.trace(
+ Matrices.mul(xa_0,
- gs.matmul(
- inv_surface_metrics_bp,
- Matrices.transpose(xb_0),
- ),
+ ginv_bp,
+ Matrices.transpose(xb_0),),)*areas_bp,
@@ -1230,13 +1227,15 @@
Source code for geomstats.geometry.discrete_surfaces
r"""Compute term of order 2 within the inner-product. Denote h and k the tangent vectors a and b respectively.
- Denote q the base point, i.e. the surface.
+ Denote q the base point.
+
+ This method computes :math:`G_{a_2} = a_2 <\Delta_q h, \Delta_q k>`, i.e.
- The equation of the inner-product is:
- :math:`\int_M (G_{a_0} + G_{a_1} + G_{b_1} + G_{c_1} + G_{d_1} + G_{a_2})vol_q`.
+ .. math::
+
+ G_{a_2}(h, h) = \sum_{i=1}^N\left\|\left(\Delta_q h\right)_{v_i}
+ \right\|^2 \operatorname{vol}_{x_i}
- This method computes :math:`G_{a_2} = a_2 <\Delta_q h, \Delta_q k>`,
- with notations taken from [HSKCB2022]_. Parameters ----------
@@ -1244,26 +1243,19 @@
Source code for geomstats.geometry.discrete_surfaces
Tangent vector at base point.
tangent_vec_b : array-like, shape=[..., n_vertices, 3] Tangent vector at base point.
- base_point : array-like, shape=[n_vertices, 3]
+ base_point : array-like, shape=[..., n_vertices, 3] Base point, a surface i.e. the 3D coordinates of its vertices.
- vertex_areas_bp : array-like, shape=[n_vertices, 1]
+ vertex_areas_bp : array-like, shape=[..., n_vertices, 1] Vertex areas for each vertex of the base_point. Returns -------
- _ : array-like, shape=[...,]
+ inner_prod_a2 : array-like, shape=[...,] Term of order 2, and coefficient a2, of the inner-product.
-
- References
- ----------
- .. [HSKCB2022] "Elastic shape analysis of surfaces with second-order
- Sobolev metrics: a comprehensive numerical framework".
- arXiv:2204.04238 [cs.CV], 25 Sep 2022. """laplacian_at_base_point=self._space.laplacian(base_point)returnself.a2*gs.sum(
- gs.einsum(
- "...bi,...bi->...b",
+ gs.dot(laplacian_at_base_point(tangent_vec_a),laplacian_at_base_point(tangent_vec_b),)
@@ -1283,7 +1275,10 @@
Source code for geomstats.geometry.discrete_surfaces
We denote q the base point, i.e. the surface.
The six terms of the inner-product are given by:
- :math:`\int_M (G_{a_0} + G_{a_1} + G_{b_1} + G_{c_1} + G_{d_1} + G_{a_2})vol_q`
+
+ .. math::
+
+ \int_M (G_{a_0} + G_{a_1} + G_{b_1} + G_{c_1} + G_{d_1} + G_{a_2})vol_q where:
@@ -1294,8 +1289,6 @@
Source code for geomstats.geometry.discrete_surfaces
Emmanuel Hartman, Yashil Sukurdeep, Eric Klassen,
+Nicolas Charon, and Martin Bauer.
+“Elastic shape analysis of surfaces with second-order Sobolev metrics:
+a comprehensive numerical framework”. arXiv:2204.04238 [cs.CV], 25 Sep 2022
+
+
+[HPBDC2023]
+
Emmanuel Hartman, Emery Pierson, Martin Bauer,
+Mohamed Daoudi, and Nicolas Charon.
+“Basis Restricted Elastic Shape Analysis on the Space of Unregistered Surfaces.”
+arXiv, November 7, 2023. https://doi.org/10.48550/arXiv.2311.04382.
Marshall, Albert W., and Ingram Olkin.
“Scaling of Matrices to Achieve Specified Row and Column Sums.”
Numerische Mathematik 12, no. 1 (August 1, 1968): 83–90.
https://doi.org/10.1007/BF02170999.
Johnson, Charles R., and Robert Reams.
“Scaling of Symmetric Matrices by Positive Diagonal Congruence.”
Linear and Multilinear Algebra 57, no. 2 (March 1, 2009): 123–40.
@@ -5917,7 +5950,7 @@
Following [Chikuse03], \(n\_samples * n * p\) scalars are sampled
+
Following [Chikuse03], \(n\_samples * n * p\) scalars are sampled
from a standard normal distribution and reshaped to matrices,
the projectors on their first \(p\) columns follow a uniform
distribution.
Following [Chikuse03], \(n\_samples * n * p\) scalars are sampled
+
Following [Chikuse03], \(n\_samples * n * p\) scalars are sampled
from a standard normal distribution and reshaped to matrices,
the projectors on their first \(p\) columns follow a uniform
distribution.
Class for bi-invariant metrics on compact Lie groups.
Compact Lie groups and direct products of compact Lie groups with vector
-spaces admit bi-invariant metrics [Gallier]. Products Lie groups are not
+spaces admit bi-invariant metrics [Gallier]. Products Lie groups are not
implemented. Other groups such as SE(3) admit bi-invariant pseudo-metrics.
Kolev, Boris. “Lie Groups and Mechanics: An Introduction.”
Journal of Nonlinear Mathematical Physics 11, no. 4, 2004:
480–98. https://doi.org/10.2991/jnmp.2004.11.4.5.
Calculate the Baker-Campbell-Hausdorff approximation of given order.
-
The implementation is based on [CM2009a] with the pre-computed
-constants taken from [CM2009b]. Our coefficients are truncated to
+
The implementation is based on [CM2009a] with the pre-computed
+constants taken from [CM2009b]. Our coefficients are truncated to
enable us to calculate BCH up to order 15.
This represents Z = log(exp(X)exp(Y)) as an infinite linear combination
of the form Z = sum z_i e_i where z_i are rational numbers and e_i are
@@ -8519,13 +8552,13 @@
F. Casas and A. Murua. An efficient algorithm for
computing the Baker–Campbell–Hausdorff series and some of its
applications. Journal of Mathematical Physics 50, 2009