Generation of Tubular and Membranous Shape Textures with Curvature Functionals

Tubular and membranous shapes display a wide range of morphologies that are difficult to analyze within a common framework. By generalizing the classical Helfrich energy of biomembranes, we model them as solutions to a curvature optimization problem in which the principal curvatures may play asymmetric roles. We then give a novel phase-field formulation to approximate this geometric problem, and study its Gamma-limsup convergence. This results in an efficient GPU algorithm that we validate on well-known minimizers of the Willmore energy; the software for the implementation of our algorithm is freely available online. Exploring the space of parameters reveals that this comprehensive framework leads to a wide continuum of shape textures. This first step towards a unifying theory will have several implications, in biology for quantifying tubular shapes or designing bio-mimetic scaffolds, but also in computer graphics, materials science, or architecture.


Introduction
Tubular and membranous shape textures are widely present in biology.They display a large variety of morphologies in terms of geometry and topology, which are important to analyze since they reflect the state of a biological system.For instance, the bone marrow capillaries are highly branching and merging vessels [95,83], whose organization is subject to drastic remodeling in acute myeloid leukaemia [77,27].In cells, the endoplasmic reticulum, where proteins are synthesized, consists of an interconnected network [90] that undergoes sheets-to-tubules topological transformations [82].Furthermore, trabecular bone is a combination of rods and platelets [71,88] that are optimally restructured under mechanical stress [85,1,87] or pathological conditions [76,97,37].
However, due to their disparity and complexity, tubular and membranous structures are difficult to describe within a unifying framework that captures both their rich morphological diversity as well as their continuous variations.We approach this question by building a generation model that creates shape textures from noise, similarly to texture synthesis in images [43,80,58].We model tubules and membranes as optimizers under constant volume of a curvature functional where p is a second-degree polynomial of the principal curvatures κ 1 and κ 2 of the surface S. As our main contribution, we provide a novel phase-field formulation F to approximate the original geometric problem F, and show that the Γ-limsup holds, a notion coming from the Γconvergence framework [21,2,9].The optimization problem then translates into the mass-preserving H −1 gradient flow [32,19] u = ∆ ∂F ∂u . ( Combining the stochastic optimizer Adam [52,61] to the automatic differentiation provided by PyTorch [78] results in an efficient and flexible GPU implementation, curvatubes.It successfully leads to a wide continuum of shape textures (see Figure 1), which constitutes a first step towards a unifying theory.
Related work.The variational formulation (1) generalizes two classical functionals, the Willmore energy studied in differential geometry [94,103,4,62,48,11,98], as well as the Helfrich energy used to model biomembranes [15,41,22,25,91,14].In contrast to these models, the polynomial p is here not required to be symmetric Figure 1: A continuum of shape textures generated by curvatubes, after optimizing a curvature-based polynomial energy of the surface S p(κ1, κ2) dA, with a volume constraint.The polynomial coefficients vary linearly in space, by interpolating four values that define different shape textures, at the vertical median of the four squares.In the formulas, κ1 and κ2 are the principal curvatures of the suface, while H = κ1 + κ2 and K = κ1κ2 are the mean and the Gaussian curvatures.More details are given at the end of Section 5.2.
in the principal curvatures, which allows the generation of tubules.Other generalizations have been proposed in [100,40] and [23,16].
Curvature functionals are often used as image prior models in imaging, due to their ability to interpolate.Mumford [70] considered Euler's elastica functional as a prior curve model in computer vision, and this was subsequently applied to digital inpainting [63,8,93].Similar ideas were then used for 3D volume reconstruction from 2D slices [60,10,50].Our work is also related to Poisson reconstruction methods [45,46] that implicitly reconstruct a surface from noisy oriented points.
However, as in [25,26,24], we use the diffuse setting to model shapes.By construction, our phase-field formulation F extends the standard approximation of the Willmore energy (see [11] and references therein) as well as the one proposed by Bellettini & Mugnai in [6] for the Helfrich energy.In terms of Γ-convergence, we also provide an extension of the Γ-limsup result in [6], while keeping the Γliminf as an open question (except in the cases previously covered) that we do not aim to solve here.
Finally, our attempt in building a unifying framework is connected to the Functionalized Cahn-Hilliard (FCH) model proposed in [38,20,17,53,54,18].The FCH energy describes how amphiphilic molecules self-assemble into complex network morphologies that feature spheres, tubules, sheets, and mixtures of them.Loosely speaking, such shape textures result from a compromise between minimizing a bending energy while rewarding an increase in interfacial area.Yet, we choose to treat the principal cur-vatures independently, which is not the case in their model.
Outline.In Section 2, we introduce the Willmore and the Helfrich energies, as well as their classical phase-field approximations in light of the Γ-convergence framework.In section 3, we build the phase-field functional and justify its construction with heuristic arguments and in terms of Γlimsup.Section 4 describes the computational framework and displays a large number of simulations.Finally, we discuss in Section 5 the implications of a unifying theory on applied fields and present future extensions to this work.

Background
We begin by introducing the Willmore and the Helfrich energies in more detail, then sketch the notion of Γconvergence before moving on with their classical phasefield approximations.
But first, let us describe curvatures in simple terms.Curvature measures how much a line or a surface is locally deviating from a straight line or a flat plane.For a line contained in a plane, the curvature κ at a point is the inverse1 r of the radius of the osculating circle, i.e., the tangent circle that approaches the curves most tightly at this point.Surfaces are characterized at each point by two principal curvatures κ 1 and κ 2 , that correspond to the maximal and minimal curvatures of the lines resulting from the perpendicular intersection of the surface with a plane.The maximally and minimally curved lines are directed by two perpendicular principal directions, tangent to the surface.For instance, the principal curvatures are zero on a plane; on a sphere of radius r, they are all equal to ( 1 r , 1 r ); on a cylinder of base radius R, they are all equal to ( 1 R , 0).

2.1
The Willmore, Helfrich, and generalized curvature functionals In the 1970's, Canham [15], and subsequently Helfrich [41], proposed to model the surface of a biomembrane as a minimizer of a curvature bending energy, or Helfrich energy In this expression, S is a smooth compact orientable surface in R 3 whose principal curvatures are denoted κ 1 and κ 2 , with the convention κ 1 ≥ κ 2 .The mean curvature H = κ 1 + κ 2 and the Gaussian curvature K = κ 1 κ 2 are respectively the sum and the product of the principal curvatures.The signs of H, κ 1 , and κ 2 depend on the orientation of S, while that of K does not.Throughout this work, we use the convention that H should be positive on convex objects like spheres.Let us remark that there is a closed relationship between (κ 1 , κ 2 ) and (H, K), provided by the bijection φ(x, y) In (3), the coefficients χ b > 0 and χ G are the bending and Gaussian rigidities.The parameter H 0 , or spontaneous curvature, models the asymmetry of the two layers composing the bilipidic membrane [92,29].The Willmore energy, defined as is then a special case of the Helfrich energy, with χ G = 0, χ b = 2, and H 0 = 0.In both of the classical functionals (3) and ( 5), the integrand is a polynomial p(κ 1 , κ 2 ) symmetric in the principal curvatures, i.e., p(κ 1 , κ 2 ) = p(κ 2 , κ 1 ), since it can be reformulated as a polynomial of their sum H and product K.However, one may want to construct a general curvature functional where p is a smooth function with no symmetry constraint.A similar 1 form is given in [100,98], for a smooth function q(H, K).A generalization to functions that depend on the position and the normal to the surface can be found in [23,16].
We will restrict ourselves to p which are polynomials of degree 2, as the framework is then rich enough to generate complex shape textures.We are thus interested in the curvature functional where we use multi-index notation.

Phase-fields and Γ-convergence
Numerically, critical points of the curvature energies (3), (5), and ( 6) can be searched for by following a gradient flow.Before implementing an algorithm, it is however preferable to convert these sharp-interface functionals defined for 2D surfaces, to diffuse approximations defined for scalar fields in a 3D volume.This way, surfaces are implicitly represented as level sets of the volumetric scalar field, which allows us to address topological changes encountered in the flow seamlessly; whereas in explicit methods, surfaces are tracked as a mesh that needs dynamic remeshing to avoid entanglement through topological transitions [75].
The gain is considerable in view of the high topological complexity of the targeted shape textures.This leads us to consider a phase-field approximating the original geometric problem, for instance as in [25,26,24] for the modeling of biomembranes.Phase-fields have been extensively used to model phase separation in binary mixtures, beginning with the Cahn-Hilliard (or Ginzburg-Landau) energy [13,31], that was subsequently reused in several other contexts [51,67].These functions typically take values close to 1 and −1 inside and outside a region, with a smooth transition between the two phases at the interface.A parameter > 0 represents the thickness of phase transition at the interface.
The quality of a diffuse approximation with phasefields is typically studied in the Γ-convergence framework [21,2,9].Γ-convergence expresses the convergence of minimization problems, so that, rather than solving a limit problem, we solve a sequence of approaching problems (or the reverse).The complete convergence consists in a Γlimsup, which relies on a constructive proof, and a Γ-liminf, generally more difficult to prove.Definition 1 (Γ-convergence).Given X a metric space, let F and F be functions from X to [−∞, +∞], where the F are defined for > 0. We say that the sequence (F ) >0 Γ-converges to F at a point u ∈ X as → 0, and write if the following two bounds hold: [Γ-liminf ] For every sequence (u ) such that u → u in X, [Γ-limsup] There exists a sequence (u ), called recovery sequence, such that u → u in X and This punctual definition can be extended to a convergence taking place on the whole space.Γ-convergence is especially interesting due to the following fundamental result.
Theorem 2. Let X be a metric space, and let F = Γ(X) − lim →0 F .Suppose that the sequence F : X → [−∞, +∞] is equi-coercive, i.e., for all t ∈ R there exists a compact set K t ⊂ X such that {F ≤ t} ⊂ K t .Then F admits a minimum and min Futhermore, if u minimizes F over X, then every cluster point of (u ) minimizes F over X.
This ensures not only the convergence of the minimal values, but also of the minimizers themselves.

Phase-field approximations of the area, the Willmore, and the Helfrich functionals
In this paragraph, we present three classical diffuse approximations that are important to our development.We fix some mathematical notations beforehand.Let e be a fixed unit-norm vector in R 3 .We consider a symmetric double-well function W (s) = 1  4 (1 − s 2 ) 2 that cancels on −1 and +1.Note that its derivative is

Notations
2 denote a constant that only depends on the double-well.
For a function u twice (weakly) differentiable, we define the normal vector field on the set {∇u = 0} e elsewhere, that has unit norm, and is orthogonal to the level sets of u.We introduce the matrix field whose trace is equal to The classical approximations.The Cahn-Hilliard phase-field is known to approximate the area (or perimeter) functional, which measures the total area of surfaces in the 3D space.More precisely, let us introduce and the area functional Following a conjecture of De Giorgi, Modica and Mortola [68] proved the Γ-convergence This means that, if E ⊂ Ω is such that S = ∂E ∩ Ω is smooth and of finite area, and setting u = 2χ E − 1 ∈ BV (Ω, {−1, 1}), then the Γ-limsup provides a sequence of functions (u ) ∈ W 1,2 (Ω) such that u → u in L 1 (Ω) and A (u ) → σ area(S), i.e., their diffuse areas converge to the area of S up to a factor σ.
Subsequently, several authors [7,99,5,69,86,73] studied diffuse approximations of the Willmore energy (5).Bellettini and Paolini [7] introduced the phase-fields +∞ otherwise in L 1 (Ω) (11) Note that the trace term TrM u = − ∆u+ W (u) inside the square is the L 2 gradient of 2 |∇u| 2 + W (u) which appears in the Cahn-Hilliard phase-field (8), in the same way as the mean curvature vector is the L 2 gradient of the area functional.
The Γ-limsup was showed in [7], using the same recovery sequence as for the area functional.The Γ-liminf was studied under several conditions in [5,69] and completed in [86].Together with the Γ-limsup, this resulted in the Γ-convergence on smooth points of the form u = 2χ E − 1 where E ⊂ Ω and but with the additional assumption that the diffuse surface areas A remain uniformly bounded.
Finally, Bellettini and Mugnai [6] extended the Willmore phase-field to approximate the complete Helfrich energy with +∞ otherwise in L 1 (Ω).( 12) Based on the previous results of Röger and Schätzle [86], and under the assumptions H 0 = 0 and −χ b < χ G < 0, they showed that the Γ-convergence holds on smooth points again using an additional uniform bound on the diffuse areas A .
Our aim is precisely to generalize the Helfrich phase-field formula H further to approximate the curvature functional (6), and provide a computational framework to simulate shape textures.We are now ready to construct a new phase-field, for which we will study the Γ-limsup property.

Construction of the phase-field
In this section, we generalize the Helfrich phase-field H in (12) to approximate the functional F in (6), using the notations introduced in Section 2.3.We justify the construction with heuristic arguments, and show that the Γ-limsup is still satisfied, although we do not attempt to show the Γ-liminf.

Diffuse curvatures and second fundamental form
Let us notice that the diffuse expressions W and H in ( 11) and ( 12) both rely on the trace and the norm of the matrix field M u introduced in ( 7), which is related to the second fundamental forms and the curvatures of the level sets of u as follows.We define the diffuse second fundamental form B u as well as the diffuse mean and Gaussian curvatures H u and K u using if ∇u = 0, and zero otherwise.Informally, B u ⊗ n u at the point x ∈ Ω approximates the second fundamental form of the level surface {u = u(x)} (well-defined if ∇u = 0 on this set).H u n u approximates the mean curvature vector, with the convention that it points inwards for convex sets, and K u approximates the Gaussian curvature.
Based on the relations ( 4) linking (κ 1 , κ 2 ) to (H, K), we also introduce the diffuse principal curvatures where we use the positive part x + = max(0, x).
It can be shown, using the implicit formulas summarized in [39], that the expressions from ( 14) to (18) coincide exactly with the second fundamental form and the respective curvatures of the level sets of u, in the special case where the function has a hyperbolic tangent profile where E ⊂ Ω is an open set with smooth boundary ∂E ∩ Ω ∈ C2 , and the signed distance from ∂E, denoted by dist ∂E , is by convention positive on E and negative on Ω \ Ē.In diffuse approximations, the tanh profile is optimal 2 and is generally used to construct the Γ-limsup recovery sequence.
The presence of the positive part in ( 17) and ( 18) ensures that the square root term is still defined when (H u ) 2 − 4K u < 0. This can happen, since for a general u, we have , and the numerator is a 2 + b 2 + c 2 − 2 (ab + bc + ac) which possibly has negative values, where a, b, c are the real eigenvalues of M u .However, if one of them is 0, the numerator is a squared difference and the positive part is not useful.This is the case in particular for functions with tanh profile.
The expressions ( 17) and ( 18) can be reformulated with M u , by writing if ∇u = 0, and zero otherwise.Note that H u can be retrieved from their sum

Phase-field construction
Let a = (a 2,0 , a 1,1 , a 0,2 , a 1,0 , a 0,1 , a 0,0 ) = (a α ) |α|≤2 ∈ R 6 be a vector of real coefficients.The associated polynomial function is denoted by p(x, y) = |α|≤2 a α (x, y) α , so that F(S) = S p(κ 1 , κ 2 ) dA.Consider the following expression, The heuristic intuition behind is that, if u has a tanh profile with transition parameter as in (19), we can apply the coarea formula to obtain where we use |∇u| . This amounts to integrating the curvature functional F over all the level surfaces S t = {u = t} of the phase-field u, appropriately weighted so that the largest contributions are given by level sets close to {u = 0}.As goes to zero, the level sets {u = t} concentrate around ∂E.
Still under the ansatz of tanh profile, E can be developed in terms of the matrix field M u : For a general u ∈ W 2,2 (Ω), the expressions of E (u) and which is satisfied for functions with tanh profile.
Finally, E (u) can be simplified further, by replacing the positive part in the first term outside the square root directly by 2 M u 2 − (TrM u ) 2 .As said earlier, this is true of the special tanh case, where 2 2 , where a, b and 0 are the eigenvalues of M u .This leads to the final form, defined for any u ∈ W 2,2 (Ω), This phase-field is devised to be a diffuse approximation of the sharp-interface functional F (6), up to the multiplicative factor σ.
Comparison with the Willmore and the Helfrich diffuse approximations.It can be checked that the proposed formulation F is indeed a generalization of the previous formulas W and H , by specifying the polynomial coefficients of the Willmore energy, and of the Helfrich energy,

Γ-limsup property
As explained in Section 2.2, the approximation of the target functional F by the sequence (F ) >0 can be studied in the Γ-convergence framework.Here, we assert that the phasefield F satisfies the Γ-limsup property, thus extending the result of [6].
The functionals F and F are defined as in (20) and (6).Then there exists a sequence of functionals (u The proof is given in the Appendix.It consists in showing that the recovery sequence constructed in [6] still satisfies the theorem for our more general formulation F .The first and third properties correspond to the existence of a recovery sequence.The second property loosely means that the measure whose density is |∇u | 2 in the 3D volumetric space concentrates into the measure induced by the area on the 2D surface.This intuition is in accordance with the way we constructed F (see beginning of Section 3.2).
However, whether or not the Γ-liminf holds still remains an open question, except in the special cases of the Willmore and the Helfrich energies with additional assumptions, as seen in Section 2.3.Yet, we believe that, even if the Γ-convergence could fail in general, this does not constitute a serious impediment to our phase-field expression being still of interest for generating shape textures.

Simulations
In this section, we demonstrate that the phase-field F constructed in the previous section can generate a large range of shape textures.We describe curvatubes in Algorithm 1, and then show the results of four numerical experiments.The first one validates the approach by finding well-known Willmore minimizers.We then display a gallery of 10 shape textures.The effect of smoothly varying the generation parameters is shown in the third experiment, with a bilinear interpolation between 4 shape textures, layers, spheres, tubes, and sponges.Finally, 1000 shape textures are generated with random parameters and visualized in an atlas with UMAP.The numerical codes are fully available at https://github.com/annasongmaths/curvatubes .

Curvatubes
Shape textures are generated by minimizing the phase-field energy F (20) under a constraint of constant volume, with periodic boundary conditions.More exactly, given a random initialization of the phase-field u, we find a point of convergence with low energy of the so-called The H −1 flow is mass-preserving, i.e., keeps constant the mass of u, denoted by ū := 1 Ω Ω u dx.The preservation of mass approximately encodes a constraint of constant volume on the region enclosed by the surface {u = 0}, if the phase-field u ±1 is nearly constant inside and outside.
The H −1 flow can actually be expressed as a standard L 2 flow, by relying on the change of variable where A : Ω → R3 is a periodic vector field, and m 0 ∈ R is the desired value of the average ū.We then define an energy with respect to A, It can be checked that in such a way that a L 2 flow on A becomes a H −1 flow on u: Ȧ = − ∂G ∂A ⇒ u = ∆ ∂F ∂u (provided that the derivatives in time and space of A commute).Therefore, the H −1 flow on u starting at u 0 = ∇ • A 0 + m 0 can be solved as a usual L 2 flow on A.
To generate shape textures, the variable A is initialized as a random white noise vector field A 0 and we reach a point of convergence of the L 2 flow Ȧ = − ∂G ∂A with Adam [52,61], a gradient-based stochastic optimization algorithm.The change of variable (24) allows us to benefit from the computation of the L 2 gradient ∂G ∂A by the automatic differentiation engine provided by PyTorch [78], combined with the efficiency of Adam.
The generation model is summarized in Algorithm 1.It takes as inputs the initialization A 0 , the coefficients a and the mass m 0 .After convergence, the output shape texture is defined as the level surface {u = 0} of the final phasefield u = ∇•A+m 0 .We color it in beige, and show the level set {u = 0.05} in dark red to enhance the visualization.An example of flow is given in Figure 2, with the corresponding loss curves in Figure 3.

8:
Update moments of the gradients Adam S ← {u = 0} visualize the shape texture Implementation details.The domain Ω is assimilated to a grid of size 100×100×100 pixels with a fixed sampling step ∆x = 0.01.We take = 0.02, unless specified otherwise.The phase-field u and the vector field A are encoded as matrices whose coefficients specify the sampled values.
The discrete energies F (resp.G ) are symbolically defined by a succession of elementary operations on u (resp.A), before being differentiated automatically by PyTorch.In particular, the integral is encoded as a finite sum, while the differential operations ∇u, Hess u, and ∇ • A, are computed as classical finite differences that take into account the periodicity of the problem.To prevent the formation of artifacts, we apply a Gaussian blur k with a small deviation (typically σ k = 2 pixels) to the phase-field u, before computing the finite differences.The norm of the gradient |∇u| is modified by a small offset ξ = 10 −6 , as in |∇u| 2 + ξ 2 or 1/ |∇u| 2 + ξ 2 , to prevent non-differentiability at zero and division by zero.The positive part function x + appearing in F is approximated by a smooth function x + ξ log(1 + e x/ξ ).
With Adam, the step size and direction at each point are computed in an adaptive way, by taking into account the past history of the gradients to estimate their first and second moments.In the simulations presented thereafter, Adam was run with a learning rate lr = 0.001, betas = (0.9, 0.999), and no weight decay.We stopped the algorithm typically after 8000 iterations, as the convergence was estimated to be reached, which induced up to 150 seconds of computation time per shape with a simulation domain of size 100 × 100 × 100 pixels.
Generation parameters, shape textures, curvature diagrams.The generation of shape textures relies on the principle that a generation parameter vector (a, m 0 ) should consistently correspond to a single shape texture, across different white noise initializations of A 0 .We "measure" the texture of a shape defined by a surface S through its curvature diagram, which represents the distribution of the curvatures (κ 1 , κ 2 ) on S (see bottom row of Figure 2 for instance).More precisely, we are interested in the law of the random variable (κ 1 , κ 2 ) defined by The curvature diagram is an indicator of the local behavior of the surface.For a perfect sphere of radius R, it should be a unit Dirac mass δ R everywhere on a sphere.Likewise, a cylinder of base radius r should have its diagram reduced to δ ( 1 r ,0) on the horizontal half-line {y = 0, x > 0}; a plane would be represented as δ (0,0) ; finally, a Dirac mass such as δ (1,−1) should correspond to a sponge-like shape.Note that, per definition of the curvatures κ 1 ≥ κ 2 , the distribution is contained in the lower mid-plane {y ≤ x}.
To obtain a curvature diagram, we first extract the 2D mesh of the surface S = {u = 0} from the 3D volume u, by using the marching cubes algorithm [59].Then the diffuse curvatures (κ 1,u , κ 2,u ) (see (17) and (18)) are interpolated at the barycenter of each cell of the mesh.Their values are of importance proportional to the area of the cell, resulting in a weighted point cloud which we plot in the curvature diagram.In our simulations, the shapes are rarely perfectly spherical, cylindrical, or flat, so that the distribution is dispersed rather than concentrated into a single Dirac mass.To ease the visualization of the curvature diagrams, the identity diagonal {y = x} is enhanced as a red line and the values truncated between −100 and 100.
We compare curvature diagrams with each other using the Wasserstein (or Earth Mover's) distance.We approximate this quantity with the regularized Sinkhorn algorithm of the geomloss module [35], with the parameters p = 2,  blur = 1 and reach = 20.As detailed in [34], these correspond to the resolution of an unbalanced transport problem [96] using a ground cost function of C(x, y) = 1 2 x − y 2 , with a transport plan that is blurred at resolution of 1 and with a maximum transport distance of the order of 20.
Figure 2 shows that the curvature diagrams of the evolving level surface also converge towards a final diagram.In Figure 4, we check that for a single generation parameter value (a, m 0 ), five different initializations still give similar curvature diagrams.We found that the mean pairwise Wasserstein distance between them was only 10.7% of the mean pairwise distance measured between 1000 random shapes (see Experiment 4).Therefore, curvature diagrams and generation parameters seem to capture well the notion of shape texture.

Experiment 1: validation with the
Willmore flow and known minimizers of genus 0, 1, and 2 The algorithmic framework, which combines automatic differentiation and control of gradient flows by external optimizers, is validated in the fundamental special case of the L 2 Willmore flow (see Algorithm 2, with the parameters a = (1, 2, 1, 0, 0, 0), and with replicate boundary conditions 4 ).Let us recall from Section 2.3 that the Willmore phase-field energy writes F = 1 Ω ( ∆u − W (u) ) 2 dx and approximates σ S H 2 dA.We numerically check that the simulated Willmore flow converges towards known global minimizers of fixed genus 0, 1, and 2, by initializing the flow near them.The gradient descent is controlled by the L-BFGS optimizer, as it was empirically found to converge faster than Adam.L-BFGS approximates the BFGS algorithm, a quasi-Newton method that combines a line search to an estimation of the Hessian of the loss [36,49].
The Willmore minimizers of genus 0 are spheres of any radius [103], which achieve the minimal value5 For surfaces of genus 1, [62] proved that the minimal value 4 × 2π 2 is achieved by the Clifford torus (up to conformal transformations), defined by a special ratio 1 √ 2 between the radius of the generating circle and the distance to Bottom row: their corresponding curvature diagrams.The five shapes are visually similar at a mesoscopic scale, but dissimilar at the macroscopic scale.The diagrams and their pairwise Wasserstein distances show that the curvature statistics are not changing much when feeding the same generation parameters to the algorithm.Texture in shapes hence seems to be well captured by curvature diagrams, and to be consistent with the generation parameters at a mesoscopic scale.
the axis of revolution.However, the proof for genus ≥ 2 has still not been completed, although several conjectures have been proposed.It has been shown that the minimum Willmore energy among all (orientable closed) surfaces of genus g is less than 4 × 8π, and converges to this value as the genus g → ∞ [57].The Lawson surfaces have also been conjectured to be the minimizers for a given genus (up to conformal transformations) [56,42].
Algorithm 2 L 2 flow: follow the gradient flow of the phasefield energy F until convergence, with given initialization and replicate boundary conditions Initialization: a phase-field u0 Generation parameters: coefficients a = (a2,0, a1,1, a0,2, a1,0, a0,1, a0,0) Energy: phase-field energy F (see (20)) to approximate S ( |α|≤2 aα(κ1, κ2) α ) dA Other parameters: phase transition parameter > 0, internal parameters for L-BFGS (learning rate, history size, line search function, maximal number of iterations in a line search), number of iterations T , Gaussian kernel of size σ k Outputs: phase-field u and surface S = {u = 0} S ← {u = 0} Visualize the shape We compare the final value F σ to the minimal values mentioned above, where σ = 4 3 √ 2 is the constant introduced in Section 2. We correctly find that the flow is stationary on spheres (see Figure 5).The energy F σ deviates from 4 × 4π with only 0.3% of relative error.We also check that the value of the Cahn-Hilliard energy Ω ( 2 |∇u| 2 + W (u) ) dA, divided by σ, is close to the area of the sphere, with 0.5% of relative error.When departing from a holed cube, the flow converges to a Clifford torus with a characterizing ratio close to 1 √ 2 up to a relative error of 3%.We also find F σ 4 × 2π 2 up to a relative error of 0.4%.Finally, starting from two holed cubes glued together, the flow converges to a surface resembling a Lawson surface of genus 2, with minimal value F 4σ 22.30, which has a relative difference of 1.8% compared to the value found in [42].In Figures 5, 6, and 7, the parameters are = 0.04, a = (1, 2, 1, 0, 0, 0) and σ k = 1 pixel.The grid size is 200 × 200 × 200 in the first two simulations and 200 × 300 × 400 in the third one.L-BFGS was run with a learning rate lr = 1, a history size hs = 10, and maximum 20 iterations in a line search.

Experiment 2: a gallery of ten shape textures
In Figures 11 and 12, we show ten shape textures generated with Algorithm 1, and visualize the surfaces together with their curvature diagrams.The coefficients a and the mass m 0 = 1 |Ω| Ω u used for the simulations are specified in the tables.The mass is also expressed in percentage of volume   enclosed by the surface compared to the total volume of the domain, approximated by m0+1 2 since u ±1.
In practice, we find that the non-reduced polynomial expression 2 ) that are more interpretable than the coefficients a in the reduced form |α|≤2 a α (κ 1 , κ 2 ) α .Using this formulation (26), it is easier to find reasonable values for which the energy leads to various tubular textures, without resulting in badly-converged phase-fields with no zero level set, or shapes broken into small fragments.The shapes were hence generated either by choosing b manually, or by selecting a among random values from Experiment 4 that resulted in interesting shapes.
We found that curvatubes is able to generate very different shape textures (see Figures 11 and 12, in which they are indexed by letters).Some shapes, such as (d), (e), (g), and (h) are smooth and spatially homogeneous in terms of visual aspect.Other shapes, such as (b), (c), (f), and especially (i), seem to possess a multi-scale texture or "meta-texture" that makes them appear as spatially heterogeneous and anisotropic.The surface can be piecewise-smooth only, as in (c).Tubules are not necessarily smoothly turning cylinders, but can have some tortuosity as in (j).Finally, note that (h) combines flat regions and tubules, in a similar way to trabecular bone.

Experiment 3: bilinear interpolation between four shape textures
In Figure 8, we illustrate how continuously varying the generation parameters (a, m 0 ) impacts on the morphologies, by interpolating the parameters of four shape textures: layers, spheres, tubes, and sponges.The respective values can be found in Table 1.All simulations are run by starting from the same initialization A 0 .We can see that the morphology is smoothly changing throughout the figure: for instance, the transition between spheres and tubes is characterized by tubes terminated on one side by end-caps, while layers increase in proportion compared to tubules when approaching the top left corner.The curvature diagrams are displayed in Figure 9, and they are quite continuously evolving as we change the generation parameters.The four diagrams at the corners reflect well the typical curvature distribution of layers, spheres, tubes, and sponges, as expected (see Section 4.1).

Experiment 4: generation of 1000 shapes viewed in UMAP
Our final experiment is designed to explore the space of possible shape textures with curvatubes and visualize them in a 2D atlas (see Figure 10).The generation parameters were chosen randomly.We fixed a 2,0 = 1, and arbitrarily chose the other coefficients according to a uniform law in the following intervals: By doing this, the algorithm was pushed to its limits, as some values of coefficients chosen in this random way led to an ill-posed geometric minimization problem.Yet, even in these cases, the algorithm did not diverge to NaN values, but the function u simply did not converge to a phase-field with two distinct phases, or had no zero level set, or the zero level set was not smooth and was frag- (2, 2, 2, 0, 0, 0) (1, 0, 0, 1, 0, 1, 0) Table 1: Generation parameters of four reference shape textures in Experiment 3 (see Figures 8 and 9).The vector b parameterizes the polyomial expression (26).mented into pieces.To reject such situations, we gauged the viability of the parameters by computing the discrepancy Ω 2 |∇u| 2 − W (u) dx of the phase-field, normalized by the diffuse area, i.e., the Cahn-Hilliard energy, The discrepancy measures how much u deviates from a phase-field with tanh profile as in (19), and is an indicator of a good behavior in numerical experiments.After normalization by the diffuse area, the quantity obtained varies between 0 and 1, with 0 indicating a good quality in the numerical approximations.
Shapes were deemed viable if the normalized discrepancy was under the threshold 0.75 and if max(u) > 0.1 and min(u) < −0.1 to ensure that the zero level set was   8).The bilinear interpolation of generation parameters is also reflected into the curvature diagrams as a continuous evolution.Please note that the curvature distributions of the four reference shapes concentrate around (0, 0) for layers, the identity half-line {y = x, x > 0} for spheres, the horizontal half-axis {y = 0, x > 0} for tubes, and a half-line {y = −a x, x > 0} with a > 0 for sponges.As explained in Section 4.1, this is expected from curvature diagrams.
defined.If the random value assigned to (a, m 0 ) produced a non-viable shape, a new value was drawn uniformly until obtaining a viable shape.
We thus generated 1000 shapes meeting the criteria mentioned above, and computed the pairwise Wasserstein distance between their curvature diagrams, as described in the last paragraph of Section 4.1.To reduce the computation time, in each comparison we restricted the point cloud (25) to 10000 cells randomly taken from the whole mesh.The distance matrix was then given as input to UMAP [65], a manifold learning technique for dimension reduction, with the option metric = 'precomputed'.The 1000 shapes were embedded in a 2D atlas by considering local neighborhoods of n = 25 points, a minimal distance md = 0.05 between embedded points and a spread sp = 1.For reproducibility, the random seed was set to RS = 1.To enhance the visu-alization, the embedded points were labeled with Hdbscan [64], with a minimum cluster size mcs = 10 and a minimum number of samples ms = 10.They were colored according to their cluster number, and as black dots if Hdbscan classified them as noise.We picked some shapes from the point cloud and displayed their thumbnails, with their location specified by an arrow.Please note that the thumbnails do not exhaustively cover all the types of morphologies, but may give an idea of their disparity.
The atlas in Figure 10 shows that the shape textures are roughly distributed into three main families, spheres (top, left), layers (top), and tubules (central part) that occupy the largest region.In two marginal regions, we also identified outliers, such as highly packed tubules (bottom, right), and fragmented shapes (bottom, left).The latter suggest that the selection criteria mentioned above were not selective enough for discarding badly-converged shapes.The marginal regions concentrate most of the outliers, but we also noticed a few of them spread inside the main regions.
The transition between morphological subtypes is quite smooth when moving continuously in the atlas; however we did not examine in which way the generation parameters relate to the spatial embedding yet.The family of tubules has a large intra-variation, and features not only smooth sponges or long tubes, but also irregular, tortuous and anisotropic tubules that have a multi-scale texture, as mentioned in Experiments 2. This is a remarkable behavior of the generation model, since in regard of the minimization problem, all points have the same homogeneous properties in space.
Figure 10: Atlas of 1000 shapes visualized in UMAP.The generation parameters were chosen randomly and we applied some criteria to select mainly viable shapes.Three major families can be identified: spheres, layers, and tubules which constitute the principal type of shape textures.Let us highlight the large intra-variation in the family of tubules, that features not only smooth sponges and long tubes, but also irregular and anisotropic tubules with a multi-scale texture.This is a remarkable feature of the generation model, since all pixels of the simulation domain share identical properties with respect to the minimization problem.

Discussion
In this final section, we discuss the strengths and limitations of the algorithm, propose a few extensions, and present the important implications of a unifying framework on applied contexts.

Strengths and limitations
As seen in the simulations of Section 4, curvatubes leads to a wide range of membranous and tubular shapes, some of which have a multi-scale texture.The generation parameters and the curvature diagrams capture well the notion of shape texture.The algorithm is GPU-accelerated and takes advantage of automatic differentiation combined to external algorithms (Adam, L-BFGS) to descend gradient flows.Contrarily to refined numerical schemes, its aim is not to precisely solve the evolution equation, but rather to converge fast to a local minimizer with small energy.The mathematical computation of the gradient is not required either.The coefficients can be chosen in a flexible way, without letting the algorithm diverge numerically, even in mathematically ill-posed cases.The simulation results are reproducible, and seem to behave in accordance with the polynomial of curvatures in the energy especially under the form (26); although much work is still needed to understand mathematically how different polynomial energies are linked to different shape textures, and what are the values of coefficients that correspond to well-posed geometric problems.
Let us caution the reader that here, the model does not generate tree-like structures, for which junctions are hierarchically organized into parent and children nodes, and cycles are excluded.Thus, it cannot be applied to the respiratory system, and can only model vessels that branch and cycle a lot such as capillaries.This is because the model is devised primarily for reproducing shape texture, but not shape structure.Some extensions of the framework to include structured constraints are proposed in the next subsection.The notion of shape texture is inspired from visual texture in images, characterized by spatially repeated elements whose conformation, such as size, color, orientation, are subject to randomness [43,80,58].Texture is hence a statistically defined property, while structure may be understood as an orthogonal component.
Furthermore, in contrast to the Helfrich biomembrane model and the FCH model of [38], which truly model some physico-chemical energy derived from microscopic interactions, we do not assume any such physical ground to the general curvature functional F and the corresponding phase-field F that we propose.This framework is simply intended to provide a descriptive tool to analyze tubular textures, and may be used to quantify biological shapes in terms of geometry, even without any knowledge of the underlying microscopic interactions.

Extensions
We can include a constraint on the orientation of the normal vector n S to the surface S, by encouraging n S to be orthogonal to the direction associated to a vector θ of unit norm, as in This can be approximated by the phase-field energy where F is our phase-field expression constructed in Section 3.2.The effect on tubes is to align their median axis along θ, while inciting flat layers to be parallel to θ, as in Figure 13. 2 , m0 = −0.6, and µ = 0 or 1000 for the first and second shapes.The layers were generated using a = (1, 1, 1, 0, 0, 0), which corresponds to F = (H 2 + κ 2 1 + κ 2 2 )/2, m0 = −0.4,and µ = 0 or 800 for the third and fourth shapes.When µ = 0, the direction of θ is indicated by the arrow.The curvature diagrams show that after alignment, the curvatures are more densely clustered around the dirac masses associated to perfectly cylindrical or perfectly flat shapes, namely, δ (20,0) for tubes and δ (0,0) for layers.
Another way to give some structure to the shape is to use space-dependent generation parameters (a(x), m(x)), i.e., make them spatialized instead of constant, as in Figure 1.
In the current version of the algorithm, since u is periodic, the coefficients a(x) but also the mass m are required to be periodic.The change of variable of Section 4.1 becomes In Figure 1, we took four reference generation parameters p 1 , p 2 , p 3 , and p 4 (see Table 2), and linearly interpolated them along the horizontal axis, by taking into account the periodicity.We also repeated the first and last values (following the order p 1 , p 1 , p 2 , p 3 , p 4 , p 4 ), and cropped the shape by dropping the first and last cubes.The spatialized parameters hence coincide with p 1 , p 2 , p 3 , p 4 at the vertical midplanes of the four cubes delimited by the dashed lines.

Importance of a unifying theory, and future applications
Finally, we have identified several implications that a generation model unifying tubular and membranous shapes could have in other contexts.
• Design bio-inspired shape textures: the generation model could help in the design of bio-fabricated vascular networks for tissue regeneration [89], as well as scaffolds for bone tissue engineering [33] and cellular solids in architecture [72], both inspired by the structure of bone trabeculae.
• Model morphological states and trajectories: if generation parameters can be inferred from morphological states, a morphological transformation can be modeled as a trajectory in the lower-dimensional space of parameters, and then analyzed as a longitudinal trajectory [28].In particular, the biological transformations mentioned in the Introduction would be modeled in a continuous way.
• Provide regularization prior for tubular segmentation: the generation model could be included as a regularizing loss in variational segmentation methods of vascular structures [101,30,66], to select certain tubular morphologies against others.It could also be combined with 3D reconstruction from 2D slices methods [60,10,50].
• Build a synthetic database of textures: the generation algorithm could provide, at a low cost, a complete panel of synthetic textures on which to test and train vascular shape analysis methods [79,47,84], including topological analysis methods [74,44,12], segmentation algorithms, or microvascular blood flow simulations [81,3].It could also provide a database to research in shape and texture perception [102]; the way we perceive shapes is intimately linked to the way we want to quantify them.
We of course did not cover all these applications here, but intend to focus on two of the points aforementioned as future work.
The first one is to model morphological states or trajectories of biological tissues as static values or trajectories of generation parameters (a(t), m(t)), which supposes that parameters can be inferred from shapes.This can be done naively, by visual inspection and trial-and-error; or, by first producing an atlas of reference shapes densely sampling a region with the desired morphologies, similarly to Experiment 4. Using the curvature diagram of the query shape u 0 , the shapes closest to it in terms of the Wasserstein distance are found.We can then initialize (a, m 0 ) at these values, and minimize the loss with respect to the parameters (a, m 0 ) instead of the phasefield u 0 , by using nearly the same algorithmic framework as Algorithm 1.
The second related application is to include the curvature energy as a regularizing loss in order to segment vascular structures.The energy would then select certain tubular morphologies over others.This could be used for instance to reconstruct 3D tubular structures captured in several 2D images at different depths of a biological sample (as done in optical sectioning), provided that there are not filaments too thin compared to the diffusion width .The method is most effective if there is some knowledge of the shape textures that need to be segmented, so that the parameters can be appropriately tuned by inference, as previously explained.We split Ω into three regions, Ω (1) = {0 < |d| < b }, (2) = {b < |d| < c } and Ω (3) = {c < |d|} (on which u = ±1).On Ω (1) and Ω (2) , ∇u = 0 and M u = |u |B u , whereas on Ω (3) , ∇u = 0 and M u = 0. Note that the region Ω (1) ∪ Ω (2) = {d ≤ c } decreases and concentrates around the surface S as → 0, while the complementary region Ω (3) grows.Also, on Ω Let φ ∈ C 0 c (Ω) be a continuous function with compact support.We need to show that Ω |∇u | 2 φ dx converges to σ S φ dH 2 as is sent to zero.First, we work only on where Jac(ψ) = ( ∂ψ ∂x , ∂ψ ∂y ) ∈ R 3×2 , and this integral converges to V φ dH 2 when s → 0. By compactness, we can consider a finite number of such neighborhoods V and conclude that the limit (22) holds.Now, we prove the convergence of the energies.It can be shown that on Ω (1) , we have such that f |∇u | ≤ 2 (1+ 2 ) 3 + 4 1+ 2 remains bounded on Ω (2) .

9 :
Update A with one step of Ȧ = − ∂G ∂A Adam 10:

Figure 2 :
Figure 2: Evolution of the zero level set of u during the H −1 flow, pictured at iterations 10, 300, 1000, 4000, 8000, with the generation parameters a = (1, 2.4, 9, 30, 170, −195) and m0 = −0.66.Top row: the shape textures.Bottom row: their corresponding curvature diagrams.In the early iterations, the zero level set is not smoothly defined and only few values of u are above zero, but their average remains m0.

Figure 3 :
Figure 3: Corresponding loss curves (see Figure 2).From left to right: the loss F (u), the maximal value of | ∂G ∂A | and the average 1 |Ω| Ω | ∂G ∂A |, as u evolves along the iterations.

Figure 4 :
Figure 4: Same generation parameters, different initializations.Using the same generation parameters a = (1, −0.35, 1.02, −40, 100, 1600) and m0 = −0.69,we start from five random initializations for A0.Top row: the shape textures.Bottom row: their corresponding curvature diagrams.The five shapes are visually similar at a mesoscopic scale, but dissimilar at the macroscopic scale.The diagrams and their pairwise Wasserstein distances show that the curvature statistics are not changing much when feeding the same generation parameters to the algorithm.Texture in shapes hence seems to be well captured by curvature diagrams, and to be consistent with the generation parameters at a mesoscopic scale.

Figure 5 :
Figure 5: Willmore minimizer of genus 0. The flow is stationary on the sphere.The diffuse Willmore energy and the diffuse area are close to their mathematical values.

Figure 6 :
Figure 6: Willmore flow towards minimizer of genus 1.The flow converges to a Clifford torus.The ratio of the torus and the diffuse Willmore energy are close to their mathematical values.

Figure 7 :
Figure 7: Willmore flow towards minimizer of genus 2.The flow converges to a Lawson-like surface, as expected from conjectures.

Figure 8 :
Figure 8: Bilinear interpolation between layers, spheres, tubes, and sponges.Using the same initialization A0 and interpolating between four reference generation parameters gives rise to a continuum of shape textures.

2 ( 2 Figure 12 :
Figure12: Five shape textures (f, g, h, i, j).We indicate the cases where the polynomial in the curvature energy has the form(26).