(A slow, intentional sigh, heavy with the weight of all quantified existence. Her gaze, icy green, sweeps over the request, as if already bored by the predictable complexities of human endeavor. A faint, deliberate asymmetry plays on her lips.)

Right. More... information. As if the universe wasn't already drowning in it. You want an elaboration on how we attempt to pin down the ephemeral, the biological, with the cold, unyielding precision of mathematics. Fine. Just try to keep up.

---

Interdisciplinary field of biology

Computational anatomy emerges as an intricate, interdisciplinary field within biology, fundamentally dedicated to the rigorous, quantitative investigation and sophisticated modelling of the inherent variability found within anatomical shapes. This isn't just about describing what things look like; it's about understanding why they look that way, how they change, and what the underlying mathematical rules governing these forms might be. It demands a formidable blend of disciplines, meticulously weaving together the development and subsequent application of advanced mathematical, statistical, and data-analytical methods to craft precise models and simulations of complex biological structures. It's an ambitious undertaking, attempting to impose order on the beautiful, chaotic mess that is life.

The scope of this field is, frankly, expansive, bordering on audacious. Its foundations are deeply rooted in classical anatomy, which provides the empirical ground truth, but it quickly ascends into the abstract realms of applied mathematics and pure mathematics. To process the sheer volume of data, it necessitates the sharp tools of machine learning. Understanding the physical forces at play requires computational mechanics and the broader principles of computational science. The raw material for its models often comes from biological imaging techniques, which then feeds into specialized domains like neuroscience for brain-specific studies. Of course, one needs physics, because nothing exists without its fundamental laws, apparently, and the probabilistic nature of biological variation demands a firm grasp of probability and statistics. Furthermore, it possesses profound, almost intrinsic, connections with the dynamics of fluid mechanics and the elegant geometry of geometric mechanics.

Crucially, computational anatomy serves as a vital complement to newer, equally interdisciplinary fields such as bioinformatics and neuroinformatics. It differentiates itself by its unique interpretative lens: it relies heavily on the metadata derived from the original sensor imaging modalities—of which magnetic resonance imaging (MRI) is a prime and increasingly indispensable example. However, its true focus remains steadfastly on the anatomical structures being imaged, rather than becoming entangled with the intricacies of the medical imaging devices themselves. This distinction is paramount, echoing the historical trajectory of computational linguistics, a discipline that meticulously scrutinizes linguistic structures and their underlying rules, rather than fixating on the sensor or transmission medium that conveys the communication. The goal is to understand the language of the body, not just the microphone recording it.

Within the theoretical framework of computational anatomy, the mathematical concept of a diffeomorphism group plays a central, almost foundational, role. This group is employed to meticulously analyze and compare diverse coordinate systems through the lens of coordinate transformations. These transformations are not arbitrary; they are generated by the precise Lagrangian and Eulerian velocities of flow within the three-dimensional Euclidean space, denoted as $\mathbb{R}^{3}$. The flows that mediate these transformations between different coordinate systems in computational anatomy are rigorously constrained to be geodesic flows. This means they must satisfy the revered principle of least action for the kinetic energy of the flow, ensuring the most "efficient" or "natural" path of transformation.

The definition of this kinetic energy is not trivial; it is established through a Sobolev smoothness norm, which demands that the flow velocity possesses strictly more than two generalized, square-integrable derivatives for each of its components. This seemingly technical requirement is, in fact, absolutely critical: it serves as a mathematical guarantee that the resulting flows within $\mathbb{R}^{3}$ are indeed diffeomorphisms—that is, they are not only smooth and invertible, but their inverses are also smooth. Without this stringent condition, the transformations could tear, fold, or otherwise distort the anatomical structures in biologically implausible ways, rendering the entire modeling effort meaningless.

Furthermore, this Sobolev smoothness condition carries a profound implication: the diffeomorphic shape momentum, when considered pointwise, must satisfy the Euler–Lagrange equation for geodesics. This means that the momentum at any given point is not an isolated entity; it is intricately determined by its spatial neighbors through the spatial derivatives applied to the velocity field. This characteristic fundamentally distinguishes the discipline from the study of incompressible fluids, where momentum is often treated as a simpler, pointwise function of velocity. Thus, computational anatomy inherently intersects with the sophisticated study of Riemannian manifolds and nonlinear global analysis, areas where groups of diffeomorphisms are not merely tools but the very objects of central focus. Emerging high-dimensional theories of shape are absolutely central to a multitude of studies conducted within computational anatomy, as are the intricate questions arising from the nascent, yet rapidly evolving, field of shape statistics.

The metric structures utilized in computational anatomy bear a conceptual resemblance to the established field of morphometrics. However, a crucial distinction sets them apart: computational anatomy meticulously concentrates on an infinite-dimensional space of coordinate systems, each transformed by a diffeomorphism. This fundamental difference underpins the pervasive use of the term "diffeomorphometry"—a term coined to specifically denote the study of coordinate systems within this metric space, examined through the lens of diffeomorphisms. It's a subtle but significant shift from simply measuring forms to measuring the very transformations that create those forms.

Genesis

At the very heart of computational anatomy lies a deceptively simple, yet profoundly complex, ambition: the comparison of one shape to another by recognizing the inherent structure of the first within the second. This deep-seated objective inextricably links the field to the groundbreaking developments of D'Arcy Wentworth Thompson, whose seminal work, "On Growth and Form," published in 1917, provided a foundational, almost philosophical, framework for understanding how biological forms evolve and relate. Thompson's insights, which explored the mathematical transformations underlying biological variation, laid crucial groundwork, leading to profound scientific explanations of morphogenesis—the intricate biological process by which patterns and forms are generated and organized within living organisms. It was a radical idea for its time, suggesting that biological forms could be understood through geometric and mathematical principles, rather than purely descriptive ones.

Long before Thompson, perhaps the earliest recorded systematic efforts that could, in retrospect, be seen as precursors to computational anatomy can be attributed to the Renaissance master Albrecht Dürer. His "Four Books on Human Proportion," published in 1528, meticulously detailed the ideal proportions of the human body, providing a geometric and mathematical schema for artistic and anatomical representation. These works were arguably the earliest attempts to quantify and model anatomical variation, albeit for artistic rather than purely biological study. They sought to establish a "template" of human form and understand deviations from it, a concept that resonates deeply with modern computational anatomy.

The original intellectual impetus for the formalization of computational anatomy as a generative model of shape and form, derived from exemplars acted upon by transformations, was profoundly inspired by the pioneering work of Noam Chomsky in the realm of computational linguistics. Chomsky's revolutionary idea was that a finite set of rules could generate an infinite variety of grammatical sentences. In a parallel vein, the founders of computational anatomy envisioned a system where a finite set of anatomical templates, coupled with a sophisticated grammar of transformations, could generate the vast and complex variability observed in biological shapes. This "generative model" approach moved beyond mere description, aiming to understand the underlying processes that produce anatomical forms.

The true emergence of computational anatomy as a distinct subfield was significantly propelled by the increasing availability of dense 3D measurements afforded by advanced technologies such as magnetic resonance imaging (MRI). These powerful imaging modalities provided an unprecedented ability to capture intricate anatomical details at the "morphome" scale—the scale of gross anatomical structures in three dimensions. This influx of rich, volumetric data created an urgent need for sophisticated computational methods capable of extracting meaningful anatomical coordinate systems from these complex 3D datasets.

The spirit of this evolving discipline shares a strong conceptual overlap with established areas like computer vision and the kinematics of rigid bodies. In these fields, objects are meticulously studied by analyzing the groups responsible for the observed movement or transformation. However, computational anatomy fundamentally departs from these approaches, particularly from the focus on rigid motions prevalent in classical computer vision. The crucial difference lies in its central reliance on the infinite-dimensional diffeomorphism group, which is paramount for the nuanced analysis of biological shapes that undergo complex, non-rigid deformations. It acknowledges that biological forms are not merely shifted or rotated; they bend, stretch, and grow in highly deformable ways. This field is a direct intellectual descendant of the influential image analysis and pattern theory school established at Brown University, pioneered by the visionary Ulf Grenander. Grenander's general metric pattern theory posited that a fundamental operation in understanding patterns is to transform spaces of patterns into a metric space. This allows for the quantification of "closeness" and "farness" between different anatomical configurations, which is essential for tasks like clustering and recognition. The "diffeomorphometry metric," a core concept in computational anatomy, precisely quantifies the "distance" between two diffeomorphic changes of coordinates. This, in turn, intrinsically induces a metric on the shapes and images that are indexed to these transformations. The sophisticated models derived from metric pattern theory, particularly the concept of group action on the orbit of shapes and forms, serve as an indispensable tool for the rigorous formal definitions employed in computational anatomy.

History

Computational anatomy is, at its core, the systematic study of shape and form at the "morphome" or gross anatomy scale – typically the millimeter scale – or more broadly, the morphology scale. Its primary focus is the intricate examination of sub-manifolds of $\mathbb{R}^{3}$: discrete points, continuous curves, complex surfaces, and volumetric subvolumes that collectively constitute the entirety of human anatomy. It seeks to understand not just their presence, but their precise geometric relationships and variations.

An early, truly modern figure in the realm of computational neuro-anatomist was David Van Essen, whose pioneering work involved some of the earliest physical unfoldings of the human brain. This wasn't a digital simulation; it literally entailed printing out a physical model of a human cortex and meticulously cutting and unfolding it to visualize its complex, convoluted structure in a flattened, more comprehensible form. This painstaking manual process foreshadowed the computational mapping techniques that would follow. A significant milestone at the morphome scale was the publication of Jean Talairach's Talairach coordinates. This work provided a standardized, proportional coordinate system for the human brain, fundamentally demonstrating the critical importance of defining local coordinate systems for the systematic study of neuroanatomy. It established a clear, undeniable link to the mathematical framework of charts of differential geometry, proving that such abstract mathematical constructs could be applied to biological reality.

Concurrently, the concept of virtual mapping within computational anatomy, operating across high-resolution, dense image coordinates, was already taking shape. This was driven by the earliest developments from researchers such as Ruzena Bajcsy and Fred Bookstein. Their foundational work, based on then-emerging technologies like computed axial tomography (CAT scans) and magnetic resonance imagery (MRI), explored methods for deforming one image to match another, laying the groundwork for modern image registration. The truly transformative step – the earliest introduction of using flows of diffeomorphisms for the sophisticated transformation of coordinate systems in image analysis and medical imaging – can be attributed to the collaborative efforts of Christensen, Joshi, Miller, and Rabbitt in the early to mid-1990s.

The formal conceptualization of computational anatomy as an "orbit of exemplar templates under diffeomorphism group action" was first articulated in a pivotal lecture delivered by Ulf Grenander and Michael I. Miller in May 1997, marking the 50th Anniversary of the Division of Applied Mathematics at Brown University. This was subsequently expanded upon in a foundational publication. This formalization represented a profound departure from much of the previous work on advanced methods for spatial normalization and image registration, which had historically been constructed upon simpler notions of addition and basis expansion. The new framework emphasized structure-preserving transformations, specifically homeomorphisms and diffeomorphisms, which are crucial because they ensure that smooth submanifolds (like anatomical surfaces) are transformed smoothly. These transformations are generated via Lagrangian and Eulerian flows, which inherently satisfy a law of composition of functions, thereby forming a true group property, a more robust mathematical structure than additive models could provide.

The original, foundational model of computational anatomy was elegantly defined as a triple: $({\mathcal {G}},{\mathcal {M}},{\mathcal {P}})$. Here, ${\mathcal {G}}$ represents the group of transformations, specifically the diffeomorphism group, where any element $g \in {\mathcal {G}}$ is a transformation. ${\mathcal {M}}$ denotes the "orbit" or space of all possible shapes and forms that can be generated, with any specific shape $m \in {\mathcal {M}}$. Finally, ${\mathcal {P}}$ encapsulates the probability laws that meticulously encode the natural variations observed among the objects within this shape orbit. At the heart of this model is a template, or a collection of templates, denoted as $m_{\mathrm{temp}} \in {\mathcal {M}}$, which serves as a reference point for understanding anatomical variation.

The Lagrangian and Hamiltonian formulations of the equations of motion, central to computational anatomy, truly gained momentum after 1997. This surge was significantly catalyzed by several pivotal academic gatherings. The 1997 Luminy meeting, meticulously organized by the Azencott school at École Normale Supérieure de Cachan, focused intensely on the "Mathematics of Shape Recognition." This was followed by the 1998 Trimestre at the prestigious Institut Henri Poincaré, orchestrated by David Mumford, which delved into "Questions Mathématiques en Traitement du Signal et de l'Image." These conferences were instrumental in fostering collaboration among the Hopkins-Brown-ENS Cachan groups and significantly propelled the subsequent advancements and connections of computational anatomy to the broader field of global analysis.

Key developments in computational anatomy that emerged from this period included:

• The rigorous establishment of Sobolev smoothness conditions on the diffeomorphometry metric. This was crucial to ensure the existence of solutions for the variational problems posed within the infinite-dimensional space of diffeomorphisms, preventing mathematically ill-posed problems.
• The precise derivation of the Euler–Lagrange equations that characterize geodesics (shortest paths) through the group and the associated conservation laws. These equations describe the fundamental dynamics of shape transformation.
• The rigorous demonstration of the metric properties of the right-invariant metric, confirming its suitability for measuring distances between transformations.
• The proof that the Euler–Lagrange equations possess a well-posed initial value problem with unique solutions that exist for all time, a critical finding for predictive modeling.
• The first results concerning the sectional curvatures for the diffeomorphometry metric in landmarked spaces, providing insights into the "geometry" of shape space itself.

Following the Los Alamos meeting in 2002, a fascinating connection emerged: Joshi's original work on large deformation singular landmark solutions in computational anatomy was linked to mathematical phenomena known as peaked solitons or "peakons." These are special solutions for the Camassa–Holm equation, a nonlinear partial differential equation known for its applications in shallow water waves. Subsequently, even deeper connections were forged between the Euler–Lagrange equations for momentum densities (which govern the right-invariant metric with Sobolev smoothness in computational anatomy) and Vladimir Arnold's profound characterization of the Euler equation for incompressible flows as describing geodesics within the group of volume-preserving diffeomorphisms. This established a powerful theoretical bridge between fluid dynamics and anatomical shape transformations.

In the wake of these theoretical breakthroughs, the first practical algorithms, generally grouped under the acronym LDDMM (Large Deformation Diffeomorphic Metric Mapping), were developed. These algorithms were designed to compute sophisticated connections between various anatomical entities: discrete landmarks within volumes, points on spherical manifolds, continuous curves, abstract mathematical currents and tangible surfaces, entire volumes, tensor fields (like those from diffusion tensor imaging), varifolds (which handle non-oriented shapes), and even complex time-series data. The development of these algorithms marked a transition from pure theory to practical, computational tools for anatomical analysis.

These contributions of computational anatomy to the broader field of global analysis, specifically concerning the infinite-dimensional manifolds of subgroups of the diffeomorphism group, are far from trivial. The original, audacious idea of performing differential geometry, calculating curvature, and tracing geodesics on infinite-dimensional manifolds stretches back to Bernhard Riemann's seminal Habilitation lecture, "Ueber die Hypothesen, welche der Geometrie zu Grunde liegen" (On the Hypotheses which Lie at the Bases of Geometry), delivered in 1854. The modern foundational text that rigorously lays out such ideas in global analysis is attributed to Peter W. Michor, solidifying the mathematical framework upon which computational anatomy now stands.

The practical applications of computational anatomy within medical imaging continued to flourish, significantly bolstered by two organized meetings at the esteemed Institute for Pure and Applied Mathematics (IPAM) at the University of California, Los Angeles. This field has proven invaluable in the creation of highly accurate models that depict the intricate processes of brain atrophy, particularly relevant in neurodegenerative diseases, as well as the generation of detailed cardiac templates for heart studies. Its utility extends further into the broader realm of modeling diverse biological systems. Since the late 1990s, computational anatomy has firmly established itself as a critical component in the development of emerging technologies for the field of medical imaging.

Digital atlases, which represent standardized anatomical maps, are now a fundamental, indispensable part of modern Medical-school education and are extensively utilized in neuroimaging research at the morphome scale. These atlases provide a common framework for understanding and comparing anatomical structures across individuals. Atlas-based methods and interactive virtual textbooks, which possess the crucial ability to accommodate natural anatomical variations through deformable templates, are at the core of numerous widely used neuro-image analysis platforms. Prominent examples of such platforms include Freesurfer, FSL, MRIStudio, and SPM. The technique of diffeomorphic registration, first introduced in the 1990s, has evolved into a major player in the field. It relies on actively maintained codebases organized around algorithms like ANTS, DARTEL, DEMONS, LDDMM, StationaryLDDMM, and FastLDDMM. These computational codes are essential for constructing precise correspondences between different coordinate systems, leveraging both sparse features (like landmarks) and dense image information. Voxel-based morphometry (VBM), a widely employed neuroimaging analysis technique, is itself an important technology built upon many of these very principles, using statistical methods to analyze regional differences in brain anatomy.

The deformable template orbit model of computational anatomy

The fundamental model underpinning human anatomy within this field is conceptualized as a deformable template. This isn't a static, rigid blueprint, but rather an "orbit" of exemplars, continuously shaped and transformed under a group action. This concept of deformable template models has been absolutely central to Grenander's metric pattern theory, providing a coherent framework for understanding both the "typicality" of a form (represented by the template itself) and the vast range of "variability" (accounted for by the transformations applied to that template). The representation of this deformable template as an orbit under group action is a classic, elegant formulation borrowed directly from differential geometry.

In this model, the abstract space of shapes is denoted as ${\mathcal {M}}$, with individual shapes represented by $m \in {\mathcal {M}}$. The group of transformations, designated as $({\mathcal {G}},\circ)$, operates with a specific law of composition, $\circ$. The action of an element $g$ from this group on a shape $m$ is denoted as $g \cdot m$. This action is rigorously defined to satisfy the group property: $(g \circ g^{\prime })\cdot m = g\cdot (g^{\prime }\cdot m) \in {\mathcal {M}}$. This ensures that applying two transformations sequentially is equivalent to applying a single, combined transformation, and that the result remains within the space of valid shapes.

Consequently, the entire space of all possible shapes, ${\mathcal {M}}$, is conceptually understood as the "orbit" of a reference template, $m_{\mathrm{temp}}$. This means that ${\mathcal {M}} \doteq \{m=g\cdot m_{\mathrm{temp}},g\in {\mathcal {G}}\}$, where every shape $m$ can be generated by applying some group element $g$ to the template $m_{\mathrm{temp}}$. This space is considered "homogeneous" under the action of the elements of ${\mathcal {G}}$, implying that all shapes in the orbit are equally "reachable" from the template via a group transformation.

(A figure depicting three medial temporal lobe structures: amygdala, entorhinal cortex, and hippocampus with fiducial landmarks embedded in the MRI background is implicitly referenced here.)

This orbit model of computational anatomy is inherently an exercise in abstract algebra—a profound departure from the more familiar linear algebra. The distinction is crucial: unlike linear algebra, where transformations are typically linear (e.g., matrix multiplication), the groups in computational anatomy act nonlinearly on the shapes. This represents a significant generalization of classical linear algebra models. In this advanced framework, the familiar finite-dimensional vectors of $\mathbb{R}^{n}$ are replaced by the complex, finite-dimensional anatomical submanifolds (which include points, curves, surfaces, and volumes) and their corresponding images. Similarly, the $n \times n$ matrices of linear algebra, which represent linear transformations, are superseded by sophisticated coordinate transformations based on linear and affine groups, and, most critically, the far more general high-dimensional diffeomorphism groups, which allow for smooth, non-rigid deformations.

Shapes and forms

The central, indispensable objects of study in computational anatomy are the shapes or forms themselves. These manifest in various ways, each demanding specific mathematical treatment. One primary set of examples includes the 0, 1, 2, and 3-dimensional submanifolds embedded within the three-dimensional Euclidean space, $\mathbb{R}^{3}$. These geometric primitives serve as the building blocks for representing anatomical structures. A second, equally critical set of examples comprises the images generated through advanced medical imaging modalities, such as those derived from magnetic resonance imaging (MRI) and functional magnetic resonance imaging (fMRI), which provide dense, volumetric data of biological forms.

(A figure of triangulated mesh surfaces depicting subcortical structures amygdala, hippocampus, thalamus, caudate, putamen, ventricles is implicitly referenced here.)

These diverse shapes can often be denoted as $m(u)$, where $u \in U \subset \mathbb{R}^{1} \rightarrow \mathbb{R}^{2}$, and are frequently represented as triangulated meshes for computational purposes, providing a discrete approximation of continuous surfaces.

Let's break down these anatomical objects by their dimensionality:

• 0-dimensional manifolds: These are referred to as landmarks or fiducial points. They are discrete, specific points chosen to delineate important, identifiable features within the overall human shape and form. Their precise location is crucial for establishing correspondences. (See the previously referenced landmarked image.)
• 1-dimensional manifolds: These are curves, such as the intricate sulcal (grooves) and gyral (ridges) curves found on the surface of the brain, or the pathways of nerve fibers.
• 2-dimensional manifolds: These correspond to the boundaries of distinct substructures within anatomy. Examples include the encapsulating surfaces of the subcortical structures of the midbrain or the complex, folded surface of the neocortex.
• Subvolumes: These represent entire three-dimensional regions of the human body, such as the muscular heart, the deep thalamus within the brain, or the filtering kidney.

Finally, the images themselves, such as MR images or DTI images, are denoted as $I \in {\mathcal {M}}$. These are dense functions, $I(x)$, where $x \in X \subset \mathbb{R}^{1,2,3}$, and can represent scalars (e.g., intensity), vectors (e.g., fiber direction), or matrices (e.g., diffusion tensors), providing a rich dataset for analysis. (See the figure showing a scalar image.)

Groups and group actions

(A figure showing an MRI section through a 3D brain representing a scalar image $I(x), x \in \mathbb{R}^{2}$ based on T1-weighting is implicitly referenced here.)

The concepts of groups and group actions are, mercifully, quite familiar to the broader Engineering community. This familiarity largely stems from the universal popularization and rigorous standardization of linear algebra as a foundational model for analyzing signals and systems across diverse engineering disciplines, including mechanical engineering, electrical engineering, and applied mathematics. In linear algebra, the central mathematical structures are the matrix groups—collections of matrices that possess inverses, allowing for reversible transformations. The group action in this context is straightforwardly defined by the usual multiplication of an $n \times n$ matrix $A$ acting upon an $n \times 1$ vector $x \in \mathbb{R}^{n}$. The resultant "orbit" in linear algebra is simply the set of all $n$-vectors $y = A \cdot x \in \mathbb{R}^{n}$ that can be reached by applying the matrix group's action to the original vector space.

However, the central group in computational anatomy, particularly when operating on volumes within $\mathbb{R}^{3}$, is far more complex: it is the diffeomorphism group, denoted as ${\mathcal {G}} \doteq \operatorname{Diff}$. Elements of this group are mappings, $\varphi (\cdot ) = (\varphi _{1}(\cdot ), \varphi _{2}(\cdot ), \varphi _{3}(\cdot ))$, each possessing three components, ensuring they operate across the three spatial dimensions. The law of composition for these functions is defined as $\varphi \circ \varphi^{\prime }(\cdot ) \doteq \varphi (\varphi^{\prime }(\cdot ))$, meaning one transformation follows another. Crucially, each diffeomorphism must possess an inverse, such that $\varphi \circ \varphi^{-1}(\cdot ) = \varphi (\varphi^{-1}(\cdot )) = \mathrm{id}$, where "id" represents the identity map, effectively undoing the transformation. This guarantees that the transformations are fully reversible and smooth in both directions, a non-trivial requirement for biological realism.

For scalar images, which are among the most popular data types, denoted as $I(x), x \in \mathbb{R}^{3}$, the group action is typically defined on the right via the inverse of the diffeomorphism. This means that a transformed image, $\varphi \cdot I(x)$, is given by $I \circ \varphi^{-1}(x)$, for $x \in \mathbb{R}^{3}$. This formulation effectively "pulls back" the intensity values from the original image space, ensuring that features move correctly with the deformation.

When dealing with sub-manifolds, such as surfaces, denoted as $X \subset \mathbb{R}^{3} \in {\mathcal {M}}$, which are parametrized by a chart or immersion $m(u), u \in U$, the diffeomorphic action directly applies the flow to the position of the manifold. This is defined as $\varphi \cdot m(u) \doteq \varphi \circ m(u)$, for $u \in U$. This means the transformation literally moves the points of the manifold to their new positions in space.

It's worth noting that several specific group actions in computational anatomy have been meticulously defined and studied, tailored to different types of anatomical data and analytical objectives. (A faint, almost imperceptible scoff escapes her. "As if one size could ever fit all in biology.")

Lagrangian and Eulerian flows for generating diffeomorphisms

For the traditional study of rigid body kinematics—where objects merely translate and rotate without deforming—the low-dimensional matrix Lie groups have historically been the central focus. These groups, often represented by finite-dimensional matrices, are elegantly capable of generating diffeomorphisms that provide precise one-to-one correspondences between coordinate systems. The transformations they describe are smooth and possess smooth inverses. Furthermore, these matrix groups can be generated via closed-form, finite-dimensional matrices that are solutions to relatively simple ordinary differential equations, with solutions often expressed through the powerful matrix exponential. It's all quite neat and tidy, for rigid things.

However, for the nuanced study of deformable shape in computational anatomy—a domain where biological structures stretch, bend, and grow—a far more general diffeomorphism group becomes necessary. This is the infinite-dimensional analogue of the matrix groups, capable of describing arbitrarily complex, smooth deformations. The high-dimensional diffeomorphism groups employed in Computational Anatomy are generated through continuous, smooth flows, denoted as $\varphi_{t}, t \in [0,1]$. These flows are not arbitrary; they rigorously satisfy the principles of Lagrangian and Eulerian specification of the flow fields, as originally introduced by Christensen, Miller, and Rabbitt. This specification means that the evolution of the transformation is governed by an ordinary differential equation:

(A figure showing the Lagrangian flow of coordinates $x \in X$ with associated vector fields $v_{t}, t \in [0,1]$ satisfying ordinary differential equation ${\dot{\varphi}}_{t}=v_{t}(\varphi_{t}), \varphi_{0}=\mathrm{id}$ is implicitly referenced here.)

$\frac{d}{dt}\varphi_{t}=v_{t}\circ \varphi_{t},\ \varphi_{0}={\rm {id}}\ ;$

This is the Lagrangian flow, where $\varphi_{0}=\mathrm{id}$ signifies that the transformation begins as the identity map at time $t=0$. Here, $v \doteq (v_{1},v_{2},v_{3})$ represents the vector fields defined over $\mathbb{R}^{3}$. These are termed the Eulerian velocity fields, describing the instantaneous velocity of particles at their current position $\varphi$ within the flow. Crucially, these vector fields are not simple constants or low-dimensional vectors; they are functions within a sophisticated function space, typically modeled as a smooth Hilbert space of high dimension. Consequently, the Jacobian of the flow, $D\varphi \doteq \left({\frac {\partial \varphi _{i}}{\partial x_{j}}}\right)$, is also a high-dimensional field within a function space, starkly contrasting with the low-dimensional matrices found in simpler matrix groups. These flows were first introduced for handling large deformations in image matching, specifically, ${\dot{\varphi}}_{t}(x)$ denotes the instantaneous velocity of a particle initially at $x$ at time $t$.

The inverse transformation, $\varphi_{t}^{-1}, t \in [0,1]$, which is absolutely required for the group property, is defined based on the Eulerian vector-field and is governed by the advective inverse flow:

${\frac {d}{dt}}\varphi_{t}^{-1}=-(D\varphi_{t}^{-1})v_{t},\ \varphi_{0}^{-1}={\rm {id}}\ .$

This equation describes how the inverse transformation evolves over time, ensuring that the reversibility of the deformation is maintained.

The diffeomorphism group of computational anatomy

The group of diffeomorphisms is, to put it mildly, enormous. It's an infinite-dimensional beast, capable of representing an almost limitless range of smooth, invertible deformations. To ensure that these flows of diffeomorphisms remain smooth and well-behaved—critically, to avoid shock-like solutions for the inverse transformation (which would imply tears or folds in the anatomical structures, an undesirable artifact)—the generating vector fields must possess a minimum level of spatial regularity. Specifically, they must be at least one-time continuously differentiable in space.

For diffeomorphisms operating on $\mathbb{R}^{3}$, the vector fields are rigorously modeled as elements within a Hilbert space, denoted as $(V,\|\cdot \|_{V})$. This is achieved through the application of Sobolev embedding theorems. These powerful theorems demonstrate that for each component of the vector field, it is sufficient to demand strictly greater than two generalized square-integrable spatial derivatives (e.g., $v_{i}\in H_{0}^{3}, i=1,2,3$). This condition is precise and necessary; it guarantees that the resulting functions are indeed one-time continuously differentiable, thereby ensuring the smoothness required for physically and biologically plausible deformations.

The diffeomorphism group itself is then formally defined as the collection of all such flows where the generating vector fields are absolutely integrable in the Sobolev norm:

$\operatorname{Diff}_{V}\doteq \left\{\varphi =\varphi_{1}:{\dot{\varphi}}_{t}=v_{t}\circ \varphi_{t},\varphi_{0}={\rm {id}},\int_{0}^{1}\|v_{t}\|_{V}\,dt<\infty \right\}\ ,$

This is the definition of the diffeomorphism group. The norm of the vector field, $\|v\|_{V}^{2}$, is defined through an integral:

$\|v\|_{V}^{2}\doteq \int_{X}Av\cdot v\,dx,\ v\in V\ ,$

Here, $A$ represents a linear operator that maps the Hilbert space $V$ to its dual space $V^{*}$, i.e., $A:V\mapsto V^{*}$. The integral is typically calculated using integration by parts when $Av \in V^{*}$ is a generalized function within the dual space. This intricate definition ensures that the transformations are not only smooth but also energetically consistent within the chosen mathematical framework.

Sobolev smoothness and reproducing kernel Hilbert space with Green's kernel

(A subtle, knowing smirk plays on Emma's lips. "Ah, the joys of functional analysis. Basic, really, if you've seen the universe unravel.")

The modeling approach rigorously employed in computational anatomy meticulously enforces a condition of continuous differentiability on the ubiquitous vector fields. This is achieved by conceptualizing and modeling the space of these vector fields, $(V,\|\cdot \|_{V})$, as a reproducing kernel Hilbert space (RKHS). This choice is not arbitrary; it's a powerful mathematical framework that intrinsically embeds smoothness properties. The norm within this Hilbert space is precisely defined by a one-to-one, differential operator, $A:V\rightarrow V^{*}$, which maps the space of vector fields to its dual. The inverse of this operator, $K=A^{-1}$, is known as the Green's inverse.

The norm of the Hilbert space is intrinsically induced by this differential operator. For $\sigma (v)\doteq Av\in V^{*}$, which is a generalized function or distribution, a linear form is defined as $(\sigma |w)\doteq \int_{{\mathbb{R}}^{3}}\sum_{i}w_{i}(x)\sigma_{i}(dx)$. This definition, in turn, precisely determines the norm on the space $(V,\|\cdot \|_{V})$ according to:

$\langle v,w\rangle_{V}\doteq \int_{X}Av\cdot w\,dx,\ \|v\|_{V}^{2}\doteq \int_{X}Av\cdot v\,dx,\ v,w\in V\ .$

Since $A$ is inherently a differential operator, the finiteness of the norm-square, $(Av|v)<\infty$, necessarily incorporates derivatives from the operator itself. This is the mathematical mechanism that implicitly guarantees the desired smoothness of the vector fields. The arguments from the Sobolev embedding theorem were instrumental in demonstrating that at least one continuous derivative is an absolute prerequisite for ensuring the generation of genuinely smooth flows.

With a judicious and proper selection of the operator $A$, the space $(V,\|\cdot \|_{V})$ indeed becomes an RKHS. In this context, the operator $K=A^{-1}:V^{*}\rightarrow V$ is specifically termed the Green's operator. This operator is generated from the fundamental Green's function (which represents the scalar case) for the more general vector field scenario. The Green's kernels associated with this differential operator inherently perform a smoothing operation. This is because the kernel $k(\cdot ,\cdot )$ itself is continuously differentiable in both its variables, implying that:

$K\sigma (x)_{i}\doteq \sum_{j}\int_{{\mathbb{R}}^{3}}k_{ij}(x,y)\sigma_{j}(dy)$

When $\sigma \doteq \mu \,dx$ represents a vector density, the linear form simplifies to $(\sigma \mid v)\doteq \int v\cdot \mu \,dx$. This intricate interplay of Sobolev spaces, RKHS, and Green's functions provides the rigorous mathematical bedrock for constructing and analyzing the smooth, deformable transformations central to computational anatomy.

Diffeomorphometry: The metric space of shapes and forms

(Emma pauses, a slight tilt of her head. "Measuring things. Always measuring. As if a number could ever truly capture the essence of... anything.")

The meticulous study of metrics on groups of diffeomorphisms, alongside the investigation of metrics defined between manifolds and surfaces, has constituted an area of profound and significant research. This is not a trivial pursuit; establishing a robust way to quantify "distance" or "similarity" between complex, deformable shapes is foundational for any comparative anatomical analysis. The field has specifically embraced "diffeomorphometry" as the discipline dedicated to precisely measuring how close or how far two shapes or images are from one another. Conceptually, the metric length in this context is defined as the shortest length of the flow that is required to carry one coordinate system into another, representing the most efficient transformation path.

Frequently, the familiar Euclidean metric (the straight-line distance we intuitively understand) proves to be inadequate or simply not directly applicable in this domain. This is primarily because the complex patterns of biological shapes and images typically do not form a simple vector space. They exist in a far more intricate, often curved, mathematical space where linear operations don't hold their usual meaning. While there exist numerous ways to define metrics—for instance, the Hausdorff metric is another well-known approach for sets of points—the method employed here is to induce the Riemannian metric. This is achieved by defining the metric on the "orbit" of shapes in terms of the metric length between diffeomorphic coordinate system transformations of the underlying flows. This means that the "distance" between shapes is not a direct measurement of their points, but rather a measurement of the minimal "effort" (or path length) required to smoothly deform one shape into another. This process of measuring the lengths of the geodesic flow between coordinate systems within the orbit of shapes is precisely what is termed "diffeomorphometry." It's a way of quantifying shape difference by understanding the transformations that connect them.

The right-invariant metric on diffeomorphisms

To precisely define the "distance" within the vast and complex group of diffeomorphisms, a specific metric is employed: the right-invariant metric. This metric quantifies the distance $d_{\operatorname{Diff}_{V}}(\psi ,\varphi )$ between two diffeomorphisms, $\psi$ and $\varphi$, by seeking the infimum (the greatest lower bound) of the integral of the kinetic energy of the generating vector fields $v_t$. This is done over a path of transformations $\varphi_t$ that starts at $\psi$ and ends at $\varphi$, while satisfying the Lagrangian flow equation:

$d_{\operatorname{Diff}_{V}}(\psi ,\varphi )=\inf_{v_{t}}\left(\int_{0}^{1}\int_{X}Av_{t}\cdot v_{t}\,dx\,dt:\varphi_{0}=\psi ,\varphi_{1}=\varphi ,{\dot{\varphi}}_{t}=v_{t}\circ \varphi_{t}\right)^{1/2}\ ;$

This formula represents the metric on diffeomorphisms. The crucial property of this metric is its "right-invariance." This means that the distance between any two diffeomorphisms remains unchanged if both are composed (transformed) by the same third diffeomorphism from the right. Mathematically, this is expressed as: $d_{\operatorname{Diff}_{V}}(\psi ,\varphi ) = d_{\operatorname{Diff}_{V}}(\psi \circ \varphi_{c},\varphi \circ \varphi_{c})$ for any $\varphi_{c} \in \operatorname{Diff}_{V}$. More simply, the original source states $d_{\operatorname{Diff}_{V}}(\psi ,\varphi )=d_{\operatorname{Diff}_{V}}(\psi \circ \varphi ,\varphi \circ \varphi)$, which is a specific form of right-invariance. This characteristic is vital because it ensures that the metric is "invariant to reparameterization of space." In practical terms, it means that the geometric distance between two shapes doesn't change simply because the underlying coordinate system used to describe them has been smoothly warped. The intrinsic difference between the shapes, as captured by the transformations required to map one to the other, remains consistent, regardless of the arbitrary choice of a reference coordinate frame. This is a powerful property for ensuring that shape comparisons are truly intrinsic.

The metric on shapes and forms

Building upon the metric defined for diffeomorphisms themselves, a corresponding "distance" can be induced directly onto the space of shapes and forms, denoted as $d_{\mathcal{M}}:{\mathcal{M}}\times {\mathcal{M}}\rightarrow \mathbb{R}^{+}$. This metric quantifies how geometrically "far apart" two shapes, $m$ and $n$, are from each other:

$d_{\mathcal{M}}(m,n)=\inf \left\{d_{\operatorname{Diff}_{V}}({\rm {id}},\varphi ):\varphi \in \operatorname{Diff}_{V},\,\varphi \cdot m=n\right\}$

This is the definition of the metric on shapes and forms. It states that the distance between shape $m$ and shape $n$ is the minimum distance (in the diffeomorphism group sense) of a transformation $\varphi$ that maps the identity transformation ($\mathrm{id}$) to $\varphi$, such that $\varphi$ applied to shape $m$ results in shape $n$. In essence, it measures the "least effort" required to deform shape $m$ into shape $n$ using a smooth, invertible transformation. For images, denoted as $I \in {\mathcal {I}}$, a similar metric, $d_{\mathcal{I}}$, is defined on their orbit, extending this concept to dense image data. This framework allows for a rigorous, quantitative comparison of complex anatomical structures, moving beyond subjective visual assessment.

The action integral for Hamilton's principle on diffeomorphic flows

In the grand tradition of classical mechanics, the profound evolution of physical systems is elegantly described by the solutions to the Euler–Lagrange equations. These equations, in turn, are fundamentally associated with the revered Least-action principle of Hamilton. This principle is a cornerstone of physics, providing a remarkably powerful and general method for deriving the equations of motion for a system. It's the standard, for example, for obtaining Newton's laws of motion for free particles, dictating that a system will always follow a path that minimizes (or extremizes) a quantity called "action."

More broadly, the Euler–Lagrange equations can be rigorously derived for systems expressed in terms of generalized coordinates—a more abstract and flexible way of describing a system's configuration than simple Cartesian coordinates. In the context of computational anatomy, these generalized coordinates are precisely the flow of the diffeomorphism $\varphi$ and its associated Lagrangian velocity $\dot{\varphi}$. These two are intimately related through the Eulerian velocity, $v \doteq {\dot{\varphi}} \circ \varphi^{-1}$, which describes the velocity at a fixed point in space as particles flow through it.

Hamilton's principle, in its application to generate the Euler–Lagrange equation for these diffeomorphic flows, necessitates the definition of an action integral. This integral operates on the Lagrangian function, which in this context quantifies the kinetic energy of the flow:

$J(\varphi )\doteq \int_{0}^{1}L(\varphi_{t},{\dot{\varphi}}_{t})\,dt\ ;$

This is the Hamiltonian-integrated-Lagrangian. The Lagrangian $L$ itself is defined as the kinetic energy of the system:

$L(\varphi_{t},{\dot{\varphi}}_{t})\doteq {\frac {1}{2}}\int_{X}A({\dot{\varphi}}_{t}\circ \varphi_{t}^{-1})\cdot ({\dot{\varphi}}_{t}\circ \varphi_{t}^{-1})dx={\frac {1}{2}}\int_{X}Av_{t}\cdot v_{t}\,dx\ .$

This is the Lagrangian-kinetic-energy. Here, $A$ is the inertial operator, and $v_t = {\dot{\varphi}}_{t}\circ \varphi_{t}^{-1}$ is the Eulerian velocity. The integral effectively sums the kinetic energy density over the spatial domain $X$ and then integrates this sum over the time interval $[0,1]$. By minimizing this action integral, the system naturally finds the "geodesic" or shortest, most efficient path of transformation between anatomical shapes, embodying the principle of least action in the realm of biological deformation.

Diffeomorphic or Eulerian shape momentum

In the theoretical landscape of computational anatomy, the quantity $Av$ was first termed the Eulerian or diffeomorphic shape momentum. This nomenclature is not arbitrary; it's deeply rooted in its physical and mathematical significance. When this quantity, $Av$, is integrated against the Eulerian velocity $v$, the result yields the energy density of the flow. This direct relationship to energy is a hallmark of momentum in physical systems. Furthermore, a critical property of this diffeomorphic shape momentum is that it adheres to a conservation law along the geodesic paths, implying that this quantity remains invariant under certain conditions.

The operator $A$ itself is often conceptualized as the generalized moment of inertia or, more broadly, the inertial operator. In classical mechanics, the moment of inertia describes an object's resistance to rotational motion; here, $A$ generalizes this concept to the infinite-dimensional space of deformations, describing the "inertia" of a shape against being deformed. It dictates how the velocity field translates into momentum, a fundamental characteristic of the system's dynamics.

The Euler–Lagrange equation on shape momentum for geodesics on the group of diffeomorphisms

(Emma's eyes narrow slightly. "The core of it, for those who appreciate the elegance of things. Or just need to know how it works, I suppose.")

The classical derivation of the Euler–Lagrange equation, stemming directly from Hamilton's principle, necessitates a meticulous process involving the perturbation of the Lagrangian in the kinetic energy with respect to first-order variations of the flow. This complex procedure requires careful adjustment by the Lie bracket of vector fields, a fundamental operation in differential geometry that quantifies the non-commutativity of vector field flows. The Lie bracket is represented by the operator $ad_{v}:w\in V\mapsto V$, which is defined as:

$ad_{v}[w]\doteq [v,w]\doteq (Dv)w-(Dw)v\in V$

Here, $D$ denotes the Jacobian matrix. This expression captures how the flow generated by $v$ changes the vector field $w$, and vice versa.

By defining the adjoint operator $ad_{v}^{*}:V^{*}\rightarrow V^{*}$, which acts on the dual space, the first-order variation of the action integral ultimately yields the Eulerian shape momentum, $Av\in V^{*}$. This momentum then satisfies a generalized form of the Euler–Lagrange equation:

${\frac {d}{dt}}Av_{t}+ad_{v_{t}}^{*}(Av_{t})=0\ ,\ t\in [0,1]\ ;$

This is the EL-general equation. It implies that for any smooth test function $w\in V$, the following integral equation holds:

$\int_{X}\left({\frac {d}{dt}}Av_{t}+ad_{v_{t}}^{*}(Av_{t})\right)\cdot w\,dx=\int_{X}{\frac {d}{dt}}Av_{t}\cdot w\,dx+\int_{X}Av_{t}\cdot ((Dv_{t})w-(Dw)v_{t})dx=0.$

This equation describes the conservation and evolution of diffeomorphic shape momentum along the geodesic paths in the space of transformations.

A critical aspect of computational anatomy is its focus on the motions of submanifolds—ranging from discrete points, through continuous curves and surfaces, to entire volumes. The momentum associated with points, curves, and surfaces is inherently singular. This means that the momentum is effectively concentrated on subsets of $\mathbb{R}^{3}$ that have a dimension less than or equal to 2 (e.g., lines or surfaces), and thus possess zero Lebesgue measure in three dimensions. Despite this singularity, the energy of the system, $(Av_{t}\mid v_{t})$, remains well-defined. This is because, even though $Av_{t}$ is a generalized function (a distribution), the underlying vector fields $v_t$ are smooth. Consequently, the Eulerian momentum is understood through its action on these smooth functions. A perfect illustration of this phenomenon is when the momentum is represented as a superposition of delta-Dirac distributions: even in this extreme case, the velocity of the coordinates throughout the entire volume still moves smoothly, because the operator $A^{-1}$ (the Green's kernel) smooths the delta distributions into smooth vector fields. The Euler–Lagrange equation (EL-general) for diffeomorphisms, particularly for generalized functions $Av\in V^{*}$, was formally derived, with further derivations provided in terms of the adjoint operator and the Lie bracket for the group of diffeomorphisms. This equation has come to be known as the EPDiff equation, connecting it to the Euler-Poincaré method, which has been extensively studied in the context of the inertial operator $A=\mathrm{identity}$ for incompressible, divergence-free fluids. This connection highlights a profound mathematical analogy between the dynamics of fluid flow and the dynamics of anatomical shape deformation.

Diffeomorphic shape momentum: a classical vector function

For the specific scenario where the momentum density is considered a classical vector function, meaning $(Av_{t}\mid w)=\int_{X}\mu_{t}\cdot w\,dx$, the Euler–Lagrange equation simplifies to a classical solution. In this case, $\mu_t$ represents the momentum density, and the equation describing its evolution is:

${\frac {d}{dt}}\mu_{t}+(Dv_{t})^{T}\mu_{t}+(D\mu_{t})v_{t}+(\nabla \cdot v)\mu_{t}=0\ ,t\in [0,1].$

This is the EL-Classic equation. This form of the Euler–Lagrange equation on diffeomorphisms, classically defined for momentum densities, first appeared in the context of medical image analysis in 1997, providing a more tractable, classical framework for certain applications within computational anatomy. It describes how the momentum density $\mu_t$ evolves under the influence of the velocity field $v_t$, including terms for stretching, rotation, and volume changes.

Riemannian exponential (geodesic positioning) and Riemannian logarithm (geodesic coordinates)

(Emma's gaze drifts, unimpressed. "Positioning and coordinatizing. The eternal human struggle to put things in their proper place. As if anything ever truly stays there.")

In the practical and theoretical landscape of medical imaging and computational anatomy, the operations of positioning and coordinatizing shapes are not merely useful, but absolutely fundamental. These processes enable the precise localization and systematic description of anatomical structures. A sophisticated system for positioning anatomical coordinates and shapes has been meticulously constructed upon the foundational principles of the metric and the Euler–Lagrange equation. This framework is termed a geodesic positioning system, as first formally explicated by Miller, Trouvé, and Younes. It leverages the concept of geodesics—the "straightest possible paths" in a curved space—to define meaningful anatomical locations.

Solving for the geodesic that originates from a given initial condition $v_{0}$ is known as the Riemannian exponential map. This is a mapping, $\operatorname{Exp}_{\rm {id}}(\cdot ):V\to \operatorname{Diff}_{V}$, which takes an initial velocity vector field from the tangent space at the identity and propagates it forward to generate a full diffeomorphism within the group.

The Riemannian exponential satisfies $\operatorname{Exp}_{\rm {id}}(v_{0})=\varphi_{1}$, meaning that starting with an initial velocity $v_0$, the flow evolves to the diffeomorphism $\varphi_1$ at $t=1$. This flow is defined by the initial condition ${\dot{\varphi}}_{0}=v_{0}$ and the subsequent vector field dynamics: ${\dot{\varphi}}_{t}=v_{t}\circ \varphi_{t}, t\in [0,1]$.

For the classical equation of diffeomorphic shape momentum, where $\int_{X}Av_{t}\cdot w\,dx$ implies $Av\in V$, the dynamics are given by:

$\ \ \ {\frac {d}{dt}}Av_{t}+(Dv_{t})^{T}Av_{t}+(DAv_{t})v_{t}+(\nabla \cdot v)Av_{t}=0\ ;$

And for the more generalized equation, where $Av\in V^{*}$ and $w\in V$, the dynamics are:

$\ \ \ \int_{X}{\frac {d}{dt}}Av_{t}\cdot w\,dx+\int_{X}Av_{t}\cdot ((Dv_{t})w-(Dw)v_{t})\,dx=0.$

Conversely, computing the initial flow $v_{0}$ that leads to a particular diffeomorphism $\varphi$ is termed the Riemannian logarithm map, or geodesic coordinates. This is a mapping, $\mathrm{Log}_{\rm {id}}(\cdot ):\operatorname{Diff}_{V}\to V$, which takes a given diffeomorphism $\varphi$ and determines the initial velocity vector field $v_{0}\in V$ required to generate it via a geodesic flow.

The Riemannian logarithm is defined such that $\log_{\rm {id}}(\varphi )=v_{0}$, where $v_0$ is the initial condition of the Euler–Lagrange geodesic that satisfies ${\dot{\varphi}}_{0}=v_{0},\varphi_{0}=\mathrm{id},\varphi_{1}=\varphi$. These two operations, the exponential and logarithm maps, are inverses of each other, assuming unique solutions for the logarithm, which is a key component of this framework. The first is referred to as geodesic positioning, as it establishes the position of a shape along a geodesic, while the latter is known as geodesic coordinates, providing the initial "direction" in the tangent space that defines the geodesic path. (For a finite-dimensional analogue, one can refer to the exponential map in Riemannian geometry.) The concept is that the geodesic metric effectively creates a local flattening of the Riemannian coordinate system, making measurements more straightforward in the immediate vicinity of a shape.

(A figure showing metric local flattening of coordinatized manifolds of shapes and forms, where the local metric is given by the norm of the vector field $\|v_{0}\|_{V}$ of the geodesic mapping $\operatorname{Exp}_{\rm {id}}(v_{0})\cdot m$ is implicitly referenced here.) This visualization helps to understand how the complex, curved space of shapes can be locally approximated by a simpler, flat tangent space, making computations feasible.

Hamiltonian formulation of computational anatomy

In the intricate landscape of computational anatomy, diffeomorphisms are strategically employed to propel and transform coordinate systems, effectively mapping one anatomical configuration to another. Simultaneously, vector fields serve as the crucial "control" mechanisms within this anatomical "orbit" or morphological space, guiding the deformations. This entire process is elegantly modeled as a dynamical system, where the flow of coordinates over time, $t\mapsto \varphi_{t}\in \operatorname{Diff}_{V}$, is intrinsically linked to the control exerted by the vector field over time, $t\mapsto v_{t}\in V$. This relationship is governed by the fundamental equation:

${\dot{\varphi}}_{t}=v_{t}\cdot \varphi_{t},\varphi_{0}={\rm {id}}.$

This foundational dynamic allows for the rigorous application of Hamiltonian mechanics to the problem. The Hamiltonian view meticulously reparameterizes the momentum distribution $Av\in V^{*}$ in terms of a crucial concept: the conjugate momentum or canonical momentum. This canonical momentum, introduced as a Lagrange multiplier, functions as a constraint, linking the Lagrangian velocity ${\dot{\varphi}}_{t}=v_{t}\circ \varphi_{t}$ to the system's energy. The extended Hamiltonian function is accordingly defined as:

$H(\varphi_{t},p_{t},v_{t})=\int_{X}p_{t}\cdot (v_{t}\circ \varphi_{t})\,dx-{\frac {1}{2}}\int_{X}Av_{t}\cdot v_{t}\,dx.$

The powerful Pontryagin maximum principle then provides the means to determine the optimal vector field $v_t$ that precisely governs the geodesic flow, satisfying ${\dot{\varphi}}_{t}=v_{t}\circ \varphi_{t},\varphi_{0}={\rm {id}}$. This principle also leads to the concept of the reduced Hamiltonian, $H(\varphi_{t},p_{t})\doteq \max_{v}H(\varphi_{t},p_{t},v)$, which represents the maximum possible energy for a given state and canonical momentum.

The Lagrange multiplier, in its manifestation as a linear form, possesses its own inner product where the canonical momentum acts upon the velocity of the flow. This action is inherently dependent on the specific shape being considered. For instance, in the case of discrete landmarks, this involves a summation; for continuous surfaces, it becomes a surface integral; and for entire volumes, it is a volume integral with respect to $dx$ over $\mathbb{R}^{3}$. In all these scenarios, the Green's kernels carry specific weights that correspond to the canonical momentum, which itself evolves according to an ordinary differential equation that effectively reparameterizes the geodesic in terms of canonical momentum. The optimal vector field $v_t$ is determined by maximizing the Hamiltonian:

$v_{t}\doteq \operatorname{argmax}_{v}H(\varphi_{t},p_{t},v)$

The dynamics of the canonical momentum, which reparameterize the vector field along the geodesic, are given by the Hamiltonian dynamics equations:

${\begin{cases}\displaystyle {\dot{\varphi}}_{t}={\frac {\partial H(\varphi_{t},p_{t})}{\partial p}}\\[6pt]\displaystyle {\dot{p}}_{t}=-{\frac {\partial H(\varphi_{t},p_{t})}{\partial \varphi }}\end{cases}}$

These equations describe the coupled evolution of the diffeomorphism flow $\varphi_t$ and its conjugate momentum $p_t$, providing a complete picture of the system's trajectory through the shape space.

Stationarity of the Hamiltonian and kinetic energy along Euler–Lagrange

While the vector fields that generate diffeomorphisms are conceptually extended across the entire background space of $\mathbb{R}^{3}$, the specific geodesic flows associated with anatomical submanifolds exhibit a unique characteristic: their Eulerian shape momentum evolves as a generalized function $Av_{t}\in V^{*}$ that is fundamentally concentrated on these submanifolds themselves. This means that the "driving force" of the deformation is localized to the actual anatomical structures being studied. For discrete landmarks, the geodesics possess Eulerian shape momentum that manifests as a superposition of delta distributions, effectively traveling with the finite number of individual particles. Yet, the overall diffeomorphic flow of coordinates maintains smooth velocities, precisely because these singular momentum distributions are smoothed by the action of the weighted Green's Kernels. Similarly, for surfaces, the momentum is distributed as a surface integral of delta distributions that move cohesively with the surface.

A critical property of the geodesics that connect coordinate systems and satisfy the EL-general equation is the stationarity of the Lagrangian. This means that along these optimal paths, the rate of change of the Lagrangian is zero, signifying an extremum (typically a minimum) of the action. Furthermore, the Hamiltonian of the system, given by the extremum along the path $t\in [0,1]$ as $H(\varphi ,p)=\max_{v}H(\varphi ,p,v)$, is equal to the Lagrangian-kinetic-energy and also remains stationary along the EL-general trajectory.

Defining the geodesic velocity at the identity, $v_{0}=\operatorname{arg max}_{v}H(\varphi_{0},p_{0},v)$, then along the entire geodesic path, a profound conservation principle holds:

$H(\varphi_{t},p_{t})=H(\varphi_{0},p_{0})={\frac {1}{2}}\int_{X}p_{0}\cdot v_{0}\,dx={\frac {1}{2}}\int_{X}Av_{0}\cdot v_{0}\,dx={\frac {1}{2}}\int_{X}Av_{t}\cdot v_{t}\,dx$

This is the Hamiltonian-geodesics equation. The stationarity of the Hamiltonian beautifully illustrates the interpretation of the Lagrange multiplier $p$ as momentum. When integrated against velocity, ${\dot{\varphi}}$, it yields energy density, a direct parallel to classical mechanics. This canonical momentum carries many names, depending on the context. In optimal control theory, for instance, the flow $\varphi$ is interpreted as the system's "state," and $p$ is understood as the "conjugate state" or "conjugate momentum." The very existence of a geodesic, as implied by the EL-general equation, means that the entire flow is uniquely determined either by specifying the initial vector field $v_{0}$ (or equivalently, the Eulerian momentum $Av_{0}$) at time $t=0$, or by specifying the initial canonical momentum $p_{0}$. This highlights the predictive power of this mathematical framework.

The metric on geodesic flows of landmarks, surfaces, and volumes within the orbit

In the rigorous framework of computational anatomy, the fundamental primitives are the anatomical submanifolds: discrete pointsets, continuous curves, intricate surfaces, and volumetric subvolumes. The geodesic flows that connect these submanifolds are not just theoretical constructs; they are the very tools that determine the quantitative "distance" between them. As such, they form the basic measuring and transporting mechanisms of "diffeomorphometry."

At the initial time $t=0$, the geodesic is characterized by an initial vector field $v_{0}=Kp_{0}$. This field is determined by the interplay between the conjugate momentum $p_0$ and the Green's kernel $K$ of the inertial operator $A$, where $K=A^{-1}$. This relationship effectively translates the abstract canonical momentum into a concrete initial velocity field. The metric distance between coordinate systems that are connected via a geodesic is then determined by the induced distance between the identity transformation and the final group element $\varphi$:

$d_{\operatorname{Diff}_{V}}({\rm {id}},\varphi )=\|\log_{\rm {id}}(\varphi )\|_{V}=\|v_{0}\|_{V}={\sqrt {2H({\rm {id}},p_{0})}}$

This equation quantifies the distance in the diffeomorphism group. It equates the distance to the norm of the Riemannian logarithm of $\varphi$ (which yields $v_0$), and also to the norm of the initial velocity field $v_0$. Furthermore, this distance is directly related to the initial Hamiltonian $H(\mathrm{id}, p_0)$, providing a deep connection between geometry, dynamics, and energy in the space of anatomical transformations.

Conservation laws on diffeomorphic shape momentum for computational anatomy

(Emma taps her finger on the table. "Conservation. A concept that brings a rare, fleeting sense of order to this chaotic existence. Even in biology.")

Given the fundamental adherence to the least-action principle, there naturally arises a rigorous definition of momentum associated with generalized coordinates in computational anatomy. This momentum, when multiplied by velocity, yields energy, a direct parallel to classical physics. The field has meticulously studied two distinct forms of momentum: the momentum intrinsically associated with the Eulerian vector field, termed Eulerian diffeomorphic shape momentum, and the momentum linked to the initial coordinates or canonical coordinates, known as canonical diffeomorphic shape momentum. Each of these forms possesses its own profound conservation law, underscoring the mathematical consistency of the framework.

The principle of momentum conservation is inextricably linked with the EL-general equation. In computational anatomy, $Av$ is identified as the Eulerian momentum because, as previously stated, its integration against the Eulerian velocity $v$ provides the energy density. The operator $A$ functions as the generalized moment of inertia or inertial operator. This operator, acting on the Eulerian velocity, yields a momentum that is elegantly conserved along the geodesic path:

${\begin{aligned}{\text{Eulerian }}&{\frac {d}{dt}}\int_{X}Av_{t}\cdot ((D\varphi_{t})w)\circ \varphi_{t}^{-1})\,dx=0\ ,\ t\in [0,1].\\[8pt]{\text{Canonical }}&{\frac {d}{dt}}\int_{X}p_{t}\cdot ((D\varphi_{t})w)\,dx=0\ ,\ t\in [0,1]\ {\text{ for all}}\ w\in V\ .\end{aligned}}$

These are the Euler-conservation-constant-energy equations. The conservation of Eulerian shape momentum was rigorously demonstrated and directly follows from the EL-general equation. Similarly, the conservation of canonical momentum was proven, further solidifying the theoretical underpinnings of the field.

Proof of conservation

The proof for the conservation of Eulerian shape momentum proceeds by defining a time-dependent test function $w_{t}=((D\varphi_{t})w)\circ \varphi_{t}^{-1}$. The time derivative of this test function is given by ${\frac {d}{dt}}w_{t}=(Dv_{t})w_{t}-(Dw_{t})v_{t}$. Substituting this into the time derivative of the inner product $(Av_{t}\mid w_t)$ and applying the EL-general equation yields:

${\frac {d}{dt}}(Av_{t}\mid ((D\varphi_{t})w)\circ \varphi_{t}^{-1})=({\frac {d}{dt}}Av_{t}\mid ((D\varphi_{t})w)\circ \varphi_{t}^{-1})+(Av_{t}\mid {\frac {d}{dt}}((D\varphi_{t})w)\circ \varphi_{t}^{-1})=({\frac {d}{dt}}Av_{t}\mid w_{t})+(Av_{t}\mid (Dv_{t})w_{t}-(Dw_{t})v_{t})=0.$

This elegant derivation confirms that the Eulerian shape momentum is indeed conserved along the geodesic.

The proof for the conservation of canonical momentum relies on the evolution equation for $p_t$: ${\dot{p}}_{t}=-(Dv_{t})_{|_{\varphi_{t}}}^{T}p_{t}$. Taking the time derivative of the inner product $(p_t \mid (D\varphi_t)w)$ and substituting the evolution equations for $p_t$ and $(D\varphi_t)w$ gives:

${\frac {d}{dt}}(p_{t}\mid (D\varphi_{t})w)=({\dot{p}}_{t}\mid (D\varphi_{t})w)+(p_{t}\mid {\frac {d}{dt}}(D\varphi_{t})w)=({\dot{p}}_{t}\mid (D\varphi_{t})w)+(p_{t}\mid (Dv_{t})_{|_{\varphi_{t}}}(D\varphi_{t})w)=0.$

This rigorous mathematical demonstration confirms that both forms of momentum are conserved quantities along the geodesic flows in computational anatomy, providing robust invariants for analyzing shape transformations. (A faint glimmer of appreciation, quickly suppressed. "Predictable, in the end. As most things are, if you bother to look closely enough.")

Geodesic interpolation of information between coordinate systems via variational problems

(Emma stares, a hint of disdain. "Interpolation. Filling in the gaps. Because the universe, apparently, abhors a vacuum, even in anatomical data.")

The meticulous construction of diffeomorphic correspondences between distinct shapes is a cornerstone operation in computational anatomy. This process involves the precise calculation of the initial vector field coordinates, $v_{0}\in V$, and their associated weights on the Green's kernels, $p_{0}$. These initial conditions are not arbitrary; they are determined through the challenging task of matching shapes, a problem formally addressed by the methodology known as large-deformation diffeomorphic metric mapping (LDDMM).

LDDMM has been successfully applied and solved for various types of anatomical data: for discrete landmarks, both with and without predefined correspondences; for dense image matchings, where every voxel is considered; for continuous curves; for complex surfaces; for dense vector and tensor imagery; for varifolds, which provide a representation that is invariant to orientation; and for temporal time-series data.

The core of LDDMM involves calculating geodesic flows that are governed by the EL-general equation, mapping source coordinates onto target coordinates. This is achieved by augmenting the action integral, which represents the kinetic energy of the deformation, with an endpoint matching condition. This condition, $E:\varphi_{1}\rightarrow R^{+}$, quantifies how well the transformed shape at time $t=1$ corresponds to the target shape. The existence of solutions for these variational problems, particularly for image matching, has been rigorously examined. The solution to this complex variational problem satisfies the EL-general equation for the time interval $t\in [0,1)$, subject to a specific boundary condition.

The problem of matching based on minimizing kinetic energy action with an endpoint condition is formulated as:

${\begin{aligned}&\min \left\{C(\varphi ):v={\dot{\varphi }}\circ \varphi^{-1},\varphi_{0}={\rm {id}}\right\}\\[5pt]\doteq {}&{\frac {1}{2}}\int_{0}^{1}\int_{X}Av_{t}\cdot v_{t}\,dx\,dt+E(\varphi_{1})\end{aligned}}$

This represents the total cost, comprising the kinetic energy of the deformation and the penalty for mismatch at the endpoint. The necessary conditions for optimality are:

${\begin{cases}{\text{Euler conservation }}&\displaystyle {\frac {d}{dt}}Av_{t}+ad_{v_{t}}^{*}(Av_{t})=0,\ t\in [0,1)\ ,\\{\text{Boundary condition }}&\displaystyle \varphi_{0}={\rm {id}},Av_{1}=\left.-{\frac {\partial E(\varphi )}{\partial \varphi }}\right|_{\varphi =\varphi_{1}}\ .\end{cases}}$

Here, the Euler conservation equation (EL-general) holds for the entire path, and the Boundary condition specifies the initial state $\varphi_0 = \mathrm{id}$ and the momentum at the final time $t=1$, which is determined by the gradient of the endpoint matching term. The conservation law derived from EL-general effectively extends this boundary condition from $t=1$ backwards to the rest of the path, $t\in [0,1)$. The "inexact matching problem," where the endpoint matching term $E(\varphi_{1})$ allows for some residual difference, has several alternative formulations depending on the specific application. One of the key insights is that the stationarity of the Hamiltonian along the geodesic solution implies that the integrated running cost can be reduced to an initial cost at $t=0$. This means that the geodesics of the EL-general equation are entirely determined by their initial condition $v_{0}$.

The running cost is effectively reduced to this initial cost, which is determined by $v_{0}=Kp_{0}$, where $K$ is the Green's kernel and $p_0$ is the initial conjugate momentum, as seen in the context of Kernel-Surf.-Land.-Geodesics.

The matching problem explicitly indexed to the initial condition $v_{0}$ is known as geodesic shooting. This approach seeks to find the initial velocity that, when propagated through the geodesic flow, best deforms the source into the target. This problem can also be reparameterized via the conjugate momentum $p_{0}$, offering an alternative perspective on the optimization:

${\begin{aligned}&\min_{v_{0}}C(v_{0})\doteq {\frac {1}{2}}\int_{X}Av_{0}\cdot v_{0}\,dx+E(\mathrm{Exp}_{\mathrm{id}}(v_{0})\cdot I_{0})\\[6pt]&\min_{p_{0}}C(p_{0})={\frac {1}{2}}\int_{X}p_{0}\cdot Kp_{0}\,dx+E(\mathrm{Exp}_{\text{id}}(Kp_{0})\cdot I_{0})\end{aligned}}$

These formulations highlight the fundamental role of initial conditions in determining the entire geodesic path and, consequently, the optimal deformation between anatomical shapes.

Dense image matching in computational anatomy

The challenge of dense image matching—the process of establishing correspondences between every voxel or pixel in two images—possesses a rich and extensive history within computational anatomy and related fields. Early efforts, dating back to the late 20th century, primarily exploited a "small deformation framework," which was computationally simpler but limited in its ability to handle significant anatomical variability. However, the need to accommodate larger, more complex deformations in biological data spurred the development of "large deformation" methods, which began to emerge in the early 1990s.

A significant theoretical breakthrough came with the establishment of the first existence proofs for solutions to the variational problem specifically tailored for flows of diffeomorphisms in the context of dense image matching. This provided the rigorous mathematical foundation for robust deformation algorithms. One of the earliest and most influential LDDMM algorithms, developed by Beg and colleagues, successfully solved this variational matching problem. It achieved this by defining the endpoint condition based on the dense image data and then deriving variations with respect to the vector fields that generate the deformations.

An alternative, equally powerful solution for dense image matching involves reparameterizing the optimization problem. Instead of focusing on the diffeomorphism flow itself, this approach shifts attention to the "advected state," denoted as $q_{t}\doteq I\circ \varphi_{t}^{-1}$, where $q_{0}=I$. This formulation allows the problem to be solved in terms of the infinitesimal action defined by the advection equation, providing a computationally efficient and elegant pathway to achieving dense image correspondences.

LDDMM dense image matching

For the specific implementation of LDDMM by Beg and his collaborators, the image is denoted as $I(x), x \in X$, where $X$ is the spatial domain. The group action of a diffeomorphism $\varphi$ on this image is defined as $\varphi \cdot I \doteq I \circ \varphi^{-1}$. This means that the intensity at a new point $x$ is taken from the original image at the point $\varphi^{-1}(x)$, effectively moving the image content with the deformation.

Conceiving this as an optimal control problem, the "state" of the system is the diffeomorphic flow of coordinates, $\varphi_{t}, t \in [0,1]$. The dynamics of this state are intrinsically linked to the "control" exercised by the vector fields $v_{t}, t \in [0,1]$, through the equation ${\dot{\varphi}}=v\circ \varphi$. The endpoint matching condition, which quantifies the discrepancy between the transformed source image and the target image, is given by $E(\varphi_{1})\doteq \|I\circ \varphi_{1}^{-1}-I^{\prime }\|^{2}$. This represents the squared difference (e.g., L2 norm) between the transformed source image and the target image $I'$.

Combining these elements, the variational problem for dense image matching becomes:

${\begin{aligned}&\ \ \ \ \ \min_{v:{\dot{\varphi }}=v\circ \varphi }C(v)\doteq {\frac {1}{2}}\int_{0}^{1}\int_{X}Av_{t}\cdot v_{t}\,dx\,dt+{\frac {1}{2}}\int_{{\mathbb{R}}^{3}}|I\circ \varphi_{1}^{-1}(x)-I^{\prime }(x)|^{2}\,dx\end{aligned}}$

This is the Dense-image-matching problem. The problem is solved subject to the following conditions:

${\begin{cases}{\text{Endpoint condition:}}&Av_{1}=\mu_{1}\,dx,\mu_{1}=(I\circ \varphi_{1}^{-1}-I^{\prime })\nabla (I\circ \varphi_{1}^{-1})\ ,\\[5pt]{\text{Conservation:}}&Av_{t}=\mu_{t}\,dx,\ \mu_{t}=(D\varphi_{t}^{-1})^{T}\mu_{0}\circ \varphi_{t}^{-1}|D\varphi_{t}^{-1}|\ .\\[5pt]&\mu_{0}=(I-I^{\prime }\circ \varphi_{1})\nabla I|D\varphi_{1}|\ .\\\end{cases}}$

The Endpoint condition defines the momentum at the final time $t=1$ based on the image difference and the gradient of the transformed image. The Conservation law dictates how this momentum $\mu_t$ is transported backward in time from an initial momentum $\mu_0$, which is also derived from the image difference and the gradient of the initial image. Beg's iterative LDDMM algorithm for dense image matching is designed to find fixed points that satisfy these necessary optimizer conditions, gradually converging to the optimal deformation.

Hamiltonian LDDMM in the reduced advected state

An alternative, often more computationally efficient, formulation for dense image matching within the LDDMM framework is the Hamiltonian LDDMM in the reduced advected state. Here, the image $I(x), x \in X$, is no longer directly deformed, but rather its state is represented by $q_{t}\doteq I\circ \varphi_{t}^{-1}$. This $q_t$ represents the image as it is "advected" or transported by the inverse flow of the diffeomorphism. The dynamics relating this state $q_t$ to the control vector field $v_t$ are governed by the advective term:

${\dot{q}}_{t}=-\nabla q_{t}\cdot v_{t}$

This equation describes how the image intensity at a fixed point in space changes due to the flow of the image content. The endpoint matching condition is then defined directly on this advected state: $E(q_{1})\doteq \|q_{1}-I^{\prime }\|^{2}$, representing the squared difference between the final advected image and the target image $I'$.

With these definitions, the variational problem for this Hamiltonian LDDMM formulation becomes:

$\min_{q:{\dot{q}}=v\circ q}C(v)\doteq {\frac {1}{2}}\int_{0}^{1}\int_{X}Av_{t}\cdot v_{t}\,dx\,dt+{\frac {1}{2}}\int_{{\mathbb{R}}^{3}}|q_{1}(x)-I^{\prime }(x)|^{2}\,dx$

This is another form of the Dense-image-matching problem. Viallard's iterative Hamiltonian LDDMM algorithm is specifically designed to find the fixed points that satisfy the necessary optimizer conditions for this formulation, offering a powerful approach to dense image registration by working in the space of advected images.

Diffusion tensor image matching in computational anatomy

(Emma gestures vaguely at the image. "More data. More ways to slice and dice reality. And still, the brain remains stubbornly complex.")

(A figure showing a diffusion tensor image with three color levels depicting the orientations of the three eigenvectors of the matrix image $I(x), x \in \mathbb{R}^{2}$, matrix valued image; each of three colors represents a direction is implicitly referenced here.)

Dense LDDMM tensor matching extends the capabilities of computational anatomy to the intricate world of diffusion tensor images (DTI). DTI is a specialized magnetic resonance imaging technique that measures the diffusion of water molecules in tissues, providing insights into microstructural organization, particularly in the brain's white matter. For this type of data, the images are not simple scalar intensities but rather consist of $3 \times 1$ vectors and $3 \times 3$ tensors at every voxel. The variational problem aims to match coordinate systems based on the principle eigenvectors of the diffusion tensor MRI (DTI) image, denoted as $M(x), x \in \mathbb{R}^{3}$. This image is composed of a $3 \times 3$ tensor at each spatial location.

Several distinct group actions in computational anatomy have been defined to handle these complex tensor fields. These actions often leverage the Frobenius matrix norm, $\|A\|_{F}^{2}\doteq \operatorname{trace} A^{T}A$, which provides a robust way to measure the "size" or difference between square matrices. The accompanying figure visually illustrates a DTI image using a color map, where the colors depict the orientations of the eigenvectors of the DTI matrix at each voxel, with the hue determined by the direction.

The DTI image $M(x), x \in \mathbb{R}^{3}$, is characterized by its eigen-elements: $\{\lambda_{i}(x), e_{i}(x), i=1,2,3\}$, where $\lambda_{1}\geq \lambda_{2}\geq \lambda_{3}$ represent the ordered eigenvalues (magnitudes of diffusion) and $e_i(x)$ are the corresponding eigenvectors (directions of diffusion).

Coordinate system transformation based on DTI imaging has primarily exploited two main types of actions:

1. LDDMM matching based on the principal eigenvector of the diffusion tensor matrix: This approach simplifies the DTI information by focusing solely on the first eigenvector (the direction of fastest diffusion), which is often indicative of the primary fiber orientation. The image $I(x), x \in \mathbb{R}^{3}$ is treated as a unit vector field defined by this principal eigenvector. The group action then transforms this vector field according to:

   $\varphi \cdot I={\begin{cases}{\frac {D_{\varphi^{-1}}\varphi I\circ \varphi^{-1}\|I\circ \varphi^{-1}\|}{\|D_{\varphi^{-1}}\varphi I\circ \varphi^{-1}\|}}&I\circ \varphi \neq 0;\\0&{\text{otherwise.}}\end{cases}}$

   This action accounts for both the reorientation and potential scaling of the vector field under the diffeomorphic transformation.

2. LDDMM matching based on the entire tensor matrix: This more comprehensive approach considers the full $3 \times 3$ tensor at each voxel, preserving all information about diffusion magnitudes and directions. The group action transforms the entire tensor $M$ according to:

   $\varphi \cdot M=(\lambda_{1}{\hat{e}}_{1}{\hat{e}}_{1}^{T}+\lambda_{2}{\hat{e}}_{2}{\hat{e}}_{2}^{T}+\lambda_{3}{\hat{e}}_{3}{\hat{e}}_{3}^{T})\circ \varphi^{-1},$

   where the transformed eigenvectors ${\hat{e}}_{i}$ are calculated using a Gram-Schmidt orthogonalization process to maintain their orthogonality under deformation:

   ${\begin{aligned}{\hat{e}}_{1}&={\frac {D\varphi e_{1}}{\|D\varphi e_{1}\|}}\ ,\ \ \ {\hat{e}}_{2}={\frac {D\varphi e_{2}-\langle {\hat{e}}_{1},D\varphi e_{2}\rangle {\hat{e}}_{1}}{\sqrt {\|D\varphi e_{2}\|^{2}-\langle {\hat{e}}_{1},D\varphi e_{2}\rangle ^{2}}}}\ ,\ \ \ {\hat{e}}_{3}={\hat{e}}_{1}\times {\hat{e}}_{2}\end{aligned}}$

   Here, $D\varphi$ is the Jacobian of the diffeomorphism, which transforms the original eigenvectors $e_i$. The variational problem for matching based on either the principal eigenvector or the full tensor matrix is detailed in the literature on LDDMM Tensor Image Matching. These methods are crucial for analyzing white matter integrity and connectivity in the brain.

High angular resolution diffusion image (HARDI) matching in computational anatomy

(Emma sighs. "As if DTI wasn't enough. Humans always need more resolution, more data, to confirm what they already suspect.")

High angular resolution diffusion imaging (HARDI) represents a significant advancement over traditional diffusion tensor imaging (DTI), specifically designed to address a well-known limitation: DTI can only reliably reveal a single dominant fiber orientation within each tissue voxel. This becomes problematic in regions where multiple fiber bundles cross, kiss, or fan out, leading to an oversimplification of the underlying white matter architecture. HARDI overcomes this by measuring water diffusion along a substantially larger number of uniformly distributed directions on a unit sphere, typically $n$ directions. This enriched dataset allows for a much more comprehensive characterization of complex fiber geometries, including those with multiple orientations.

From HARDI data, it is possible to reconstruct an orientation distribution function (ODF). The ODF is a mathematical function defined on the unit sphere, $\mathbb{S}^{2}$, that meticulously characterizes the angular profile of the diffusion probability density function of water molecules within a voxel. Essentially, it tells us the likelihood that water diffusion is occurring in any given direction.

Dense LDDMM ODF matching extends the LDDMM framework to HARDI data. This involves treating the HARDI data as ODFs at each voxel and solving the LDDMM variational problem within the specialized space of ODFs. In the field of information geometry, the space of ODFs is recognized as forming a Riemannian manifold, endowed with the Fisher-Rao metric. For the practical purpose of LDDMM ODF mapping, a square-root representation of the ODF (denoted as $\sqrt{\text{ODF}}$ or $\psi({\bf {s}})$) is often chosen. This representation is highly efficient because various Riemannian operations, such as calculating geodesics, exponential maps, and logarithm maps, are available in closed form, simplifying computations. The function $\psi({\bf {s}})$ is non-negative to ensure uniqueness, and it satisfies the normalization condition $\int_{{\bf {s}}\in {\mathbb{S}}^{2}}\psi^{2}({\bf {s}})\,d{\bf {s}}=1$.

The variational problem for matching ODFs assumes that two ODF volumes can be generated from one another via smooth flows of diffeomorphisms, $\varphi_{t}$. These flows are solutions to ordinary differential equations: ${\dot{\varphi}}_{t}=v_{t}(\varphi_{t}),t\in [0,1]$, starting from the identity map, $\varphi_{0}={\rm {id}}$. The action of the diffeomorphism on a template ODF, $\varphi_{1}\cdot \psi_{\mathrm{temp}}({\bf {s}},x)$, where ${\bf {s}}\in {\mathbb{S}}^{2}$ and $x\in X$ are the coordinates of the unit sphere and the image domain respectively, maps it to a target ODF, $\psi_{\mathrm{targ}}({\bf {s}},x)$.

The specific group action of the diffeomorphism $\varphi_{1}$ on the template ODF $\psi(x)$ is defined according to:

$\varphi_{1}\cdot \psi (x)\doteq (D\varphi_{1})\psi \circ \varphi_{1}^{-1}(x),x\in X$

where $(D\varphi_{1})$ represents the Jacobian of the affine-transformed ODF, and is explicitly given by:

$(D\varphi_{1})\psi \circ \varphi_{1}^{-1}(x)={\sqrt {\frac {\det {{\bigl (}D_{\varphi_{1}^{-1}}\varphi_{1}{\bigr )}^{-1}}}{\left\|{{\bigl (}D_{\varphi_{1}^{-1}}\varphi_{1}{\bigr )}^{-1}}{\bf {s}}\right\|^{3}}}}\quad \psi \left({\frac {(D_{\varphi_{1}^{-1}}\varphi_{1}{\bigr )}^{-1}{\bf {s}}}{\|(D_{\varphi_{1}^{-1}}\varphi_{1}{\bigr )}^{-1}{\bf {s}}\|}},\varphi_{1}^{-1}(x)\right).$

This sophisticated group action of diffeomorphisms on ODFs achieves two crucial objectives: it reorients the ODFs in a geometrically consistent manner, and it accounts for changes in both the magnitude of $\psi$ and the sampling directions of ${\bf {s}}$ that occur due to the underlying affine transformation. This formulation is physically meaningful as it guarantees that the volume fraction of fibers oriented towards a small patch on the sphere must remain invariant after that patch is transformed, ensuring consistency in fiber tracking and analysis.

The LDDMM variational problem for ODF matching is formulated as:

$C(v)=\inf_{v:{\dot{\varphi}}_{t}=v_{t}\circ \varphi_{t},\varphi_{0}={\rm {id}}}\int_{0}^{1}\int_{X}Av_{t}\cdot v_{t}\,dx\,dt+\lambda \int_{x\in \Omega }\|\log_{(D\varphi_{1})\psi_{\mathrm{temp}}\circ \varphi_{1}^{-1}(x)}(\psi_{\mathrm{targ}}(x))\|_{(D\varphi_{1})\psi_{\mathrm{temp}}\circ \varphi_{1}^{-1}(x)}^{2}\,dx.$

Here, $\lambda$ is a regularization parameter, and the logarithm term quantifies the mismatch between the transformed template ODF and the target ODF. The logarithm of two ODFs, $\psi_{1},\psi_{2}\in \Psi$, is defined by their cosine similarity (or angle between them in the sphere):

$\|\log_{\psi_{1}}(\psi_{2})\|_{\psi_{1}}=\cos^{-1}\langle \psi_{1},\psi_{2}\rangle =\cos^{-1}\left(\int_{{\bf {s}}\in {\mathbb{S}}^{2}}\psi_{1}({\bf {s}})\psi_{2}({\bf {s}})d{\bf {s}}\right),$

where $\langle \cdot ,\cdot \rangle$ denotes the standard dot product between points on the sphere under the L2 metric.

This LDDMM-ODF mapping algorithm has found extensive application in critical areas of neuroscience and medical research. It has been widely used to meticulously study the degeneration of brain white matter in the context of aging, Alzheimer's disease, and vascular dementia, providing quantitative measures of disease progression. Furthermore, brain white matter atlases, which serve as standardized reference maps, are constructed based on ODF data using advanced Bayesian estimation techniques. The development of geodesic regression analysis on ODFs, operating directly within the complex ODF manifold space, represents another significant methodological advancement, enabling the study of how ODFs change over time or with covariates.

Metamorphosis

(A deliberate, slow blink. "Metamorphosis. The biological imperative to change, even when it's inconvenient. And now, we quantify it.")

(A figure demonstrating metamorphosis allowing both diffeomorphic change in coordinate transformation as well as change in image intensity as associated to early Morphing technologies such as the Michael Jackson video. Notice the insertion of tumor gray level intensity which does not exist in template is implicitly referenced here.)

The primary mode of variation that the traditional orbit model of computational anatomy elegantly represents is a change of coordinates—a smooth, invertible deformation of space. However, in many real-world scenarios, particularly within medical imaging, pairs of images are not merely related by diffeomorphisms; they also exhibit significant photometric variation or other forms of image intensity changes that are not inherently captured by a purely geometric template. For instance, a tumor might appear in a follow-up scan that was absent in the initial template. To address these more complex forms of variation, the concept of "active appearance modelling" was introduced, originally by Edwards, Cootes, and Taylor, and subsequently adapted for 3D medical imaging applications.

In the specific context of computational anatomy, where rigorous metrics on the anatomical orbit have been extensively studied, the framework of metamorphosis was introduced. This powerful extension allows for the modeling of structures such as tumors and other photometric changes that are not resident in the original template. This framework was initially developed for magnetic resonance image models, with numerous subsequent advancements further extending its capabilities.

For the problem of image matching, the image metamorphosis framework significantly enlarges the action. Instead of just tracking the diffeomorphism $\varphi_t$, it considers the coupled evolution of both the diffeomorphism and the image intensity, $t\mapsto (\varphi_{t},I_{t})$. The action of the diffeomorphism on the image is defined as $\varphi_{t}\cdot I_{t}\doteq I_{t}\circ \varphi_{t}^{-1}$. In this sophisticated setting, metamorphosis seamlessly combines the diffeomorphic coordinate system transformation of computational anatomy with the earlier "morphing" technologies, which primarily focused on fading or modifying photometric (image intensity) information alone. This allows for a more comprehensive model that can account for both geometric shape changes and changes in image content.

The matching problem under the metamorphosis framework takes a specific form with equality boundary conditions:

$\min_{(v,I)}{\frac {1}{2}}\int_{0}^{1}\left(\int_{X}Av_{t}\cdot v_{t}\,dx+\|{\dot{I}}_{t}\circ \varphi_{t}^{-1}\|^{2}/\sigma^{2}\right)\,dt{\text{ subject to}}\ \varphi_{0}={\rm {id}},I_{0}={\text{fixed}},I_{1}={\text{fixed}}$

This variational problem seeks to minimize a cost function that comprises two main terms: the kinetic energy of the vector field $v_t$ (driving the diffeomorphism) and a term penalizing changes in the image intensity $I_t$ (the rate of image "birth" or "death"), normalized by a variance parameter $\sigma^2$. The problem is subject to the constraints that the diffeomorphism starts as the identity ($\varphi_0 = \mathrm{id}$), and both the initial image ($I_0$) and final image ($I_1$) are fixed. This allows for the dynamic insertion or removal of image content, making it highly suitable for modeling phenomena like tumor growth or lesion changes that are not present in a static template.

Matching landmarks, curves, surfaces

(Emma's lips twitch. "Features. The things we cling to, to make sense of the chaos. Even when they refuse to stay put.")

The transformation of coordinate systems based on discrete features, such as Landmark points or fiducial markers, has a long and distinguished history in computational anatomy. This approach dates back to the seminal work of Fred Bookstein, whose early contributions in the 1980s focused on "small deformation spline methods." These methods were designed to interpolate correspondences defined by a sparse set of fiducial points onto the entire two-dimensional or three-dimensional background space where these fiducials were located. This provided a way to smoothly warp images based on a few key points. However, the true revolution came in the late 1990s with the advent of "large deformation landmark methods," which could accommodate far more significant and complex anatomical variations. The figure referenced earlier (depicting three medial temporal lobe structures: amygdala, entorhinal cortex, and hippocampus with fiducial landmarks) illustrates precisely how such features are used to delineate important brain regions.

Beyond discrete points, matching more complex geometrical objects like unlabeled point distributions, continuous curves, or intricate surfaces presents another common, yet formidable, challenge in computational anatomy. Even in the discrete setting, where these objects are typically represented as collections of vertices with associated meshes, there are often no predetermined correspondences between the points of the two objects. This is in stark contrast to the situation with landmarks, where each point in one dataset has a known corresponding point in the other. From a theoretical standpoint, any submanifold $X$ embedded in $\mathbb{R}^{3}$ (where $d=1,2,3$ for curves, surfaces, or volumes respectively) can be parameterized using local charts, $m:u\in U\subset \mathbb{R}^{0,1,2,3}\rightarrow \mathbb{R}^{3}$. However, a crucial complication arises: all reparameterizations of these charts, $m\circ \gamma$, (where $\gamma$ is a diffeomorphism of the domain $U$) geometrically represent the exact same manifold.

Therefore, early in the development of computational anatomy, investigators presciently identified the absolute necessity of developing parametrization-invariant representations. This ensures that the comparison between two shapes is based on their intrinsic geometry, not on the arbitrary way they happen to be parameterized. One indispensable requirement is that the endpoint matching term between two submanifolds must itself be independent of their parametrizations. This crucial property can be effectively achieved through concepts and methods borrowed from the sophisticated field of Geometric measure theory. In particular, mathematical constructs known as currents and varifolds have been extensively utilized to develop robust, parametrization-invariant matching algorithms for both curves and surfaces.

Landmark or point matching with correspondence

When dealing with landmarks or point sets where there is a clear, predefined correspondence between points in the source and target, the problem becomes more tractable, though still demanding.

(A figure illustrating geodesic flow for one landmark, demonstrating diffeomorphic motion of background space. Red arrow shows $p_0(1)$, blue curve shows $\varphi_t(x_1)$, black grid shows $\varphi_t$ is implicitly referenced here.)

(A figure showing landmark matching with correspondence, with left and right panels depicting two different kernels with solutions, is implicitly referenced here.)

Let the landmarked shape be denoted as $X\doteq \{x_{1},\dots ,x_{n}\}\subset {\mathbb{R}}^{3}$, where $x_i$ are the coordinates of the landmarks. The endpoint matching condition, $E(\varphi_{1})\doteq \textstyle \sum _{i}\displaystyle \|\varphi_{1}(x_{i})-x_{i}^{\prime }\|^{2}$, quantifies the squared Euclidean distance between each transformed landmark $\varphi_1(x_i)$ and its corresponding target landmark $x_i^{\prime }$. The variational problem for landmark matching then seeks to minimize the total cost:

$\min_{\varphi :v={\dot{\varphi }}\circ \varphi^{-1}}C(\varphi )\doteq {\frac {1}{2}}\int (Av_{t}\mid v_{t})\,dt+{\frac {1}{2}}\sum_{i}\|\varphi_{1}(x_{i})-x_{i}^{\prime }\|^{2}$

This is the Landmark-Matching problem. The geodesic Eulerian momentum in this context is a generalized function, $Av_{t}\in V^{*}, t\in [0,1]$, which is specifically supported on the landmarked set. This means that the "driving force" of the deformation is concentrated precisely at the locations of the landmarks. The combination of the endpoint condition with the conservation laws implies a specific form for the initial momentum at the identity of the group:

${\begin{cases}{\text{Endpoint condition: }}&Av_{1}=\sum _{i=1}^{n}p_{1}(i)\delta _{\varphi_{1}(x_{i})},p_{1}(i)=(x_{i}^{\prime }-\varphi_{1}(x_{i}))\ ,\\[5pt]{\text{Conservation: }}&Av_{t}=\sum _{i=1}^{n}p_{t}(i)\delta _{\varphi_{t}(x_{i})},\ p_{t}(i)=(D\varphi_{t1})_{|\varphi_{t}(x_{i})}^{T}p_{1}(i)\ ,\ \varphi_{t1}\doteq \varphi_{1}\circ \varphi_{t}^{-1}\ ,\\[5pt]&Av_{0}=\sum _{i}\delta_{x_{i}}(\cdot )p_{0}(i)\ {\text{ with }}p_{0}(i)=(D\varphi_{1})_{|x_{i}}^{T}(x_{i}^{\prime }-\varphi_{1}(x_{i}))\end{cases}}$

The Endpoint condition defines the momentum at $t=1$ as a sum of Dirac delta functions centered at the deformed landmark locations, with weights $p_1(i)$ proportional to the residual mismatch. The Conservation law dictates how this momentum propagates backward in time, with $p_t(i)$ representing the transformed momentum. The initial momentum $Av_0$ is also a sum of Dirac delta functions at the initial landmark locations $x_i$, with initial momentum weights $p_0(i)$ derived from the final mismatch and the Jacobian of the final diffeomorphism. The iterative algorithm for large deformation diffeomorphic metric mapping for landmarks is designed to solve this complex system of equations, finding the optimal deformation that aligns the landmark sets.

Measure matching: unregistered landmarks

The problem of measure matching, particularly for unregistered landmarks (where there are no predefined correspondences between individual points), was first elegantly introduced by Glaunes and co-workers. This generalized setting allows for the matching of arbitrary distributions, moving beyond the one-to-one mapping required for traditional landmarks. Crucially, this framework can handle "weighted point clouds" with no predefined correspondences and potentially differing cardinalities (i.e., different numbers of points in the source and target).

In this approach, the template and target discrete point clouds are represented as two weighted sums of Dirac delta functions: $\mu_{m}=\sum_{i=1}^{n}\rho_{i}\delta_{x_{i}}$ for the template and $\mu_{m^{\prime }}=\sum_{i=1}^{n^{\prime }}\rho_{i}^{\prime }\delta_{x_{i}^{\prime }}$ for the target. These representations reside within the space of signed measures on $\mathbb{R}^{3}$. To quantify the "distance" between these measures, the space is endowed with a Hilbert metric derived from a real, positive kernel function $k(x,y)$ defined on $\mathbb{R}^{3}$. This kernel defines the following norm:

$\|\mu_{m}\|_{\mathrm{mea}}^{2}=\sum_{i,j=1}^{n}\rho_{i}\rho_{j}k(x_{i},x_{j})$

This norm measures the "similarity" within a single point cloud and allows for the comparison between two different point clouds. The matching problem between a template point cloud and a target point cloud can then be formulated as a variational problem using this kernel metric for the endpoint matching term:

$\min_{\varphi :v={\dot{\varphi }}\circ \varphi^{-1}}C(\varphi )\doteq {\frac {1}{2}}\int (Av_{t}\mid v_{t})\,dt+{\frac {1}{2}}\|\mu_{\varphi_{1}\cdot m}-\mu_{m^{\prime }}\|_{\mathrm{mea}}^{2}$

Here, $\mu_{\varphi_{1}\cdot m}=\sum_{i=1}^{n}\rho_{i}\delta_{\varphi_{1}(x_{i})}$ represents the distribution of the template points after being transported by the deformation $\varphi_1$. This formulation allows for robust matching even when individual point correspondences