# A variational method for curve extraction with curvature-dependent energies

Majid Arthaud, Antonin Chambolle, Vincent Duval

<figure><img src="image_1.png" /></figure>

To cite this version:

Majid Arthaud, Antonin Chambolle, Vincent Duval. A variational method for curve extraction with curvature-dependent energies. 2025. hal-05386995

HAL Id: hal-05386995

https://hal.science/hal-05386995v1

Preprint submitted on 28 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A variational method for curve extraction with
curvature-dependent energies

Majid Arthaud<sup>1,2,4</sup>, Antonin Chambolle<sup>3,2[0000-0002-9465-4659]</sup>, and Vincent Duval<sup>2,3[0000-0002-7709-256X]</sup>

<sup>1</sup> ENPC, 6 Av. Blaise Pascal, 77420 Champs-sur-Marne, France
majid.arthaud@eleves.enpc.fr

<sup>2</sup> INRIA Mokaplan, INRIA Paris, Paris-Dauphine, CNRS, France
vincent.duval@inria.fr

<sup>3</sup> CEREMADE, CNRS and Université Paris-Dauphine, PSL, Paris, France
antonin.chambolle@ceremade.dauphine.fr

<sup>4</sup> M.A. is now at University of Michigan, Ann Arbor, United States of America

**Abstract.** We introduce a variational approach for extracting curves between a list of possible endpoints, based on the discretization of an energy and Smirnov's decomposition theorem for vector fields. It is used to design a bi-level minimization approach to automatically extract curves and 1D structures from an image, which is mostly unsupervised. We extend then the method to curvature-dependent energies, using a now classical lifting of the curves in the space of positions and orientations equipped with an appropriate sub-Riemannian or Finslerian metric.

**Keywords:** Charges · Numerical analysis · Bi-level problems · Curve extraction · Geodesic curves · Active contours.

## 1 Introduction

This paper is built upon the conference contribution [2], which is extended towards more general line energies, 3D examples and curvature-dependent energies.

We consider the problem of detecting (open) curves in images, by a variant of active contours type models. Although our final goal will be to simultaneously extract multiple curves and find their endpoints, for simplicity we first discuss the extraction of one curve given its endpoints. Active contours [24,9,14] usually rely on the minimization of a potential along a curve, with functionals of the form

$$ \int_{\Gamma} g(x, \tau) d\mathcal{H}^1 \qquad (1) $$

which are minimal when the curve $\Gamma$ passes through the lower values of the potential $g$ defined from the observed image (here, $\mathcal{H}^1$ is the length measure, given by the 1-dimensional Hausdorff measure). We consider a general setting, where the potential is also allowed to depend on the tangent direction $\tau$ to the curve, yet, to start with, the reader may discard this point. In general in
these models, the curve is assumed to be closed, but one could also specify the endpoints, in order to find $\Gamma$ as a minimal length curve (anisotropic geodesic) joining these two points [14,17]. See for instance [32] (in particular Section 3) for a general overview of these techniques.

A standard approach is then to solve an Eikonal equation (by fast-marching or fast-sweeping) to evaluate the distance to one of the endpoints (or several [13]), and compute then the geodesic to some other endpoint(s). An advantage is that the computation is fast, efficient, and can easily be adapted to many interesting frameworks (3D, surfaces, nontrivial liftings such as in [28,29] for the tracking of the width of retina vessels, or as in [19] for their local orientation — where here it is crucial to allow for a spatial and orientation dependent weight $g$). Alternatively, the problem of computing geodesics may be reformulated as a primal-dual problem as in [20].

We consider here a (seemingly) different type of approach, based on the representation of paths as measure fields whose divergence is concentrated on their endpoints, as proposed in [26,27]. Consider for instance the problem:

$$ \min_z \left\{ \int g(x, \frac{z}{|z|}) d|z|(x) : -\operatorname{div} z = \delta_B - \delta_A \right\} \quad (2) $$

where $z$ is a measure vector field of total variation measure $|z|$, and $A$ and $B$ are two given endpoints, while $z/|z|$ is the direction in the polar decomposition of the vectorial measure $z$. A typical example of such a measure field is given by the integration along curves, that is, measures of the form $\tau H^1 \llcorner \Gamma$, where $\tau$ is tangent to $\Gamma$ and $\Gamma$ is a curve. We show in Theorem 2 below (thanks to the celebrated Smirnov theorem on extreme points of such vector fields [36]) that the minimizers of (2) are superpositions of measures of this form, where each curve is geodesic between $A$ and $B$ for the cost $\int_{\Gamma} gdH^1$. One can then recover the curve(s) by integrating along the vector field $z$.

Another path to this equivalence is through convex duality, which shows that (2) is essentially equivalent to the standard approach described in [32], based on the Eikonal equation. Yet a slight difference is that in (2), one can consider more general constraints on $\operatorname{div} z$, such as, for instance, consisting of atomic measures with vanishing total mass, as we propose in Section 4. This allows to compute simultaneously many geodesics between a possibly large family of points. In that case, the equivalence with the Eikonal equation approach is broken and the setting we discuss becomes of interest.

After detailing our approach and showing how to derive an algorithm for extracting curves from an image by automatically detecting their endpoints, which was essentially the contents of our contribution in the SSVM 2025 conference proceedings [2], augmented by an extension in dimension 3, we show how it is extended to energies which additionally encompass a convex function of the curvature, following the framework in [11] where the approach is developed for a modified total variation regularizer. A possible further extension could then be to incorporate this setting in variational formulations for more general inverse problems (such as a “deconvolution” of thin structures) as in [26,27].
The paper is organized as follows. Before describing our numerical approach in Section 3 (which corresponds to finding appropriate discretizations of (2)), we explain more precisely in the next section the connection between (1) and (2). Then, in Section 4, we introduce a bi-level method for automatically placing the endpoints. In Section 5, we explain how the proposed framework may be adapted to the “roto-translational” representation of oriented curves, in order to build weights g which penalize the curvature. Some proofs, and the precise description of the optimization algorithm for solving the main convex problem, are postponed to the Appendix.

## 2 Charges and curves

To highlight the connection between the active contour functional (1) and our model (2), we consider the space of “normal charges” $\mathcal{V}$ which was recently brought to the attention of the image processing community in [26]. The idea consists in embedding finite curves in a linear space of vector valued measures, similar to the spaces of (normal) 1-currents [21] or 1-flat chains [37].

### 2.1 The space of normal charges

We work in a compact arcwise connected set $Q \subseteq \mathbb{R}^d$ (in practice $d$ will be 2 or 3, and $Q$ will be either a rectangle/cube or a cube periodic in one direction). The space of normal charges in $Q$ is the set of vector finite Radon measures supported in $Q$ and whose divergence is a finite Radon measure, that is,

$$ \mathcal{V} = \{z \in \mathcal{M}(\mathbb{R}^d)^d : \mathrm{spt}(z) \subseteq Q, \mathrm{div}(z) \in \mathcal{M}(\mathbb{R}^d)\}, \quad (3) $$

where $\mathcal{M}(\mathbb{R}^d)$ is the space of finite (signed) Radon measures over $\mathbb{R}^d$. It is a normed space when equipped with the norm

$$ \forall z \in \mathcal{V}, \quad \|z\|_{\mathcal{V}} = \|z\|_{\mathrm{TV}} + \|\mathrm{div}(z)\|_{\mathrm{TV}}, \quad (4) $$

where $\| \cdot \|_{\mathrm{TV}}$ denotes the total variation of (vector or scalar) Radon measures.

As normal charges are Radon measures, it is possible to define convex functionals of charges. We follow here the presentation of [1, Sec. 2.6]. Given a lower semi-continuous function $g: \mathbb{R}^d \times \mathbb{R}^d \to [0, +\infty]$, positively 1-homogeneous and convex in the second variable, we may define

$$ \forall z \in \mathcal{V}, \quad G(z) = \int g \left( x, \frac{dz}{d|z|}(x) \right) d|z|(x). \quad (5) $$

Then, the functional $G: \mathcal{V} \to [0, +\infty]$ is convex positively homogeneous [1, Prop. 2.37] and Reshetnyak's lower semi-continuity theorem [1, Th. 2.38] ensures that it is sequentially weak-* lower semi-continuous. We usually denote $G(z)$ by $\int g(x, z)$.

Furthermore, it is possible to define the measure $g(x, z)$ (see [7,16]) and to provide integral representations for the convex conjugate of $G$ (as in [7]).
## 2.2 Charges induced by curves

An important example of a normal charge is induced by oriented curves of finite length, that is Lipschitz functions $\gamma: [a, b] \to Q$ for some $a, b \in \mathbb{R}$.

The charge $z_\gamma$ is then the vectorial measure defined by

$$ \forall \varphi \in C_c(\mathbb{R}^d; \mathbb{R}^d), \quad \langle z_\gamma, \varphi \rangle = \int_a^b \varphi(\gamma(t)) \cdot \gamma'(t) \, dt. \qquad (6) $$

One checks that $\mathrm{div}(z_\gamma) = \delta_{\gamma(a)} - \delta_{\gamma(b)}$, so that $\mathrm{div}(z_\gamma) = 0$ if and only if the curve is closed.

It is possible to give a geometric interpretation to (6). Indeed, by the area formula and its consequence, the generalized change of variable [1, eq. 2.47], the following formula holds for every bounded Borel function $g: [a, b] \to \mathbb{R}$,

$$ \int_{\mathbb{R}^d} \left( \sum_{t \in [a, b] \cap \gamma^{(-1)}(y)} g(t) \right) d\mathcal{H}^1(y) = \int_a^b g(t) |\gamma'(t)| \, dt. $$

Setting $g(t) = 1$ if $\gamma'(t) = 0$ and $0$ otherwise, we see that for $\mathcal{H}^1$-a.e. $y \in \mathbb{R}^d$, $\gamma'(t) \neq 0$ for every $t \in \gamma^{(-1)}(y)$. Hence, given a test vector-field $\varphi \in C_c(\mathbb{R}^d; \mathbb{R}^d)$, we may set $g(t) := \varphi(\gamma(t)) \cdot \frac{\gamma'(t)}{|\gamma'(t)|}$ if $|\gamma'(t)| \neq 0$, and $0$ otherwise, and we get

$$ \int_{\mathbb{R}^d} \varphi(y) \cdot \left( \sum_{t \in [a, b] \cap \gamma^{(-1)}(y)} \frac{\gamma'(t)}{|\gamma'(t)|} \right) d\mathcal{H}^1(y) = \int_a^b \varphi(\gamma(t)) \cdot \gamma'(t) \, dt \\ = \langle z_\gamma, \varphi \rangle. $$

If the curve is simple, that is, if $\gamma$ is one-to-one on $[a, b)$, for $\mathcal{H}^1$-a.e. $y \in \mathrm{Im}\gamma$, there is a unique $t \in \gamma^{(-1)}(y)$ (and $\gamma$ is differentiable at $t$). We may thus define the tangent vector as $\tau(y) = \frac{\gamma'(t)}{|\gamma'(t)|}$, and the above equality states that the charge is equal to

$$ z_{\gamma} = \tau \mathcal{H}^1 \llcorner (\mathrm{Im}\gamma) \quad \text{and} \quad |z_{\gamma}| = \mathcal{H}^1 \llcorner (\mathrm{Im}\gamma). \qquad (7) $$

Note that, as a consequence of the above discussion, (6) is invariant by reparametrization, and, therefore, it is always possible assume that $|\gamma'(t)| \le 1$ for a.e. $t$ and $\gamma$ is defined on $[0, L]$ for some $L > 0$.

## 2.3 Smirnov's decomposition theorem

The landmark paper [36] shows that any normal charge may be described as a superposition of charges induced by curves and their generalization. In order to state the main result we are interested in [36, Theorem C], we briefly describe its setting.
Given some normal charge $z \in \mathcal{V}$, we say that $z$ decomposes into $p \in \mathcal{V}$ and $q \in \mathcal{V}$ if

$$z = p + q \quad \text{and} \quad |z| = |p| + |q|. \qquad (8)$$

The above equalities are understood in the sense of measures, and the expression $|r|$ for $r \in \mathcal{V}$ denotes its variation measure, defined as $|r|(E) = \sup \sum_i |r(E_i)|$ for all Borel set $E \subseteq \mathbb{R}^d$, where the supremum is over all finite Borel subdivisions of $E$. In particular $|r|(\mathbb{R}^d) = \|r\|_{TV}$.

Furthermore, we say that $z$ completely decomposes into $p$ and $q$ if (8) holds and

$$|\operatorname{div} z| = |\operatorname{div} p| + |\operatorname{div} q|. \qquad (9)$$

Those definitions extend straightforwardly to finite (or integral) sums of charges.

In order to decompose some charge into curves of finite length, we endow the set of curves with the structure of a compact metric space. In view of Section 2.2, up to a reparametrization, a curve and the corresponding charge may be determined by some function $\gamma: \mathbb{R} \to \mathbb{R}^d$ which is 1-Lipschitz. Smirnov defines $\text{Lip}_{1,d}^\wedge$ as the collection of all such functions together with the curve “at infinity” $f_\infty: t \mapsto \infty$. He endows it with a distance which metrizes uniform convergence on compact sets, and which makes it a compact space. It is then possible to define a Borel measure on that space as in the following theorem.

**Theorem 1 ([36, Thm. C]).** Let $z \in \mathcal{V}$. Then there exist two normal charges $p, q \in \mathcal{V}$ such that $z$ completely decomposes into $p$ and $q$, $\operatorname{div} p = 0$, and $q$ completely decomposes into simple oriented curves of finite length. In other words, there exists some nonegative Borel measure $\sigma$ on $\text{Lip}_{1,d}^\wedge$ such that

$$q = \int z_\gamma \, d\sigma(\gamma) \qquad (10)$$

$$|q| = \int |z_\gamma| \, d\sigma(\gamma), \qquad (11)$$

$$|\operatorname{div} q| = \int |\operatorname{div} z_\gamma| \, d\sigma(\gamma). \qquad (12)$$

The “conic combinations” in (10-12) are in the weak-* sense, e.g.

$$\forall \varphi \in C_c(\mathbb{R}^d), \quad \langle q, \varphi \rangle = \int \langle z_\gamma, \varphi \rangle \, d\sigma(\gamma). \qquad (13)$$

In turn, $p$ can also be completely decomposed as well, but the decomposition may involve elementary solenoids, that is charges induced by generalized curves (see [36, Def. and Thm. B]). A generalized curve may consist in an infinitely winding curve of infinite length, but since we show below that these objects do not appear in the solutions of our variational problem, we do not describe them more precisely. See [36] for more detail.
An alternative interpretation of Smirnov's results is that the extreme points of the unit ball of $||\cdot||_V$ are the measures supported on simple oriented rectifiable curves, of norm one. It has led Laville et al. [26,27] to use (4) as a regularizer in inverse problems, as they prove that some solutions of their variational problem are superpositions of a finite number of measures induced by simple curves.

## 2.4 A minimization problem with prescribed divergence

Here, we focus on solving problems of the form

$$ \min_z \left\{ \int_Q g(x, \frac{z}{|z|}) d|z|(x) : -\operatorname{div} z = \mu \right\} \quad (14) $$

where $\mu$ is a given measure with support in $Q$ and $\mu(Q) = \int_Q d\mu = 0$. This kind of problem is known in the field of branched transportation as Beckmann's problem (see for instance [34, Sec. 4.2] or [30]). We make the following assumptions on the weight $g$.

Assumptions:

(A0) There exists $D > 0$ such that for every $(x, y) \in Q^2$, there exists a Lipschitz curve $\gamma: [0, 1] \to Q$ with $\gamma(0) = x$, $\gamma(1) = y$, and

$$ \int_0^1 g(\gamma(t), \gamma'(t)) dt \le D. \quad (15) $$

(A1) The function $g: \mathbb{R}^d \times \mathbb{R}^d \to [0, +\infty]$ is lower semi-continuous, and it is convex positively 1-homogeneous in the second variable.

(A2) There exists some constant $c > 0$ such that $g(x, t) \ge c$ for all $(x, t) \in Q \times \mathbb{S}^{d-1}(\mathbb{R})$, where $\mathbb{S}^{d-1}(\mathbb{R})$ denotes the unit sphere of $\mathbb{R}^d$.

The following result is a continuous version of known and standard results in graph theory and network flows (see also the discussion below (22)).

**Theorem 2.** Let $\mu \in \mathcal{M}(\mathbb{R}^d)$ with support in $Q$ and $\mu(Q) = 0$, and assume that (A0), (A1) and (A2) hold.

Then, there exists a minimizer $z$ to (14), and for any such minimizer there exists a Borel positive measure $\sigma$ defined on $\text{Lip}^{\wedge}_{1,d}$, the set of oriented 1-Lipschitz curves such that

$$ z = \int z_{\gamma} d\sigma(\gamma), \quad (16) $$

$$ |z| = \int |z_{\gamma}| d\sigma(\gamma), \quad (17) $$

$$ |\operatorname{div}(z)| = \int |\operatorname{div}(z_{\gamma})| d\sigma(\gamma), \quad (18) $$

$$ \int g(x, z) = \int \left( \int g(x, z_{\gamma}) \right) d\sigma(\gamma). \quad (19) $$

Moreover, $\sigma$-a.e. $\gamma$ is open, simple, and
(i) $\gamma$ is a weighted geodesic curve in $Q$, i.e., minimizes $\int g(\gamma(t), \gamma'(t))dt$ for fixed beginning point $b(\gamma)$ and endpoint $e(\gamma)$,

(ii) $b(\gamma) \in \text{spt } \mu^{-}$ and $e(\gamma) \in \text{spt } \mu^{+}$ where $\mu = \mu^{+} - \mu^{-}$ is the Hahn-Jordan decomposition of $\mu = -\text{div } z$.

The proof of Theorem 2 is given in Appendix A.

In the particular case where $\mu$ is a sum of Dirac masses $\sum_{i} \pm \delta_{x_i}$, the solution is carried by geodesic curves with endpoints on the $x_i$'s (oriented from the negative to the positive Diracs).

In [27], Laville et al. exploit a representer theorem which bounds the number of curves to solve their problem in an “off-the-grid” greedy approach relying on the Frank-Wolfe algorithm. The implementation is non-trivial, and requires a non-convex step. Since our setting possibly yields many more curves, we rather consider a convex approach based on a discretization. The difficulty arises in properly discretizing singular vector fields.

*Remark 1.* Contrary to what is written in [2], in general, the measure $\sigma$ is not a probability measure.

We now introduce a discretized version of these normal charge, in order to use this representation in practical image analysis to represent 1D curves in images. This raises a few difficulties since in particular, a discrete version of Theorem 2 does not exist in general, except for very elementary anisotropic curve energies.

# 3 Discrete curves

In this section, we describe different discretization strategies and discuss their performance. To make the description easier to follow, we first consider the bi-dimensional case (the generalization to 3D will be straightforward) and weights which do not depend on the orientation. Orientation-dependent weights (Section 5) will require some additional caution.

We consider 2D images of $N \times M$ pixels, and we introduce $\mathcal{N} = \{(i, j) : 1 \le i \le N, 1 \le j \le M\}$ the set of nodes and $\mathcal{E} = \{(i + \frac{1}{2}, j) : 1 \le i \le N - 1, 1 \le j \le M\} \cup \{(i, j + \frac{1}{2}) : 1 \le i \le N, 1 \le j \le M - 1\}$ the set of edges between the neighbouring nodes (where $(i + \frac{1}{2}, j)$ denotes the edge between the nodes $(i, j)$ and $(i + 1, j)$, etc).

## 3.1 Flux on a graph

Playing the same role as the space of normal charges above, the space of discrete vector fields is defined as

$$V_d = \mathbb{R}^{\mathcal{E}} \simeq (\mathbb{R})^{(N-1) \times M} \times (\mathbb{R})^{N \times (M-1)} \quad (20)$$

For a vector field $v \in V_d$, the component $v_{i+\frac{1}{2},j}$ may be seen both as a field between $(i, j)$ and $(i+1, j)$ (on the edge), and a flux through the facet separating the two pixels.
An important case of a vector field is the “gradient” of an image $u \in \mathbb{R}^N \simeq \mathbb{R}^{N \times M}$, defined as the following finite difference operator:

$$ (Du)_{i+\frac{1}{2},j} = u_{i+1,j} - u_{i,j}, \quad (Du)_{i,j+\frac{1}{2}} = u_{i,j+1} - u_{i,j}. \qquad (21) $$

In this convention, the “horizontal” derivative $(Du)_{i+\frac{1}{2},j}$ is defined on the edge between $(i,j)$ and $(i+1,j)$, for $i = 1, \dots, N-1$, $j = 1, \dots, M$ and similarly the “vertical” derivative is defined on the edge between $(i,j)$ and $(i,j+1)$, for $i = 1, \dots, N$ and $j = 1, \dots, M-1$. Then, a discrete divergence operator is naturally defined, for $v \in \mathcal{V}_d$, as $D^*v \in \mathbb{R}^N$, given by $\langle D^*v, u \rangle_{\mathbb{R}^N} = \langle v, Du \rangle_{\mathcal{V}_d}$ for all $u$ (with the scalar products given by the canonical Euclidean products on the respective spaces). With our convention in the definition of $D$ (we only compute a difference when the two points are in the discrete domain), $D^*$ is a divergence with vanishing flux condition on the boundary.

The advantage of this setting is its similarity with the continuous one (Sec. 2). A straightforward adaptation of (14) is

$$ \min_{z \in \mathcal{V}_d, D^*z = \mu} \sum_{e \in \mathcal{E}} g_e |z_e| \qquad (22) $$

where $(g_e)_{e \in \mathcal{E}}$ is a collection of positive weights. It benefits from a discrete counterpart to Theorem 2, which comes from well known results in graph theory (see, e.g. [5, Prop. 3.10 and Rem. 3.12]).

Nevertheless, as we illustrate below (Fig. 1, center), this $\ell^1$-norm yields blurry results, due to the non-uniqueness of the corresponding geodesics. In addition, it measures the “length” of the loops or curves in a very anisotropic way, only through their horizontal and vertical projections. In the active contour community, that phenomenon is well known, and people prefer to use discretizations of the isotropic Eikonal equation together with fast marching approaches over Dijkstra’s algorithm on the graph of the image [14] in order to compute geodesic.

<figure><img src="image_3.png" /><figcaption>Fig. 1. Left: 200 × 200 “noisy comma” shape (crop from [27, Fig. 1]), middle: result with weighted $\ell^1$ norm, right: weighted $\ell^2$ norm.</figcaption></figure>
One has therefore to consider other norms, consistent (as $M, N \to \infty$) with the Euclidean norm in (14), even if in doing so we lose at the discrete level the decomposition theorem of a field as a superimposition of curves.

## 3.2 Isotropic discretizations

We introduce a family of discrete approximations of the convex curve reconstruction problem (14), with different anisotropy properties. Given a weight function $g \in \mathbb{R}^N$ with $0 \le g_{i,j} \le 1$ (derived from the gray level values of an image), and the scalar field of curve endpoints $\mu \in \mathbb{R}^N$, the general form of our problem is:

$$ \min_{z \in V_d, D^*z = \mu} \sum_{i,j} g_{i,j} \|(Az)_{i,j}\|, \qquad (23) $$

where $A$ is an operator which sends $z$ to a $\mathbb{R}^2$-valued field in $\mathbb{R}^N$ and $\|\cdot\|$ a norm in $\mathbb{R}^2$. For instance, the weighted counterpart of the $\ell^1$-norm $\|z\|_1$ (as in (22)) corresponds to choosing

$$ (Az)_{i,j} = \begin{pmatrix} z_{i+\frac{1}{2},j} \\ z_{i,j+\frac{1}{2}} \end{pmatrix}, \quad \text{and } \forall x \in \mathbb{R}^2, \quad \|x\| = |x_1| + |x_2|. \qquad (24) $$

By solving (23) we extract, through the discrete vector field $z$, dark curves (where $g$ is close to zero) on a light background (where $g$ is close to one). An experimental result with (24) is shown in Fig. 1, middle: the weighted $\ell^1$ result is quite blurry, due to the fact that several curves have almost the same energy, and that the optimization outputs a convex combination of these. In addition, horizontal and vertical lines are strongly favored by the energy.

A variant yielding sharper and more isotropic results consists in replacing the $\ell^1$-norm in (24) with the Euclidean norm $\|x\| = \sqrt{x_1^2 + x_2^2}$. This corresponds to a consistent discretization of (14) with forward differences. The result is shown in Fig. 1, right. The arbitrary choice of forward-forward differences for the gradient still induces some anisotropy at small scale. This is clearly seen when comparing the bottom left and the bottom right of the comma, the first being more blurry than the second (see also Fig. 2, top right).

A possibility to correct this anisotropic behaviour and obtain a sharper result is through averaging. As already mentioned, in the spirit of discrete calculus [23], the discrete vector field is living on the edges between the cells of the pixels, while the scalar fields of the grayscale levels or of the curve extremities should be seen as living on these cells. To re-center the norms of the $z_{i,j}$ on the cells, we replace in (24) the operator $A$ by a true averaging operator:

$$ (Az)_{i,j} = \frac{1}{2} \begin{pmatrix} z_{i-\frac{1}{2},j} + z_{i+\frac{1}{2},j} \\ z_{i,j-\frac{1}{2}} + z_{i,j+\frac{1}{2}} \end{pmatrix}, $$

and using the $\ell^2$-norm again in (23). This yields the discrete problem:

$$ \min_{z \in V_d, D^*z = \mu} \sum_{(i,j) \in \mathcal{N}} g_{i,j} \sqrt{ \left( \frac{z_{i-\frac{1}{2},j} + z_{i+\frac{1}{2},j}}{2} \right)^2 + \left( \frac{z_{i,j-\frac{1}{2}} + z_{i,j+\frac{1}{2}}}{2} \right)^2 } \qquad (25) $$
Notice that here, we have defined the coefficients outside of the image with zero padding (that is, $z_{\frac{1}{2},j} = z_{N+\frac{1}{2},j} = z_{i,\frac{1}{2}} = z_{i,M+\frac{1}{2}} = 0$), which is consistent with the vanishing flux condition on the boundary. The results are very isotropic and

<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Left: result with the weighted and averaged $\ell^2$-norm on the same “noisy comma” image. Right: zoom of the same result (bottom) and zoom without averaging (top).</figcaption></figure>

sharper than those obtained through the weighted $\ell^2$ norm, as seen in Fig. 2, left. A zoom on the bottom part compares these two versions (Fig. 2, right). In

<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Left: an angiogram (source: Wikipedia), right: geodesic computed by our method (compare for instance with [17, Fig. 4(d)]).</figcaption></figure>

Fig. 3, we illustrate the use of this method on an angiogram. The result is a geodesic similar to the results obtained in the literature on minimal paths [17].

### 3.3 Optimization

All of the experimental results exposed in this document have been implemented with a primal-dual algorithm [10], or “PDHG”, accelerated by a reprojection on the constraint $D^*z = \mu$ computed with FFTW3 [22], which solves the Lagrangian saddle point problem:

$$ \min_{z \in \mathcal{V}_d} \max_{p \in (\mathbb{R}^N)^2} \langle Az, p \rangle + \chi_{\{z: D^*z=\mu\}} - \chi_{\{p: \|p_{i,j}\|_* \le g_{i,j} \ \forall i,j\}} \quad (26) $$
Here, $\|\cdot\|_*$ is the dual norm of the norm $\|\cdot\|$ in (23) (defined by $\|q\|_* = \sup_{\|p\|\le 1} q\cdot p$), and $\chi$ denotes a characteristic function in the classical sense of convex analysis (0 if the condition is satisfied, $+\infty$ else). The PDHG algorithm described in Appendix B is elementary to implement and only requires the matrix-vector products with $A$ and its adjoint $A^*$, and the knowledge of the projections onto the sets whose characteristic function appear in (26). The main advantage of this approach is its versatility: it is straightforward to adapt to any linear operator $A$, or to change the penalization of $(Az)_{i,j}$ by suitably modifying the constraint set for $p$ (here, $\|p_{i,j}\|_* \le g_{i,j}$; more complex convex, one-homogeneous penalizations will be considered in Section 5). See the full description in Appendix B.

The idea of using an averaging operator comes from the context of discretizing the total variation, and was introduced in [15] and further developed in [12]. An interesting topic for research would be to analyze and reproduce more closely their findings, or develop variants in the context of curve extraction.

<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Result on a 3D potential: the left images correspond to the 70×70×30 potential, and the right images are the retrieved curves. Now the result is obtained after 2000 primal-dual steps, with fixed endpoints, for a wall time of about 6 min.</figcaption></figure>

The extension to a 3D volume of this method is straightforward, with the same formulas, the same averaging and the same algorithm. We show an example (on synthetic data and as before with given endpoints) in Figure 4.

# 4 An iterative discrete curve reconstruction algorithm

## 4.1 Finding the endpoints

The discrete convex curve reconstruction problem of the previous section assumes that the endpoints (represented by discrete Dirac masses) of the curves to
reconstruct were given. Now, we propose an iterative algorithm which attempts to retrieve discrete curves in images with no prior knowledge of their positions. The only parameter are the maximal number of endpoints and a threshold $g_{max}$, chosen so that the curves should roughly describe the set $\{g \le g_{max}\}$. Formally, given $g$ and $g_{max}$, and a maximum number of source/sink pairs $n_0$, we aim at solving the non-convex bi-level optimization problem (here $Q$ is a rectangle representing the domain of the image):

$$ \min_{\substack{1 \le n \le n_0 \\ (s_i, t_i)_{i=1}^n \in Q^{2n}}} \left\{ \int_Q \left( g\left(x, \frac{dz}{d|z|}(s)\right) - g_{max} \right) d|z|(x) : z \text{ solves (14) for } \mu = \sum_{i=1}^n \delta_{t_i} - \sum_{i=1}^n \delta_{s_i} \right\}. \quad (27) $$

The idea is to try to find the endpoints of the curves present in the image, by maximizing the length of the minimal curves between these endpoints inside the sublevel set $\{g \le g_{max}\}$.

To implement this minimization, we need to describe the gradient with respect to a source/sink point $s_i$ or $t_i$ of the energy in (27). Since the optimal measure is expected to be carried by geodesic curves joining each $s_i$ to some point $t_j$, the largest change for the energy corresponds to moving $s_i$, $t_j$ along the geodesic (that is, in direction $\pm z(s_i)$ or $\pm z(t_i)$), towards the direction which shortens the curve if $g(x) > g_{max}$, and the direction which lengthens the curve if $g(x) < g_{max}$. At equilibrium, we also check whether moving the sources/sinks in the direction of $\nabla g$ improves the criterion. Also, we sometimes need to decrease the number of sources/sinks.

<figure><img src="image_7.png" /><figcaption>**Fig. 5.** For the discretization of the continuous curve drawn in black and the Dirac mass at its endpoint, we represent the four different cases of the algorithm.</figcaption></figure>

We now describe the discrete implementation in two dimensions. In what follows, for $x \in \mathcal{N}$ we denote by $\delta^x \in \mathbb{R}^N$ a discrete Dirac mass at $x$, given by $\delta^x_{i,j} = 1$ if $x = (i, j)$ and $\delta^x_{i,j} = 0$ else. The algorithm is the following:

1. Choose a set of initial points $S \subset \mathcal{N}$ and a set of final points $T \subset \mathcal{N}$ of same cardinality $n_0$, and initialize $\mu = \sum_{s \in S} \delta^s - \sum_{t \in T} \delta^t \in \mathbb{R}^N$ as a sum of discrete Dirac masses with total sum zero;

2. Solve (23) (e.g. with the primal-dual algorithm [10], see Appendix B);

3. For each Dirac mass at some $x = (i, j) \in \mathcal{N}$: estimate the discrete curve orientation near $x$, by the average of the discrete vector field $z$ in the $3 \times 3$ square around $x$. We distinguish four cases:
(a) The **shortening** case: if $g > g_{max}$ at the current pixel of the Dirac mass, move the Dirac mass to the pixel in the 3 × 3 square in the curve shortening direction;

(b) The **lengthening in the curve direction** case: if the condition of the previous case is not satisfied, and $g \le g_{max}$ on one of the pixels opposite or at 45° of the opposite of the curve shortening direction, and this pixel has not yet been visited by this Dirac mass, move the Dirac mass to the pixel verifying this condition with the lowest potential $g$;

(c) The **lengthening in an orthogonal direction** case: if the conditions of the two previous cases are not satisfied, and $g \le g_{max}$ on one of the pixels in one of the two orthogonal directions to the shortening direction, and this pixel has not yet been visited by this Dirac mass, move the Dirac mass to the pixel with this condition and the lowest potential $g$;

(d) The **shifting** case: if the conditions of the previous three cases are not satisfied, and $g \le g_{max}$ on one of the pixels in one of the two directions at 45° of the curve shortening direction, and this pixel has not yet been visited by this Dirac mass, move the Dirac mass to the pixel verifying this condition with the lowest potential $g$;

4. Reiterate steps 2 and 3 until all of the Dirac masses have converged.

The directions in the 3 × 3 square around a Dirac mass are represented in Fig. 5. The Dirac masses converge in practice, as they cannot go from the “lengthening and shifting” stage to the “shortening” stage, and in the “lengthening and shifting” stage, each pixel can only be visited once. In some instances (see next Section), we may have to merge some endpoints/curves and start again the algorithm to improve the results, ending in a number of endpoints which is always less than the initial choice. It remains unclear how to extend properly these rules in 3D and to cope with the fact that 3D curves are more likely to avoid the zones of interest, which makes trickier the initialization to reproduce our 2D results. This is a topic for future study and experiments.

<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Result of the iterative algorithm after post-processing, on a 200 × 200 image of 42 bird chromosomes. Obtained after 100 iterations of the iterative algorithm, each consisting in 60 primal-dual steps, for a total of 1 min. 05 s. wall time on an average laptop. A post-processing, which consists in 5000 primal-dual steps on the last set of Dirac masses, adds another 1 min. 12 s. wall time.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Fig. 7.** Result of the iterative algorithm on the 200 × 200 “comma” shape. The iteration number *k* for the algorithm is indicated. Top: μ initialized as a pair of Dirac masses of opposite intensities randomly chosen in the sublevel set {$g \leq g_{max}$}. Middle: μ initialized as 15 pairs of Dirac masses of opposite intensities randomly chosen in the sublevel set {$g \leq g_{max}$}. Bottom: same as “middle”, but the “comma” shape is blurred with a Gaussian filter of standard deviation 2. Each iteration had 80 primal-dual steps. The middle experiment was the longest, taking 5 min. 07 s. wall time after 390 iterations on an average laptop. A post-processing of 5000 primal-dual steps on the last set of Dirac masses added 1 min. 03 s. wall time.</figcaption></figure>

## 4.2 Numerical results

We show a few results obtained with our implementation<fnref n="5" />. The computation of the optimal *z* given *μ* relies on a primal-dual method [10], see Appendix B for details. In practice, getting a good estimate for *z* which solves (23) at each iteration requires few steps, since we initialize the variables of the algorithm as the variables from the last step of the previous iterate (and the measure *μ* does not change much from one step to the next). Moreover, we use simple projection steps and an adapted choice of the step sizes. The implementation is developed

<footnote n="5">⁵ available at https://github.com/majidarth/Discrete-curve-reconstruction.</footnote>
in Python and does not use parallelization (outside of the default Numpy multi-threading). We indicate in a few captions the wall time of the execution on an average laptop.

The success of the iterative bi-level algorithm relies upon the quality of the initialization of $\mu$: in the case of an image with many curves to reconstruct, in order to make sure to capture all of them, one may need to initialize $\mu$ with a very big set of pairs of discrete masses of intensity +1 and -1, randomly over the whole image, such that each two Dirac masses of a pair are close (e.g. in each other's 4 × 4 neighborhood). For the result presented in Fig. 6, $\mu$ was initialized as 2000 such pairs of Dirac masses, while there are only 42 chromosomes in the original image$^6$.

When $\mu$ is initialized as a very large number of Dirac pairs, some additional processing is required. First, at each iteration, whenever two Dirac masses of opposite coefficients are in each other's 3 × 3 neighborhood, it is considered that they should merge and they are removed from $\mu$. Furthermore, after the Dirac masses have converged, it is important to post-process them. In practice, many pairs of Dirac masses may converge close to the endpoints of the same acquisition of a curve in the original image: they need to be merged. In order to do this, the curves need to be deduced from the discrete vector field $z \in V_d$ (i.e., $z$ is decomposed into a superposition of open curves whose endpoints are the Dirac masses). Then, every time two pairs of opposite Dirac masses correspond to curves that superpose and are roughly of the same orientation, these curves are essentially merged, keeping only the two Dirac masses (of opposite signs) which correspond to the endpoints of the global curve, and discarding the other two. With the resulting set of Dirac masses, (23) is then computed one last time, and the resulting $z$ is the final result. In addition to Fig. 6, we show three results on the “comma” image on Fig. 7: one with a random initialization with two endpoints, and two with a random initialization with 15 pairs of Dirac masses in the region of low potential $g$. The last one, obtained after slightly blurring the potential, allows to recover nicely and completely the main curve in the image.

# 5 Curvature penalization

## 5.1 Lifting of curvature-dependent energies

We now detail how the method is extended to incorporate curvature penalization, using the celebrated representation in the group of “roto-translations” introduced and developed in [35,6,33,18,3,4,19,25] (and many other works by the same groups). Yet as before, instead of relying on oriented edge detection [4] or solving an anisotropic eikonal equation to minimize the length of curves [19], we build a convex energy of vector fields with divergence constraints. We follow the framework in [11], which addresses the implementation of a “total roto-translational

<footnote>$^6$Initially from [31, Fig. 2], with a rescaling, a change of contrast and to which noise was added</footnote>
variation”, that is an energy which enforces a curvature penalization of the gradient of a function. Up to a 90° rotation, in 2D, this is equivalent to penalizing the mass of a zero-divergence field, and we can re-use much of the numerical approach of [11] for our implementation. Adaption to higher dimension would be possible (and not very difficult), yet computationally hardly tractable.

We want to penalize a planar curve, parameterized by a Lipschitz map $\gamma: [0, 1] \to R$ where $R$ is the image domain (a closed, arcwise connected set, in general a rectangle), by a curvature dependent energy of the form:

$$ \int_{0}^{1} g(\gamma(t))f(\kappa_{\gamma}(t))|\gamma'(t)|dt, \qquad (28) $$

for $f$ an even, convex function, with $f \ge 1$, and $g$ as before a (positive, continuous) varying weight. Here, $\kappa_{\gamma}(t)$ is the (absolute) curvature, that is, the length of the component of $\gamma''/|\gamma'|$ orthogonal to $\gamma'$ (it is given by $|\gamma''|$ if $\gamma$ is parameterized so that $|\gamma'| = 1$ a.e. in $[0, 1]$).

In practice, we lift the image domain $R$ by adding a variable corresponding to the orientation of the curves, setting $Q = R \times \mathbb{S}^1$ where $\mathbb{S}^1 := \mathbb{R}/(2\pi\mathbb{Z})$ is the periodic circle of length $2\pi$. By convention, we will denote $(x, \theta)$ a point in $R \times \mathbb{S}^1$ and $\theta = (\cos\theta, \sin\theta)$ the point in the circle defined by the angle $\theta$.

In [11], it is shown that — following the framework in [35,6,33,18,3,4,19] — one can represent (28) by lifting the curve $\gamma: [0, 1] \to R$ as a curve $\Gamma(t) := (\gamma(t), \theta(t)) : [0, 1] \to Q$ where $\theta$ is the angle of $\gamma(t)'$ with the direction $(1, 0)$ (so that $\theta = \gamma'/|\gamma'|$ and the curvature is $\kappa_{\gamma} = \theta'/|\gamma'|$), and weighting $\Gamma$ with an appropriate “length”, as we describe now. We introduce $f^{\infty}(t) = \lim_{a \to +\infty} f(at)/a$, the recession function of $f$ at infinity. Then, we define $\bar{h}$ as the perspective function of $f$: (cf [11, eq. (3)]):

$$ \bar{h}(s, t) = \begin{cases} sf(t/s) & \text{if } s > 0, \\ f^{\infty}(t) & \text{if } s = 0, \\ +\infty & \text{else.} \end{cases} \qquad (29) $$

Equivalently, $\bar{h}$ is the support function of $\{(a, b) \in \mathbb{R}^2 : f^*(a) + b \le 0\}$, that is:

$$ \bar{h}(s, t) = \sup \{as + bt : f^*(a) + b \le 0\}. \qquad (30) $$

Then, for $\theta$ an angle (with $\theta$ the associated unit vector), and $v = (v^x, v^\theta) \in \mathbb{R}^2 \times \mathbb{R}$, we let:

$$ h(\theta, v) = \begin{cases} \bar{h}(\lambda, v^\theta) & \text{if } v^x = \lambda\theta \text{ for some } \lambda \in \mathbb{R}, \\ +\infty & \text{else.} \end{cases} $$

In particular,

$$ h(\theta, v) = \sup \{a \cdot v^x + bv^\theta : (a, b) \in \mathbb{R}^2 \times \mathbb{R}, f^*(a \cdot \theta) + b \le 0\}. \qquad (31) $$

One can check that with such a definition,

$$ \int_{0}^{1} g(\Gamma^{x}(t))h(\Gamma^{\theta}(t), \Gamma'(t))dt = \int_{0}^{1} g(\gamma(t))f(\kappa_{\gamma}(t))|\gamma'(t)|dt. \qquad (32) $$
The curvature naturally pops out because the vertical component $\Gamma^{\theta'}$ of $\Gamma'$ is precisely the derivative of the direction of the curve. We refer to [11] for the computational details.

Now, to adapt our setting to this framework, we consider normal charges $z = (z^x, z^\theta)$ supported in $Q$, and their horizontal projection (or marginalization) $\tilde{z}$ onto $R$ defined by $\tilde{z}(B) := z^x(B \times \mathbb{S}^1)$ for any borel set $B \subset R$. Observe that if

$$-\operatorname{div} z = \delta_b - \delta_a$$

(3D divergence) in $Q$, then, denoting $a = (a^x, a^\theta) \in R \times \mathbb{S}$, etc, we find that:

$$-\operatorname{div} \tilde{z} = \delta_{b^x} - \delta_{a^x},$$

in particular a curve in $Q$ between the endpoints $a, b$ is projected onto a curve in $R$ with endpoints $a^x, b^x$. As before, given $\mu^+, \mu^-$ two measures in $Q$ with $\mu^+(Q) = \mu^-(Q)$ (in practice, we use sums of Dirac masses), we can solve the problem:

$$\min \left\{ \int_Q g(x)h(\theta, z) : -\operatorname{div} z = \mu^+ - \mu^- \right\}. \qquad (33)$$

Thanks to Theorem 2 (provided $g$ and $R$ are such that (A0-2) hold), there exists a measure $\sigma$ which decomposes $z$ into curves $\Gamma$ of finite length, with properties (16-19). In particular, $\sigma$-a.e. curve $\Gamma$ is minimal for the energy (32), so that its horizontal projection is minimal for (28). If $\mu^\pm$ are given by $n$ Dirac masses, $n \ge 1$, we find that in general (up to a possible non-uniqueness of the geodesic curves), the horizontal projection $\tilde{z}$ consists therefore in $n$ (possibly overlapping) curves minimizing the curvature dependent energy (28) among all curves joining the points in the support of $\mu^-$ to the points in the support of $\mu^+$.

## 5.2 Numerical implementation

In this new context, the discretization of (33) has the form

$$\min_{z \in \mathcal{V}_d^{RT}, D^* z = \mu} \sum_{i,j,k} g_{i,j} h(\theta_k, (Az)_{i,j,k})$$

where $\mathcal{V}_d^{RT} \simeq (\mathbb{R})^{(N-1) \times M \times K} \times (\mathbb{R})^{N \times (M-1) \times K} \times (\mathbb{R})^{N \times N \times K}$ is periodic in the last component and represents the space of charges in $Q = R \times \mathbb{S}^1$, discretized on a $N \times M \times K$ grid, $N, M, K \ge 1$, and $\theta_k = 2k\pi/K$. The operator $D$ is still a discrete differentiation operator (with adjoint $D^*z$), with the third component given by:

$$Du_{i,j,k+\frac{1}{2}} = \begin{cases} u_{i,j,k+1} - u_{i,j,k} & k = 1, \dots, K-1 \\ u_{i,j,1} - u_{i,j,K} & k = 0. \end{cases}$$

The implementation relies again on optimizing the 3D version of (26):

$$\min_{z: D^*z = \mu} \max_{p \in C_{g,h}} \langle p, Az \rangle \qquad (34)$$
and in order to implement our algorithm, one now needs to be able to project on the constraint set $C_{g,h}$ whose support function is the energy. Thanks to (31), this is given by:

$$C_{g,h} = \{ p_{i,j,k} = g_{i,j}(\boldsymbol{q}_{i,j,k}^x, q_{i,j,k}^\theta) \in (\mathbb{R}^3)^{N \times M \times K} : f^*(\boldsymbol{q}_{i,j,k}^x \cdot \boldsymbol{\theta}_k) + q_{i,j,k}^\theta \le 0 \forall i, j, k \}.$$

Here as before, we have defined the horizontal vector $\boldsymbol{\theta}_k = (\cos \theta_k, \sin \theta_k)$. Since each component $(i, j, k)$ is independent, the projection onto $C_{g,h}$ is built upon the projection onto the set whose $\bar{h}$ is the support function, cf (30):

$$C_{\bar{h}} := \{(a, b) \in \mathbb{R}^2 : f^*(a) + b \le 0\}.$$

Then, to obtain the projection $q = \Pi_{C_{g,h}}(p)$ of a dual variable $p$ onto $C_{g,h}$, we compute for each coordinate $(i, j, k)$:

$$\begin{align}
(a, b) &= g_{i,j} \Pi_{C_{\bar{h}}} \left( \frac{1}{g_{i,j}} (\boldsymbol{p}_{i,j,k}^x \cdot \boldsymbol{\theta}_k, \boldsymbol{p}_{i,j,k}^\theta) \right) \nonumber \\
\boldsymbol{q}_{i,j,k}^x &= \boldsymbol{p}_{i,j,k}^x + (a - \boldsymbol{p}_{i,j,k}^x \cdot \boldsymbol{\theta}_k) \boldsymbol{\theta}_k \tag{35} \\
\boldsymbol{q}_{i,j,k}^\theta &= b. \nonumber
\end{align}$$

As in [11, Sec. 4.2], we considered the following choices for $f$ (where each time $\alpha$ is a positive parameter):

- $f_1(t) = 1 + \alpha|t|$ ("TAC", total absolute curvature);

- $f_2(t) = \sqrt{1 + \alpha^2 t^2}$ ("TRL", total roto-translational length—TRV in [11]);

- $f_3(t) = 1 + \alpha^2|t|^2$ ("EL", Elastica energy—TSC in [11]).

The first and second allow for brutal change of directions in the curves (which are penalized by the turning angle), while the third enforces smoothness of the curves. The corresponding sets $C_{\bar{h}}$ are given by:

$$C_{\bar{h}_1} = \{(a, b) \in \mathbb{R}^2 : a \le 1, |b| \le \alpha\},$$

$$C_{\bar{h}_2} = \{(a, b) \in \mathbb{R}^2 : \max(0, a)^2 + (b/\alpha)^2 \le 1\},$$

$$C_{\bar{h}_3} = \{(a, b) \in \mathbb{R}^2 : a + b^2/(2\alpha)^2 \le 1\}.$$

The projections are implemented as in [11], actually we re-used the c++ programs developed by Thomas Pock for this purpose. While the projection onto $C_{\bar{h}_1}$ is straightforward, the projections onto $C_{\bar{h}_2}$ (for $\alpha \ne 1$) and $C_{\bar{h}_3}$ require to solve a non-linear problem and rely on a few iterations of a Newton method. We refer to [11, Sec. 4.2] for details.

## 5.3 Rules for moving the endpoints

We also adapted the method of Section 4 to find the endpoints of a family of curves. In that case, the endpoints are pairs $(x, \theta)$ with $x = (i, j)$ and $\theta =$
<figure><img src="image_10.png" /><figcaption>**Fig. 8.** Minimization of the EL energy with $\alpha = 1$, on a 100×100 version of the “comma” image (note that the image is half the size here compared to previous examples). Left: projection on the image domain, center and right: result in the lifted image domain. This was obtained after 1500 primal-dual steps, with fixed endpoints, for a total of 4 min. 50s. wall time on an average laptop.</figcaption></figure>

$2\pi k/K$ for some integers $(i, j, k)$. The update of $(i, j)$ follows the same rules as in Section 4.1. The update of $k$ is simpler: at the free endpoints, we would like the curvature to vanish and the curves to satisfy $\Gamma^{\theta'}(0) = \Gamma^{\theta'}(1) = 0$. Hence if $z_{i,j,k \pm \frac{1}{2}}$ is positive (up to some threshold), we move up or down (depending whether the point is in $\mu^+$—leaving or $\mu^-$—entering, which is given by the sign of the Dirac mass) the Dirac mass at $(i, j, k)$, while if negative we move it in the other direction.

<figure><img src="image_11.png" /><figcaption>**Fig. 9.** Left: TAC, right: EL, on the 100 × 100 “comma” image, with parameter $\alpha = 5$, 1500 primal-dual steps with fixed endpoints. Observe that the TAC allows for jumps in the direction, which are prevented by the EL.</figcaption></figure>

## 5.4 Numerical experiments

We show a few numerical results for curvature-dependent energies. These results require more processing time than curvature-independent implementations, since the underlying computations are in the 3D lifted variables. In every example in this section, the angle $\theta$ is discretized into 30 possible values (with the notation of Section 5.2, $K = 30$). We consider the energies described earlier, namely the “TAC”, the “TRL” and the “EL”. We also show some results in the 3D lifted
domain, where the length and width of the volume represent the domain of the image, and the height represents the angle $\theta$, as detailed in Section 5.1. We start by looking at a simple example, recovering the previous “comma” shape with the EL energy, see Figure 8. This energy yields a very smooth curve. Then in Figure 9 we compare the TAC and the EL, with a stronger curve penalization ($\alpha = 5$): the curve obtained with EL is clearly smoother than the one obtained with TAC, which allows for jumps in the direction.

We then compare in Figure 10 the three different curvature-dependent energies, on a simple image with two crossing curves (with $\alpha = 1$ and fixed endpoints). The smaller image on the bottom right is the result obtained with the curvature-independent energy (25) (after 1800 primal-dual steps). Observe that in the latter case, the crossing is obviously not recovered, and the method finds a combination of many geodesic curves with almost same length which lie in the wide the low-potential area. On the other hand, the curvature penalized results are much more stable. The resulting curves are smoother for the TRL energy

<figure><img src="image_12.png" /><figcaption>**Fig. 10.** Top left: 80 × 50 pixel input image. Top right: TAC. Bottom left: TRL, center: EL. Small image: without curvature penalization. We used $\alpha = 1$, and 1800 primal-dual steps (with fixed endpoints). Wall time of roughly 4 min. 30 s. for every experiment, on an average laptop. The crossing is better reconstructed with the EL, which totally prevents part of the flow to make an abrupt turn towards a wrong direction.</figcaption></figure>

than for the TAC energy, and smoothest for the Elastica EL energy. Also, the two curves are not perfectly identified with TRL and TAC, since curves turning abruptly still have bounded energy—while they have infinite EL energy. This is better seen in the lifted image domain, see Figure 11.

We then present a result where the endpoints are not fixed, using the rules of Section 4 extended to the case of curvature-dependent energies (cf. Sec. 5.3), to evolve the endpoints; see Figure 12. Now, the curves can cross each other (unlike
<figure><img src="image_13.png" /><figcaption>**Fig. 11.** From left to right: results of Figure 10 represented in the lifted image domain for the TAC, the TRL and the EL.</figcaption></figure>

in Figure 6) so that curvature-dependent energies are more appropriate. In this experiment the iterative algorithm recovers all of the curves in the image.

<figure><img src="image_14.png" /><figcaption>**Fig. 12.** Final result of the iterative algorithm, on a 100 × 100 image, using the TRL energy with parameter $\alpha = 1$, after 10 iterations of the algorithm of Section 4, each iteration corresponding to 500 primal-dual steps, for a total of about 5 min. wall time on an average laptop. The Dirac masses were initialized as 1000 pairs of opposite intensities randomly chosen in the sublevel set $\{g \le g_{max}\}$.</figcaption></figure>

# 6 Conclusion

In this paper, we have introduced a variational method for curve extraction which differs from the usual Eikonal equation approach [32], as we represent paths as measure vector fields. It allows to simultaneously compute geodesics between a large numbers of endpoints. We also proposed a bi-level method to automatically fit the endpoints, yielding a complete, mostly unsupervised automatic curve extraction algorithm. It was extended to consider curvature penalization, allowing to reconstruct smooth or crossing curves. Further work could focus on more refined averaging operators to improve the sharpness and isotropy of the results, building upon recent findings in the context of discrete total variations.
Acknowledgements

A. Chambolle acknowledges the support of the “France 2030” funding ANR-23-PEIA-0004 (“PDE-AI”). Part of this work was done while M. Arthaud was an intern in the Mokaplan team, with the support of INRIA Paris. We thank Thomas Pock for allowing to re-use the code developed for [11].

References

1. Ambrosio, L., Fusco, N., Pallara, D.: Functions of Bounded Variation and Free Discontinuity Problems. Oxford Mathematical Monographs, Oxford University Press, Oxford, New York (Mar 2000)

2. Arthaud, M., Chambolle, A., Duval, V.: A variational method for curve extraction. In: Scale Space and Variational Methods in Computer Vision: 10th International Conference, SSVM 2025, Dartington, UK, May 18–22, 2025, Proceedings, Part II. p. 82–94. Springer-Verlag, Berlin, Heidelberg (2025). https://doi.org/10.1007/978-3-031-92369-2_7

3. Bekkers, E.J., Duits, R., Mashtakov, A., Sanguinetti, G.R.: A PDE approach to data-driven sub-Riemannian geodesics in SE(2). SIAM J. Imaging Sci. 8(4), 2740–2770 (2015), https://doi.org/10.1137/15M1018460

4. Bekkers, E., Duits, R., Berendschot, T., ter Haar Romeny, B.: A Multi-Orientation Analysis Approach to Retinal Vessel Tracking. Journal of Mathematical Imaging and Vision 49(3), 583–610 (Jul 2014). https://doi.org/10.1007/s10851-013-0488-6

5. Bonnans, F., Gaubert, S.: Recherche opérationnelle. Aspects mathématiques et applications. Ellipse (Mar 2016)

6. Boscain, U., Duits, R., Rossi, F., Sachkov, Y.: Curve cuspless reconstruction via sub-Riemannian geometry. ESAIM Control Optim. Calc. Var. 20(3), 748–770 (2014), https://doi.org/10.1051/cocv/2013082

7. Bouchitté, G., Valadier, M.: Integral representation of convex functionals on a space of measures. Journal of Functional Analysis 80(2), 398–420 (Oct 1988). https://doi.org/10.1016/0022-1236(88)90009-2, https://www.sciencedirect.com/science/article/pii/0022123688900092

8. Brown, L.D., Purves, R.: Measurable selections of extrema. The annals of statistics pp. 902–912 (1973)

9. Caselles, V., Catté, F., Coll, T., Dibos, F.: A geometric model for active contours in image processing. Numer. Math. 66(1), 1–31 (1993). https://doi.org/10.1007/BF01385685

10. Chambolle, A., Pock, T.: A first-order primal-dual algorithm for convex problems with applications to imaging. J. Math. Imaging Vision 40(1), 120–145 (2011). https://doi.org/10.1007/s10851-010-0251-1

11. Chambolle, A., Pock, T.: Total roto-translational variation. Numer. Math. 142(3), 611–666 (2019). https://doi.org/10.1007/s00211-019-01026-w, https://doi.org/10.1007/s00211-019-01026-w

12. Chambolle, A., Pock, T.: Learning consistent discretizations of the total variation. SIAM J. Imaging Sci. 14(2), 778–813 (2021). https://doi.org/10.1137/20M1377199, https://doi.org/10.1137/20M1377199

13. Cohen, L.D.: Multiple contour finding and perceptual grouping using minimal paths. J. Math. Imaging Vision 14(3), 225–236 (2001). https://doi.org/10.1023/A:1011281928379
14. Cohen, L.D., Kimmel, R.: Global minimum for active contour models: A minimal path approach. International Journal of Computer Vision **24**, 57–78 (1996)

15. Condat, L.: Discrete total variation: new definition and minimization. SIAM J. Imaging Sci. **10**(3), 1258–1290 (2017). https://doi.org/10.1137/16M1075247

16. Demengel, F., Temam, R.: Convex Functions of a Measure and Applications. Indiana University Mathematics Journal **33**(5), 673–709 (1984), https://www.jstor.org/stable/45010267, publisher: Indiana University Mathematics Department

17. Deschamps, T., Cohen, L.D.: Fast extraction of minimal paths in 3D images and applications to virtual endoscopy. Medical image analysis **5** **4**, 281–99 (2001)

18. Duits, R., Boscain, U., Rossi, F., Sachkov, Y.: Association fields via cuspless sub-Riemannian geodesics in SE(2). J. Math. Imaging Vision **49**(2), 384–417 (2014), https://doi.org/10.1007/s10851-013-0475-y

19. Duits, R., Meesters, S.P.L., Mirebeau, J.M., Portegies, J.M.: Optimal Paths for Variants of the 2D and 3D Reeds-Shepp Car with Applications in Image Analysis. Journal of Mathematical Imaging and Vision **60**(6), 816–848 (2018). https://doi.org/10.1007/s10851-018-0795-z

20. Ennaji, H., Quéau, Y., Elmoataz, A.: A primal-dual algorithm for computing Finsler distances and applications. Calcolo **61**(3), 53 (Aug 2024). https://doi.org/10.1007/s10092-024-00596-y, https://hal.science/hal-03620343

21. Federer, H.: Geometric measure theory, Die Grundlehren der mathematischen Wissenschaften, vol. Band 153. Springer-Verlag New York, Inc., New York (1969)

22. Frigo, M., Johnson, S.: The design and implementation of fftw3. Proceedings of the IEEE **93**(2), 216–231 (2005). https://doi.org/10.1109/JPROC.2004.840301

23. Grady, L., Polimeni, J.: Discrete Calculus: Applied Analysis on Graphs for Computational Science. Springer London (2010)

24. Kass, M., Witkin, A., Terzopoulos, D.: Snakes: Active contour models. International Journal of Computer Vision **1**(4), 321–331 (1988). https://doi.org/10.1007/BF00133570

25. Laville, B., Bertrand, T.: Dynamic off-the-grid untangling of curves with Reeds-Shepp metric. (Jun 2025), https://hal.science/hal-05124672, working paper or preprint

26. Laville, B., Blanc-Féraud, L., Aubert, G.: Off-the-grid curve reconstruction through divergence regularization: an extreme point result. SIAM J. Imaging Sci. **16**(2), 867–885 (2023). https://doi.org/10.1137/22M1494373

27. Laville, B., Blanc-Féraud, L., Aubert, G.: A $\Gamma$-convergence result and an off-the-grid charge algorithm for curve reconstruction in inverse problems. J. Math. Imaging Vision **66**(4), 572–583 (2024). https://doi.org/10.1007/s10851-024-01190-1

28. Li, H., Yezzi, A.: Vessels as 4d curves: Global minimal 4d paths to extract 3d tubular surfaces. In: 2006 Conference on Computer Vision and Pattern Recognition Workshop (CVPRW’06). pp. 82–82 (2006). https://doi.org/10.1109/CVPRW.2006.210

29. Li, H., Yezzi, A., Cohen, L.: 3D multi-branch tubular surface and centerline extraction with 4D iterative key points. In: Proceedings of the 12th International Conference on Medical Image Computing and Computer-Assisted Intervention: Part II. p. 1042–1050. MICCAI ’09, Springer-Verlag, Berlin, Heidelberg (2009). https://doi.org/10.1007/978-3-642-04271-3_126

30. Lohmann, J., Schmitzer, B., Wirth, B.: Formulation of branched transport as geometry optimization. Journal de Mathématiques Pures et Appliquées **163**, 739–779 (2022). https://doi.org/10.1016/j.matpur.2022.05.021
31. Malinovskaya, L., Zadesenets, K., Karamysheva, T., Akberdina, E., Kizilova, E., Romanenko, M., Shnaider, E., Scherbakova, M., Korobitsyn, I., Rubtsov, N., Borodin, P., Torgasheva, A.: Germline-restricted chromosome (GRC) in the sand martin and the pale martin (hirundinidae, aves): synapsis, recombination and copy number variation. Scientific reports 10(1) (Jan 2020). https://doi.org/10.1038/ s41598-020-58032-4

32. Peyré, G., Péchaud, M., Keriven, R., Cohen, L.D.: Geodesic methods in computer vision and graphics. Foundations and Trends in Computer Graphics and Vision 5(3-4), 197-397 (2010). https://doi.org/10.1561/0600000029

33. Prandi, D., Boscain, U., Gauthier, J.P.: Image processing in the semidiscrete group of rototranslations. In: Geometric science of information, Lecture Notes in Comput. Sci., vol. 9389, pp. 627-634. Springer, Cham (2015). https://doi.org/10.1007/ 978-3-319-25040-3_67

34. Santambrogio, F.: Optimal Transport for Applied Mathematicians: Calculus of Variations, PDEs, and Modeling. Progress in Nonlinear Differential Equations and Their Applications, Springer International Publishing (2015)

35. Sarti, A., Citti, G.: Subjective surfaces and Riemannian mean curvature flow of graphs. Acta Math. Univ. Comenian. (N.S.) 70(1), 85-103 (2000)

36. Smirnov, S.K.: Decomposition of solenoidal vector charges into elementary solenoids, and the structure of normal one-dimensional flows. Algebra i Analiz 5(4), 206-238 (1993)

37. Whitney, H.: Geometric integration theory. Princeton University Press, Princeton, NJ (1957)

A Proof of Theorem 2

We first provide the proof of Theorem 2, before stating some Lemmas on which it relies.

*Proof (of Theorem 2).* We begin by proving the existence of a solution $z$. By Lemma 1 below, Problem (14) is feasible, hence we may consider a minimizing sequence $(z_n)_{n \in \mathbb{N}}$. There exists a constant $C > 0$ such that for all $n \in \mathbb{N}$,

$$
C \geq \int g(x, z_n) \geq c |z_n| (\mathbb{R}^d),
$$

where we have used Assumption (A2) in the second inequality. Thus $(z_n)_{n \in \mathbb{N}}$ is
bounded in $M(\mathbb{R}^d; \mathbb{R}^d)$ and we may extract a subsequence $(z_n)_{n \in I}$ (with $I \subseteq \mathbb{N}$
infinite) which converges in the weak-$^*$ topology to some $z \in M(\mathbb{R}^d; \mathbb{R}^d)$. The
constraint $-\text{div} \tilde{z} = \mu$ is weak-$^*$ closed, hence $-\text{div} z = \mu$. Moreover the lower-
semi-continuity of $G: \tilde{z} \mapsto \int g(x, \tilde{z})$ for the weak-$^*$ convergence (see Section 2.1)
implies that

$$
G(z) \leq \liminf_{\substack{n \to +\infty, \\ n \in I}} G(z_n) = \inf(14).
$$

As a result, $z$ is a solution to (14).

Now, we prove the existence of a complete decomposition of for z. Smirnov's
theorem (Theorem 1) implies that z completely decomposes into two parts, p
and $q$, with $\operatorname{div} p = 0$ and $q$ is decomposed on curves of finite length in the support of $q$. By Lemma 2 below, this implies that

$$ \int g(x, z) = \int g(x, p) + \int g(x, q) \geq \int g(x, q), $$

and the minimality of $z$ for (14) implies that $p = 0$ and $z = q$. Hence, the same theorem ensures that $z$ can be completely decomposed as in (16, 17, 18) for some Borel measure $\sigma$. Moreover, $\sigma$-a.e. curve lies in $\operatorname{spt} z \subseteq Q$.

Our next step is to prove (19). This follows from [11, Lemma 5.2] (for the case where $g$ is continuous, which may be extended to $\operatorname{l.s.c.} g$ by a monotone convergence argument in [11, Cor. 5.4]). We propose below an alternative proof which relies on Lemma 2 and on Fatou's Lemma. Let us write

$$ \Gamma = \{ \gamma \in \operatorname{Lip}_{1,d}^{\wedge} : \forall t \in \mathbb{R}, \gamma(t) \in Q \}, $$

so that $\sigma(\operatorname{Lip}_{1,d}^{\wedge} \setminus \Gamma) = 0$. For every Borel partition $(\Gamma_1, \Gamma_2)$ of $\Gamma$, the decomposition formula (17) together with Lemma 2 yield

$$ \begin{aligned} \int g(x, z) &= \int g \left( x, \int_{\Gamma} z_{\gamma} \, d\sigma(\gamma) \right) \\ &= \int g \left( x, \int_{\Gamma_1} z_{\gamma} \, d\sigma(\gamma) \right) + \int g \left( x, \int_{\Gamma_2} z_{\gamma} \, d\sigma(\gamma) \right). \end{aligned} $$

By induction, if $(\Gamma_i)_{i \in I}$, with $I \subseteq \mathbb{N}$, is a Borel partition of $\Gamma$, we obtain similarly

$$ G \left( \int_{\Gamma} z_{\gamma} \, d\sigma(\gamma) \right) = \sum_{i \in I} G \left( \int_{\Gamma_i} z_{\gamma} \, d\sigma(\gamma) \right) \quad (36) $$

with the notation $G(z) = \int g(x, z)$. Now, let $n \in \mathbb{N}$. Since $\Gamma$ is a compact metric space (see Section 2.3), there exists a Borel partition $(\Gamma_i^n)_{i \in I^n}$ such that $\operatorname{diam}(\Gamma_i^n) \le 1/(n+1)$ for all $i \in I^n$. Moreover, let us define, for all $\ell \in \mathbb{N}$,

$$ B_{\ell} = \{ \gamma \in \Gamma : \operatorname{spt} \gamma' \subseteq [-\ell, \ell] \}, \quad (37) $$

$$ \text{and } \Gamma_{i, \ell}^n = \Gamma_i^n \cap (B_{\ell} \setminus B_{\ell-1}) $$

with $B_{-1} = \emptyset$. We obtain a Borel subdivision$^7$ of $\Gamma \cap (\cup_{\ell} B_{\ell})$, and possibly discarding the indices such that $\sigma(\Gamma_{i, \ell}^n) = 0$, we may rewrite (36) as

$$ G \left( \int_{\Gamma} z_{\gamma} \, d\sigma(\gamma) \right) = \sum_{i, \ell} \sigma(\Gamma_{i, \ell}^n) G \left( \frac{1}{\sigma(\Gamma_{i, \ell}^n)} \int_{\Gamma_{i, \ell}^n} z_{\gamma} \, d\sigma(\gamma) \right) = \int_{\Gamma} \tilde{G}^n(\gamma) \, d\sigma(\gamma), $$

$$ \text{where } \tilde{G}^n(\gamma) := G \left( \frac{1}{\sigma(\Gamma_{i, \ell}^n)} \int_{\Gamma_{i, \ell}^n} z_{\tilde{\gamma}} \, d\sigma(\tilde{\gamma}) \right) \text{ for } i, \ell \text{ such that } \gamma \in \Gamma_{i, \ell}^n. $$

<footnote>$^7$ There are curves in $\operatorname{Lip}_{1,d}^{\wedge}$ which are not in any $B_{\ell}$, but Smirnov's construction explicitly charges the set of curves such that $\operatorname{spt} \gamma' \subseteq [a, b]$ for some $a, b \in \mathbb{R}$ (see [36, Sec. 2.8.2]). As a result $\sigma(\Gamma \setminus \cup_{\ell} B_{\ell}) = 0$.</footnote>
Fatou's lemma then implies that

$$G\left(\int_{\Gamma} z_{\gamma} \, d\sigma(\gamma)\right) = \liminf_{n \to +\infty} \int_{\Gamma} \tilde{G}^n(\gamma) \, d\sigma(\gamma) \ge \int_{\Gamma} \liminf_{n \to +\infty} (\tilde{G}^n(\gamma)) \, d\sigma(\gamma).$$

On the other hand, arguing as in [36, Sec. 2.4.2], we note that, restricted to each $B_\ell$, the map $\gamma \mapsto z_\gamma$ is continuous for the weak-* topology. The fact that $\text{diam}(\Gamma_{i,\ell}^n) \le 1/(n+1)$ thus implies that for $\sigma$-a.e. $\gamma \in \Gamma$, if for each $n$, $\Gamma_{i,\ell}^n$ is the cell which contains $\gamma$,

$$\frac{1}{\sigma(\Gamma_{i,\ell}^n)} \int_{\Gamma_{i,\ell}^n} z_{\tilde{\gamma}} \, d\sigma(\tilde{\gamma}) \xrightarrow{*} z_{\gamma},$$

so that $\liminf_{n \to +\infty} \tilde{G}^n(\gamma) \ge G(z_\gamma)$ by the lower semi-continuity of $G$. This yields

$$G\left(\int_{\Gamma} z_{\gamma} \, d\sigma(\gamma)\right) \ge \int_{\Gamma} G(z_{\gamma}) \, d\sigma(\gamma),$$

and the converse inequality follows from Jensen's inequality. As a result, (19) holds.

To prove that $\sigma$-a.e. is a geodesic, we use a measurable selection result: there is a Borel map $\rho$ which maps any $(x, y) \in Q^2$ to some geodesic (in the sense of point (i) in the statement) $\rho_{x,y}: [0, 1] \to Q$ with $\rho_{x,y}(0) = x$, $\rho_{x,y}(1) = y$. That result follows from [30, Prop. 2.3.5] (with the minor difference that they state it for isotropic functionals $g(x, z) = g(x)$), but our Assumption (A0) makes the proof more straightforward, see Lemma 3 below.

For each $\gamma \in \bigcup_\ell B_\ell$, the map $\gamma \mapsto b(\gamma)$ (resp. $\gamma \mapsto e(\gamma)$) which associates each curve to its beginning point (resp. endpoint) is Borel. We define

$$\tilde{z} = \int z_{\rho_{b(\gamma),e(\gamma)}} d\sigma(\gamma), \quad \text{that is,} \quad \int \varphi \cdot d\tilde{z} = \int \langle z_{\rho_{b(\gamma),e(\gamma)}}, \varphi \rangle d\sigma(\gamma)$$

for all $\varphi \in C_c(\mathbb{R}^d; \mathbb{R}^d)$, one may check that $\text{div} \tilde{z} = \text{div} z$ and

$$\int g(x, \tilde{z}) \le \int \int g(\rho_{b(\gamma),e(\gamma)}(t), \rho'_{b(\gamma),e(\gamma)}(t)) dt \, d\sigma(\gamma) \\ \le \int \int g(\gamma(t), \gamma'(t)) dt \, d\sigma(\gamma) = \int g(x, z).$$

By optimality of $z$, each inequality must be an equality, hence $\sigma$-almost every $\gamma$ is a minimizing curve.

Eventually, introduce $S := \int \delta_{b(\gamma)} d\sigma(\gamma)$, $R := \int \delta_{e(\gamma)} d\sigma(\gamma)$. Both are non-negative measures, with

$$S - R = \int \text{div} z_{\gamma} \, d\sigma(\gamma) = \text{div} z.$$
Hence, the Hahn-Jordan decomposition theorem ensures that $(\operatorname{div} z)^+ \le S$ and $(\operatorname{div} z)^- \le R$. Then, Equality (18) yields

$$ (S + R)(Q) = \int \|\operatorname{div}(z_\gamma)\|_{\mathrm{TV}} d\sigma = \|\operatorname{div}(z)\|_{\mathrm{TV}} = ((\operatorname{div} z)^+ + (\operatorname{div} z)^-)(Q) $$

so that $(\operatorname{div} z)^+ = S$ and $(\operatorname{div} z)^- = R$, and point (ii) follows by considering their support. □

The following Lemma addresses the existence of a feasible point for Problem (14).

**Lemma 1.** Under Assumptions (A0), (A1), (A2), with $\mu(Q) = 0$, Problem (14) is feasible.

*Proof.* Let $\mu = \mu^+ - \mu^-$ be a Hahn-Jordan decomposition of $\mu$, and so that $\mu^+(Q) = \mu^-(Q) = |\mu|(Q)/2$. Since the case $\mu = 0$ is straightforward, we may assume, up to a rescaling, that $\mu^+(Q) = 1$. For each $n \in \mathbb{N}$, there exist discrete measures of the form $\mu_n^+ = \frac{1}{n} \sum_{i=1}^n \delta_{x_{i,n}}$, $\mu_n^- = \frac{1}{n} \sum_{j=1}^n \delta_{y_{j,n}}$ where $\{x_{i,n}\}_{i=1}^n \subseteq Q$, $\{y_{j,n}\}_{j=1}^n \subseteq Q$, and such that $\mu_n^+ \stackrel{*}{\to} \mu^+$, $\mu_n^- \stackrel{*}{\to} \mu^-$ as $n \to +\infty$.

For each pair $(x_i, y_j)$ we consider a curve $\gamma_{i,j}$ joining $x_i$ and $y_j$ satisfying (15), and we let $z^n = (\sum_{i,j} z_{\gamma_{i,j}})/n^2$. The convexity of $G: z \mapsto \int g(x, z)$ implies that

$$ D \ge \frac{1}{n^2} \sum_{i,j} G(z_{\gamma_{i,j}}) \ge G(z^n) \ge c |z^n| (\mathbb{R}^d) $$

(where the first equality follows from (A0) and the last one from (A2)). As a result we may extract a subsequence $(z_n)_{n \in I}$ (with $I \subseteq \mathbb{N}$ infinite) such that $z^n \stackrel{*}{\to} z$ for some $z \in \mathcal{M}(\mathbb{R}^d; \mathbb{R}^d)$ as $n \to +\infty$ in $I$. Then, in the weak-* sense,

$$ \operatorname{div} z = \lim_{n \to +\infty, n \in I} \operatorname{div} z^n = \lim_{n \to +\infty, n \in I} \left( \frac{1}{n} \sum_{i=1}^n \delta_{x_{i,n}} - \frac{1}{n} \sum_{j=1}^n \delta_{y_{j,n}} \right) = \mu^+ - \mu^- = \mu. $$

Moreover,

$$ |z| (\mathbb{R}^d \setminus Q) \le \liminf_{n \to +\infty, n \in I} |z^n| (\mathbb{R}^d \setminus Q) = 0, $$

so that $\operatorname{spt}(z) \subseteq Q$. Eventually, the lower semi-continuity of $G$ ensures that $G(z) \le L < +\infty$, and $z$ is thus admissible for (14). □

The next lemma shows that the decomposition of a charge in two parts implies a similar decomposition of integral functionals.

**Lemma 2.** Let $p, q \in \mathcal{V}$ such that

$$ |p+q| = |p| + |q| \qquad (38) $$

as measures, and $G(z) = \int g(x, z)$ such that (A1) holds. Then

$$ \int g(x, p+q) = \int g(x, p) + \int g(x, q). \qquad (39) $$
*Proof.* Let $z = p+q$. By assumption, $|z| = |p| + |q|$, hence $p$ and $q$ are absolutely continuous w.r.t. $|z|$. As a result of (38), for $|z|$-a.e. $x \in \mathbb{R}^d$ their densities w.r.t. $|z|$ satisfy

$$ \frac{dz}{d|z|}(x) = \frac{dp}{d|z|}(x) + \frac{dq}{d|z|}(x) \quad \text{and} \quad \left| \frac{dz}{d|z|}(x) \right| = \left| \frac{dp}{d|z|}(x) \right| + \left| \frac{dq}{d|z|}(x) \right|. $$

Hence $\frac{dp}{d|z|}(x)$ and $\frac{dq}{d|z|}(x)$ must be colinear and pointing in the same direction. Thus, there exists $\alpha \in [0, 1]$ such that

$$ \frac{dp}{d|z|}(x) = \alpha \frac{dz}{d|z|}(x) \quad \text{and} \quad \frac{dq}{d|z|}(x) = (1 - \alpha) \frac{dz}{d|z|}(x). $$

Then by the positive 1-homogeneity of $g$ (w.r.t its second argument) we get

$$ \begin{aligned} g\left(x, \frac{dz}{d|z|}(x)\right) &= \alpha g\left(x, \frac{dz}{d|z|}(x)\right) + (1-\alpha)g\left(x, \frac{dz}{d|z|}(x)\right) \\ &= g\left(x, \frac{dp}{d|z|}(x)\right) + g\left(x, \frac{dq}{d|z|}(x)\right), \end{aligned} $$

so that $\int g(x, z) = \int g(x, p) + \int g(x, q)$ (see (5)).

The next Lemma provides a measurable selection map for the geodesics.

**Lemma 3.** Assume that (A0), (A1), (A2) hold. Then, there is a Borel map $\rho: Q \times Q \to \text{Lip}^{\wedge}_{1,d}$ which maps any $(x, y) \in Q^2$ to some geodesic (for $g$) from $x$ to $y$ in $Q$. More precisely, $\rho(x, y) = \rho_{x,y}$ with $\rho_{x,y}: [0, 1] \to Q$ with $\rho_{x,y}(0) = x$, $\rho_{x,y}(1) = y$ and $\rho_{x,y}$ is a Lipschitz curve which minimizes $\int g(\gamma(t), \gamma'(t)) dt$ among the Lipschitz curves from $x$ to $y$ in $Q$).

*Proof.* The proof consists in invoking [8, Cor. 1]. By Assumptions (A0) and (A2), for each $(x, y) \in Q$, there exists some curve $\gamma$ from $x$ to $y$, with

$$ D \geq \int g(\gamma(t), \gamma'(t)) dt \geq c \int |\gamma'(t)| dt. $$

Let us fix $\ell \in \mathbb{N}$, $\ell \geq D/c$. Possibly reparametrizing $\gamma$ using arclength, we deduce that for every $(x, y) \in Q^2$, there is some $\gamma \in B_\ell$ with $\gamma(-\ell) = x$, $\gamma(\ell) = y$, and $\int g(\gamma, \gamma') \leq D$, where $B_\ell$ defined in (37). The set $B_\ell$ is compact by the Arzelà-Ascoli theorem, hence it is a complete metric space, and so is $Q^2$.

Let $\mathcal{F}: \gamma \in B_\ell \mapsto \int g(\gamma, \gamma') \in \mathbb{R}_+ \cup \{+\infty\}$. We observe that $\mathcal{F}$ is lower semi-continuous. Indeed, let $(\gamma_n)_{n \in \mathbb{N}}$ be a sequence which converges to some $\gamma \in B_\ell$ (that is, for the uniform convergence on $[-\ell, \ell]$). Possibly extracting a subsequence, we assume that $\lim_{n \to +\infty} \int g(\gamma_n, \gamma'_n) = \liminf \int g(\gamma_n, \gamma'_n)$. Then, the functions $(\gamma'_n)$ are bounded in $L^1(-\ell, \ell)$, hence we may again extract a subsequence (that is find $I \subseteq \mathbb{N}$ infinite) such that $\gamma'_n$ converges a.e. to some
$\lambda \in L^1(-\ell, \ell)$ as $n \to +\infty$ in $I$. Passing to the limit in the equality (using dominated convergence)

$$ \gamma_n(t) = \gamma_n(-\ell) + \int_{-\ell}^{t} \gamma'_n(s) \, ds $$

we see that $\lambda = \gamma'$. Then, by Fatou's lemma, $\int g(\gamma, \gamma') \le \liminf \int g(\gamma_n, \gamma'_n)$, so that $\mathcal{F}$ is lower semi-continuous.

As a result, the set

$$ \mathcal{D} = \{((x, y), \gamma) \in Q^2 \times B_{\ell} : \mathcal{F}(\gamma) < +\infty \text{ and } \gamma(-\ell) = x, \gamma(\ell) = y\} $$

is Borel, and, for each $(x, y) \in Q$, the set $\mathcal{D}_{(x,y)} = \{\gamma : ((x, y), \gamma) \in \mathcal{D}\}$ is $\sigma$-compact (it is the union of $\{\mathcal{F} \le n\} \cap \{\gamma : \gamma(-\ell) = x, \gamma(\ell) = y\}$ for $n \in \mathbb{N}$, which are closed in the compact set $B_{\ell}$). As a result, [8, Cor. 1] ensures the existence of a Borel measurable selection map $\rho$ as claimed.

## B Algorithm for the convex optimal path problem

To solve (23) we consider the saddle-point formulation (26), in which $\| \cdot \|_*$ is the dual norm of $\| \cdot \|$, defined by $\|q\|_* = \sup_{\|p\| \le 1} q \cdot p$, and $\chi$ denotes a characteristic function (0 if the condition is satisfied, $+\infty$ else). Given $\tau, \sigma > 0$ with $\tau\sigma\|A\|^2 \le 1$, the algorithm in [10] is:

1. Initialize $z^0 \in \mathcal{V}_d$, $p^0 \in (\mathbb{R}^N)^2$, set $\bar{z}^0 = z^0$

2. For each iteration $n \ge 0$, update:

$$ \left\{ \begin{aligned} p^{n+1} &= \text{Proj}_{\{p: \forall i, j, \|p_{i,j}\|_* \le g_{i,j}\}} (p^n + \sigma A \bar{z}^n) \\ z^{n+1} &= \text{Proj}_{\{D^*z=\mu\}} (z^n - \tau A^* p^{n+1}) \\ \bar{z}^{n+1} &= 2z^{n+1} - z^n \end{aligned} \right. $$

where the first step requires to project onto $\{p : \forall i, j, \|p_{i,j}\|_* \le g_{i,j}\}$, which amounts to project independently each component on a 2D disc, and

$$ \text{Proj}_{\{D^*z=\mu\}}(z) = z + D(D^*D)^{-1}D^*(z_0 - z) $$

is the projection onto the space $D^*z = \mu$, with $z_0$ being any vector field such that $D^*z_0 = \mu$. The discrete Neumann Laplacian $(D^*D)$ (which is invertible on functions with zero average) is diagonalized and inverted by means of a DCT (which turns out to be consistent with the no-flux condition on the boundary for vector fields), except in the “roto-translational” representation (Sec. 5) where the block is periodic in the third component, and one has to use a FFT for that one. We relied on the Python bindings `pyFFTW` for the `FFTW3` [22] library to compute the Fourier transforms. The extension described in Section 5 just requires to adapt the first projection step to the corresponding anisotropic energy, replacing it with the formulas (35), for the various choices of $\bar{h}$.
For bounded sets

$$B_1 \subset \{z \in \mathcal{V}_d, D^*z = \mu\} \quad \text{and} \quad B_2 \subset \{p \in (\mathbb{R}^N)^2, \forall i, j, \|p_{i,j}\|_* \le g_{i,j}\},$$

we consider the partial gap:

$$\mathcal{G}_{B_1 \times B_2}(z, p) = \max_{p' \in B_2} \langle Az, p' \rangle - \min_{z' \in B_1} \langle Az', p \rangle.$$

It is shown [10] that:

$$\mathcal{G}_{B_1 \times B_2}(\bar{z}_k, \bar{p}_k) \le \frac{D(B_1, B_2)}{k} \qquad (40)$$

where we have denoted $\bar{z}_k = \sum_{i=1}^k z^i/k$ and $\bar{p}_k = \sum_{i=1}^k p^i/k$ and where

$$D(B_1, B_2) = \max_{(z,p) \in B_1 \times B_2} \frac{\|z - z^0\|_2^2}{2\tau} + \frac{\|p - p^0\|_2^2}{2\sigma} \qquad (41)$$

Since the z variable is expected to be a superposition of curves, one expects $\|z - z^0\|_2^2 \sim \sqrt{N \times M}$ (it behaves like the characteristic length of the domain), and since the dual variable represents the (uniformly bounded) gray levels values of the whole image, $\|p - p^0\|_2^2 \sim N \times M$. Hence, optimizing for $\tau, \sigma$ in (41) under the constraint $\tau\sigma\|A\|^2 \le 1$, we set in practice $\tau = .99/(\|A\|(N \times M)^{1/4})$ and $\sigma = (N \times M)^{1/4}/\|A\|$. With this choice, $D(B_1, B_2) \lesssim \|A\|(N \times M)^{3/4}$ (the norm of the averaging operator A is of order 1).

For images of similar size as the one in Fig. 6 (200 × 200 pixels), only about 60 steps of the primal-dual algorithm are needed at each motion of the Dirac masses. Practical convergence is a lot faster than the $O(1/k)$ rate guaranteed in (40), as for a 200 × 200 image, that would require thousands of iterations (estimating $D(B_1, B_2) \sim 3000$). It seems that the “preconditioning” induced by the choice of projecting onto the divergence constraint, while not changing much the theoretical bound, improves drastically the practical convergence.