Biological morphometry historically relies on manually selected homologous landmarks, defined as anatomically unambiguous and consistent features of the object of interest (Bookstein, 1997). Relying on these landmark points provides an implicit form of sparse correspondence between objects, which then allows aligning collections of specimen via Procrustes analysis, a classical strategy consisting of removing geometrical transformations that do not affect the shape of an object (specifically, translations, scaling and rotations) so as to statistically study the extent and nature of shape differences (Rohlf and Slice, 1990).

Landmark-based morphometry, also called geometric morphometrics (Dryden and Mardia, 2016), is actively used at the macroscopic scale in the context of medical imaging (Yeh et al., 2021), anatomy (Finka et al., 2019), taxonomy (Karanovic and Bláha, 2019), and plant science (Lucas et al., 2013). However, when considered at the smaller mesoscopic or microscopic scales, biological objects such as soft tissues and cells in isolation rarely possess well-localized and unambiguously-identifiable features that could reliably act as landmarks, despite having non-trivial shapes. In geometric morphometrics, the fact that landmarks are homologous and that they correspond to the same biological structure or function across different individuals is relevant (Klingenberg, 2008). For the non-rigid, featureless objects most often encountered in microscopy, a more appropriate alternative consists of identifying a dense (continuous) correspondence that only uses geometry and does not rely on the existence of a finite set of unique, localized features derived from functional or evolutionary factors. This type of approach is successfully used in several frameworks in 2D (Laga et al., 2014; Phillip et al., 2021). Extending these methods to 3D shapes is unfortunately challenging when at all possible, and existing solutions impose strong constraints on the topology of input objects and on the way they are parameterized (Srivastava et al., 2010; Koehl and Hass, 2015). As a result, most 3D shape correspondence workflow involving objects extracted from bioimages still rely on manually annotated (pseudo-) landmarks (Boehm et al., 2011; White et al., 2019).

Several algorithms that automatically retrieve a dense correspondence between 3D surfaces have been proposed by the computational geometry and computer graphics communities. Herein it is generally assumed that the shapes to correspond are near-isometric, meaning that they do not exhibit significant deformations. This assumption alone unfortunately makes the vast majority of solutions unsuitable to problems involving collections of biological objects, where natural individual-to-individual variability may fully deviate from isometry. Most popular approaches retrieve a correspondence map that minimizes a measure of distortion, such as the degree of stretching or bending, using continuous optimisation techniques (Schreiner et al., 2004; Sorkine and Alexa, 2007; Schmidt et al., 2019). These methods however suffer from highly non-convex energy landscapes composed of many sub-optimal local minima, and thus require a relatively good initialization to converge to a good solution. This problem is circumvented by Windheuser et al. (2011), in which a global solution is found using a sophisticated technique to reduce the solution space and then solve the optimization problem using a linear programming approach. The latter approach is guaranteed to yield the lowest distortion mapping and is therefore less likely to suffer from symmetry problems and possible mis-assignments provided that an appropriate distortion function is used. Each matching however then requires hours of processing time even with GPU speedup, making it poorly scalable to large bioimage datasets.

Another class of methods approaches the problem by mapping each of the objects to correspond into intermediate domains and then recovers a mapping between these domains. One such example exploiting the relatively easier task of finding mapping between intermediate domains is provided by Lipman and Funkhouser (2009), available at https://github.com/pedrofreire/shape-matching. There, a conformal parameterisation (i.e., an angle-preserving mapping of the objects to the plane ℝ²) is used to flip the problem into searching for conformal automorphisms of the plane using Möbius transformations. This method, however, requires good landmarks for initialization and is restricted to shapes of specific topology. Finally, Ovsjanikov et al. (2012) take a more abstract approach by using the scalar function space associated with each object as intermediate domain, in a strategy referred to as functional mapping. The function space is defined as the set of all functions from the object's surface, represented as a mesh M, to the real numbers, FM = {f:M → ℝ}. This last technique is particularly interesting for biological shape correspondence for a number of reasons: firstly, it reduces 3D correspondence to a linear problem that does not require initialization; secondly, it is relatively robust to non-isometry and flexible compared to other continuous approaches as it imposes few constraints on the input shapes; and finally, it provides a computationally efficient means of determining correspondences. For these reasons, we chose to make use of this approach in the μMatch pipeline.

The above literature overview, and the design of the μMatch pipeline, focuses solely on surface-based alignment. An alternative could have been to review and rely on volume-based alignment methods. Our rationale for choosing a surface-based strategy instead is as follows. Volume-based alignment approaches, as implemented for instance in the popular elastix software (Klein et al., 2009) and commonly used in medical imaging, are unable to handle objects with significantly different orientations. Images must then be pre-aligned or a close-enough initialization must be provided (Miao et al., 2016; Yang et al., 2017). While medical image data often exhibit small enough variations in object orientation (due to the acquisition protocols) and small enough sample variation (due to the nature of the objects being imaged) for this issue to be addressed with ad-hoc methods, the extent of sample variation and orientation difference in biological experiments cannot be known a priori. In contrast, the performance of surface-based approaches is not affected by the orientation difference and, as a consequence, they do not require any pre-alignment step. A further issue with biological data lies in the anisotropic nature of microscopy image volumes, making additional steps of interpolation mandatory when relying on voxels to align and further complicating the task of volume-based strategies in the case of large orientation differences. Surface-based methods, as they rely on meshes, have the advantage of being blind to the anisotropy of the input data. In addition to these technical aspects, most alignment algorithms for bioimage analysis originate from medical imaging in general and neuroimaging in particular, in which surface-based alignment is now accepted to be superior to volume-based alignment. Surface-based matching has indeed been shown to map borders more accurately between brains than volume-based registration (Brodoehl et al., 2020), further strengthening the case for a surface-based approach. Finally, since biological images most often capture objects that purely have a surface signal (e.g., from membrane stains) and rarely have a conserved inner structure akin to human organs, the information available for alignment is overwhelmingly held in the object's surface, making it all the more relevant to rely on a surface-based method.

Through the article, we focus on the problem of establishing correspondence between 3D biological objects represented as meshes. We in particular do not study the downstream problems of object segmentation and voxel set meshing. Obtaining object segmentation from possibly noisy biological data can be challenging but is beyond the scope of this paper. This problem has been extensively studied and several robust state-of-the-art methods relying on deep learning, such as the popular 3D U-net (Isensee et al., 2021), are now freely available in user-friendly softwares (Lucas et al., 2021). Once segmented, voxel sets can reliably be meshed using classical methods such as the popular marching cubes algorithm (Lorensen and Cline, 1987), provided in widely-used 3D image processing Python libraries such as scikit-image (van der Walt et al., 2014).

Triangular meshes are discrete approximations of 3D surfaces commonly used as representations for computer graphics and computational geometry. A mesh M is composed of vertices, edges, and (triangular) faces to form a continuous but not smooth surface, due to the presence of sharp edges connecting any two neighboring faces. A mesh is generally defined by two arrays of numbers, the first one containing the spatial positions of the vertices, v∈ℝnv×3, and the second one the mesh faces, f∈ℕnf×3. Increasing the number of vertices of a mesh provides increased resolution and smoothness, at the expense of increased complexity and memory requirements.

Numerous classical Euclidean operators can be extended to manifolds, and as a consequence to discrete meshes, including the gradient ∇ and the Laplacian Δ = ∇² = ∇·∇ operators. The Laplacian is interesting in particular because its eigenfunctions (i.e., the functions that are only scaled by the action of the operator) provide a geometrically informative basis for the scalar space of the manifold. The Laplace-Beltrami operator is the generalization of the Laplacian to triangular meshes. For a mesh of n_v vertices, the Laplace-Beltrami operator is represented as a sparse n_v × n_v matrix. Because a mesh is not differentiable, such a discretisation requires care in order to reproduce the expected behavior of the Laplacian: the most common formalism for doing is the cotangent Laplace-Beltrami operator (Pinkall and Polthier, 1993), given by

where α_ij is the angle adjacent to the edge ij, N(i) the connected neighbors of the i^th vertex, and cot the cotangent operator.

The shortest distance between any two points on a mesh is given by a geodesic path, as illustrated in Figure 1. Formally, given a triangular mesh M with associated metric g, the geodesic path x parameterized by t ∈ ℝ between two points a and b on the mesh is defined as

where μand ν correspond to coordinates on the mesh.

Two meshes are equivalent if there is a map between them that preserves the metrics, and therefore the geodesic distances, on them. The degree of deformation induced by a map can thus be measured by calculating the extent to which geodesic distances are altered by the map, captured in the geodesic matrix GM. The geodesic matrix is constructed by computing all pairwise geodesic distances such that [GM]_ij contains the geodesic distance between the two vertices i and j on the mesh M.

A symmetry of a surface is defined as a self-mapping (i.e., automorphism) Ψ of the mesh that leaves the geodesic matrix unchanged, which is formally expressed as Ψ:M→M such that [GM]_x,y = [GM]_{Ψ(x), Ψ(y)} for every x, y ∈ M.

Finally, a correspondence Φ between two meshes is formally defined as a mapping that assigns each vertex on a first mesh M1 to each vertex on a second mesh M2 as Φ:M1 → M2. A functional mapping T:FM1 → FM2 maps the space of scalar functions on one mesh (FM1 ) to the scalar functional space on the other (FM2). As will be extensively discussed through the paper, the functional mapping provides a useful way of representing the correspondence between two meshes. In order to use it, we need a concrete representation in the form of a finite (k × k) correspondence matrix, denoted as C_Φ.

To make them amenable to processing with μMatch, the input objects must satisfy a number of technical requirements. First, each object's surface must be represented in the form of a triangular mesh. These meshes must additionally be both connected and manifold. Being connected implies that, given any two vertices in the mesh, it is possible to find a path consisting of a subset of the mesh edges (i.e., the sides of the faces) that links the two. The manifold condition imposes that the mesh must represent a physically realizable continuous surface. More specifically, it means that each edge in the mesh should be incident to at most two faces, and that the faces attached to a vertex should form an open disk or half-disk around the vertex. An example of non-manifold mesh may involve two parts that are connected by a single vertex, or the presence of an internal face. These different notions are illustrated in Figure 3.

In addition, the matching process assumes that the correspondence map is a bijection, meaning that it defines a way to match M1 to M2 as much as a way to match M2 to M1. It is therefore crucial that the pairs of meshes have the same coverage, which implies that for each point in M1, a corresponding point exists in M2 and vice versa. Finally, it is important for all meshes in the collection to be constructed by sampling as uniformly as possible the object's surface. This translates to vertices in the mesh being equally spaced on each surface, or equivalently to triangular faces having a constant area. μMatch makes use of a number of geometric quantities associated with each mesh, including the geodesic matrix and the Laplace-Beltrami eigen-decomposition introduced in Section 2.2, that may be adversely affected by a non-uniform sampling.

In μMatch, input meshes are cleaned using PyMeshFix 0.14.1 (Attene, 2010), which offers built-in functionalities to remove common errors in triangular meshes, including degenerate and intersecting faces. As a result, output meshes are manifold and watertight with a single connected component. Meshes are subsequently resampled to a user-defined number of vertices N that is smaller or equal than the total number of vertices in the smallest of the input meshes. Importantly, this resampling step processes all meshes to have the same number of vertices regardless of the original number of vertices they were composed of, relieving users from having to handle this constraint. The value of N should be chosen so as to aim at having fine enough meshes to capture important shape features, while limiting the number of vertices to what is necessary in order to avoid overburdening the matching process later.

A final preprocessing step in μMatch, consists of calculating the geodesic matrix of each input mesh using Lib-igl (Jacobson and Panozzo, 2018). Although Lib-igl implements the fast exact geodesic algorithm (Mitchell et al., 1985), it is still computationally demanding with a run-time of O(N²logN) and memory requirement of O(N²), with N the number of vertices in the mesh. For the sake of computational efficiency, we therefore provide a custom data preparation script that precomputes the geodesic matrices for the whole collection of objects and save them to disk prior to the correspondence pipeline.

The core of μMatch is a landmark-free dense correspondence algorithm adapted from Ovsjanikov et al. (2012), Litany et al. (2017), and Halimi et al. (2019). The correspondence matrix, describing the final mapping, is built from a collection of feature descriptors and refined through filtering steps. The overall workflow is summarized in Figure 4.

In order to quantitatively validate μMatch, we use a teeth scan dataset originally introduced by Boyer et al. (2011) and available at www.wisdom.weizmann.ac.il/~ylipman/CPsurfcomp/ and commonly exploited to assess 3D morphometry algorithms. The dataset includes 45 meshes corresponding to mandibular second molars of a variety of prosimian primates and their close non-primate evolutionary relatives, summarized in Table 2. Each specimen is annotated with a set of 16 manually-determined landmarks placed by expert morphometricians that can be used as ground-truth to assess the quality of an automatically-retrieved correspondence. Although this dataset is composed of objects toward the macroscale of the spectrum, it is one of the rare available resources of biological objects for which a ground truth sparse correspondence is available and is as such an excellent candidate for benchmarking. Performance assessment is then carried out following the Princeton benchmark protocol (Shilane et al., 2004; Kim et al., 2011) by computing the mean geodesic error (GE) expressed as

where {li1}i=116 and {li2}i=116 are the ground-truth landmarks on M1 and M2, respectively, Φ is the point-to-point mapping, and GM2 and AM2 are the geodesic matrix and the surface area of the mesh M2, respectively. The GE therefore measures how well the main features of the objects have been mapped, regardless of the amount of deformation there may be elsewhere on the surface.

Several processing steps of the μMatch correspondence module (Figure 4), namely the deep functional maps and product manifold filter, aim at improving the final correspondence but are not strictly necessary. While we do recommend including them whenever possible, these steps do increase computation time and the extent of the improvement they will bring cannot be predicted in general for any arbitrary dataset. To provide users with an example of the cost and gains involved, we have quantified the quality of the retrieved correspondence in different settings, by enabling and disabling the deep functional maps and product manifold filter in the μMatch pipeline. In Table 1, we report the average GE between ground-truth landmarks and correspondence points obtained automatically with μMatch for the whole dataset. We also indicate the average runtime to establish correspondence in any given pair. A breakdown of the average GE by species and the cumulative error curves are provided in Figure 7A and Table 2, respectively. We observe the respective merits of the deep functional maps and product manifold filter to be as follows. The deep functional maps improve the signature functions but do not impose smoothness in the final correspondence, resulting in large GE variance. While the correspondence runtime is not affected by deep functional maps since the network is trained during preprocessing, a significant extra amount of compute (~3 h for this dataset) is required ahead of the matching itself, impacting the total duration of the pipeline. In contrast, the product manifold filter ensures that the final correspondence is smooth and does not feature major tears or discontinuities, resulting in small GE variance. It is limited by the correspondence quality that can be obtained from the original signature functions alone and doubles the runtime per pair, but as a trade-off requires no prior computation or training. The combination of the two yields the best results, as it allows obtaining a smooth correspondence from improved signature functions. The cost is however a high runtime per pair for the correspondence itself, and the necessity to go through an expensive training stage during preprocessing.

Once correspondence has been established between all meshes in the collection, we carry out a further sanity check by reporting the geodesic distance, or deformation, given by

where E is the set of edges of the meshes M1 and M2, and observing how it compares between the different species present in the dataset. We hypothesize that samples from different species will exhibit more shape deviations than individuals from the same species, and investigate whether this can be observed in the geodesic distances retrieved after correspondence established by μMatch. The quantity (Equation 19) is calculated between all pairs in the collection, producing a n × n matrix δG, with n the number of meshes in the collection. Each sub-blocks of the matrix corresponding to the different species can be averaged to produce a reduced distance matrix, depicted in Figure 7B for the full μMatch pipeline (with deep functional maps and product manifold filter). The geodesic distances between species match what would be expected from their phylogeny: molar surfaces of chimpanzees and bonobos are highly similar, and also closely resemble human teeth, while gorilla and orangutan molars are observed to be morphologically more dissimilar. In addition to their differences with other species, samples from the orangutan class appear to be extremely diverse, resulting in a large intra-class distance. This experiment provides a sanity check assessing the validity of the algorithms implemented in the μMatch pipeline. We orient the reader interested in a comparative assessment of the relative performance of the selected correspondence algorithms against published alternatives to the original works introducing functional mapping (Ovsjanikov et al., 2012, 2016).

In our previous study (Martínez-Abadías et al., 2018), we performed geometric morphogenesis on a mouse model of Apert syndrome, in which the Fgfr2 gene contains a mutation which models the human syndrome. The goal was to explore whether shape analysis of a gene expression pattern could suggest the molecular basis of the phenotype. In addition to analyzing the developing anatomical changes in the limb, we also analyzed the 3D expression pattern of Dusp6, a gene whose expression reflects the activity of FGF signaling. This allowed subtle alterations in the gene expression pattern to be detected before the anatomical phenotype was apparent, thus strengthening the idea that altered FGF signaling is directly responsible for the phenotype. That previous analysis however depended on manual annotation of 3D landmarks on the data, a process which was very labor-intensive and taking many weeks to achieve. Here, we demonstrate that μMatch can be used to automate this process, jointly quantify the evolution of the shape of embryonic limb and gene expression pattern over developmental time without human bias and achieving similar results in a fraction of the time. We have chosen to reproduce the results of this study to illustrate the use of μMatch on real biological data containing objects that differ in size, orientation, and shape, and exhibit a mixture of inter-group and intra-group variations of morphology.

The dataset we refer to as MOUSE_LIMB, publicly available at dx.doi.org/10.5061/dryad.8h646s0, contains meshes extracted from real optical projection tomography scans (Sharpe et al., 2002) of the limbs and of the Dusp6 gene expression domains of Apert syndrome (Fgfr2 mutants) mouse embryos and unaffected littermates ranging from E10 to E11.5, where EN stands for N days after conception. Extracting these meshes requires the raw image data to be first segmented, a step that can be carried out by different methods. We refer readers interested in the details of the segmentation process to the original study (Martínez-Abadías et al., 2018). The meshes are grouped into forelimbs and hind limbs of early (~E10), mid (~E10.5 to ~E11) and late (~E11.5) developmental periods, according to Table 3. To process these data, we set the number of vertices N to be 2000, and use default μMatch values, namely a functional space dimension of k = 100, and 100 computed heat kernel and wave kernel signatures.

For each age and genetic background group, a template shape is chosen at random and used to initiate correspondences with the remaining samples. From the mesh correspondences obtained with μMatch, a Procrustes alignment of the limbs surfaces is obtained following the procedure described in Section 3.3.1, yielding an average limb with an associated measure of variance at each vertex. Once limb surfaces are matched, the correspondence can be further used to map internal processes between them, allowing in particular to compare Dusp6 gene expression data available for some of the limbs. The parameters obtained from the limb alignment can then be used to do a rigid registration of the gene expression patterns. This is then further refined by iterative closest point alignment (Besl and McKay, 1992; Chen and Medioni, 1992) and thin plate spline (TPS) deformation to finely warp the objects onto one another (Duchon, 1977). To implement the TPS deformation, the surfaces are first sub-sampled to obtain M control points {cm}m=1M and the warping is obtained as

where ψ(r) = r²log(r) is the thin plate spline kernel. The coefficients {ωm}m=1M are found by solving for Aω = v, where [A]_ij = ψ(||c_i−c_j||)∀i, j = 1, …, M and [v]_m is the difference vector between the source and target surface at the control points m = 1, …, M. Once the objects have been warped onto one another in this manner, a simple nearest neighbors search is used to obtain a correspondence and the shape analysis procedure described in Section 3.3.1 can be carried out again for the gene expression data. The resulting average limbs and gene expression patterns, color-coded according to local variance, are shown in Figure 8. Following Martínez-Abadías et al. (2018), we display results for hind limbs of the early and mid-early stages, and fore limbs of the mid-late and late stages.

In Martínez-Abadías et al. (2018), the limb buds of Fgfr2 mutants from the late period (~E11.5) were observed to separate from those of unaffected littermates. During the mid-early (~E10.5) and early (~E10) periods, the limbs of Fgfr2 mutant mice were however undistinguishable from those of unaffected specimen. Significant difference in Dusp6 gene expression pattern was observed between Fgfr2 mutant animals and their unaffected littermates for all groups except for the early time period. We recover the same results using μMatch, as illustrated in Figure 9, demonstrating the validity of the automatically retrieved correspondence. The difference in shapes was originally quantified by Procrustes analysis relying on sparse set of semi-landmarks along the most distal edge and on the dorsal and ventral sides of the limb (Martínez-Abadías et al., 2018). In contrast, our approach solves for a dense correspondence, taking into account the entire object's surface. The major advantage of μMicro is that it is entirely landmark-free and fully automated. Once limb surfaces are put in correspondence, their geodesic distance can be used to further quantify the extent of their morphological difference, and the recovered matching used to also align their corresponding gene expression. As such, provided that gene expression patterns have been acquired so as to be appropriately registered inside the limb, comparing the shape of gene expression patterns does neither require artificial landmarking nor a separate correspondence map: they instead get aligned based on their “container.” These results demonstrate that μMatch can reliably be used on real bioimage data, which may be noisy, in which samples may exhibit too large variations in their morphology or orientation to rely on volume-based alignment, and which may involve extracted meshes of different numbers of vertices.