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

HAL
open science

# A variational method for curve extraction

Majid Arthaud, Antonin Chambolle, Vincent Duval

► To cite this version:

Majid Arthaud, Antonin Chambolle, Vincent Duval. A variational method for curve extraction. SSVM 2025 - 10th Scale-Space and Variational Methods in Computer Vision, May 2025, Totnes, Devon, United Kingdom. hal-04980911

HAL Id: hal-04980911

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

Submitted on 6 Mar 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

Majid Arthaud¹,², Antonin Chambolle³,²[0000-0002-9465-4659], and Vincent Duval²,³[0000-0002-7709-256X]

¹ ENPC, 6 Av. Blaise Pascal, 77420 Champs-sur-Marne, France
majid.arthaud@eleves.enpc.fr

² INRIA Mokaplan, Paris, France
vincent.duval@inria.fr

³ CEREMADE, Université Paris-Dauphine, PSL, Paris, France
antonin.chambolle@ceremade.dauphine.fr

**Abstract.** We propose a functional for extracting curves between a list of possible endpoints, based on a discretization of a variational energy and Smirnov's decomposition theorem for vector fields. It is then used to design a bi-level minimization approach to automatically extract curves and 1D structures from an image, which is mostly unsupervised.

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

# 1 Introduction

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 [14,2,6] usually rely on the minimization of a potential along a curve, with functionals of the form

$$ \int_{\Gamma} g(x) 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). 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 (geodesic) joining these two points [6,8]. See for instance [21] (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 [5]), 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 [17,18] for the tracking
of the width of retina vessels, or as in [9] for their local orientation). Alternatively, the problem of computing geodesics may be reformulated as a primal-dual problem as in [11].

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 [15,16]. Consider for instance the problem:

$$ \min_{z} \left\{ \int g(x) d|z|(x) : -\operatorname{div} z = \delta_{B} - \delta_{A} \right\} \qquad (2) $$

where $z$ is a measure vector field of total variation measure $|z|$, and $A$ and $B$ are two given endpoints. A typical example of such a measure field is given by the integration along curves, that is, measures of the form $\tau H^1|\Gamma$, where $\tau$ is tangent to $\Gamma$ and $\Gamma$ is a curve. We show in Proposition 1 below (thanks to the celebrated Smirnov theorem on extreme points of such vector fields [22]) 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 [21], 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, atomic measures with vanishing total mass, as will be considered 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. In addition, one may incorporate this setting in variational formulations for more general inverse problems (such as a “deconvolution” of thin structures) as in [15,16].

Before describing our setting 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. We discuss its numerical implementation and show results in Section 5, before some concluding remarks.

## 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 image processing community in [15]. Given an compact arcwise connected set $Q \subseteq \mathbb{R}^2$, 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}^2)^2 : \operatorname{spt}(z) \subseteq Q, \operatorname{div}(z) \in \mathcal{M}(\mathbb{R}^2)\}, \qquad (3) $$
where $\mathcal{M}(\mathbb{R}^2)$ is the space of finite (signed) Radon measures over $\mathbb{R}^2$. It is a normed space when equipped with the norm

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

where $\| \cdot \|_{\text{TV}}$ denotes the total variation of (vector or scalar) Radon measures. An important example of a normal charge is induced by simple oriented curves, that is Lipschitz functions $\gamma: [0, 1] \to Q$ such that $\gamma$ is one-to-one on $[0, 1)$. The charge $T_\gamma$ is then the vectorial measure defined by

$$ \forall \varphi \in C_c(\mathbb{R}^2; \mathbb{R}^2), \quad \langle T_\gamma, \varphi \rangle = \int_0^1 \varphi(\gamma(t)) \cdot \gamma'(t) dt. \qquad (5) $$

One may check that $\text{div}(T_\gamma) = \delta_{\gamma(0)} - \delta_{\gamma(1)}$, so that $\text{div}(T_\gamma) = 0$ if and only if the curve is closed. Smirnov's decomposition theorem [22] shows that any normal charge $z$ can be decomposed as a superposition of charges induced by curves, and that a similar decomposition holds for the norm $\|z\|_{\mathcal{V}}$.

An alternative interpretation of Smirnov's results is that the extreme points of the unit ball of $\| \cdot \|_{\mathcal{V}}$ are the measures supported on simple oriented rectifiable curves, of norm one (in higher dimension, this remains mostly true up to a generalization of the notion of curves, allowing for infinitely winding curves). It has led Laville et al. [15,16] 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.

Here, we focus on solving problems of the form

$$ \min_z \left\{ \int_Q g(x) d|z|(x) : -\text{div} z = \mu \right\} \qquad (6) $$

where $\mu$ is a given measure with support in $Q$ and $\mu(Q) = \int_Q d\mu = 0$, and $g$ is a continuous weight, bounded from above and below by a positive constant. This kind of problem also appears in the field of branched transportation (see for instance [19]). The following result is a continuous version of known and standard results in graph theory and network flows (see also the discussion below (10)).

**Proposition 1.** Let $z$ be a minimizer in (6) for $\mu \in \mathcal{M}(\mathbb{R}^2)$ with support in $Q$ and $\mu(Q) = 0$. There exists a Borel probability measure $\sigma$ defined on the set of open simple oriented curves such that

$$ z = \int T_{\gamma} d\sigma(\gamma), \quad \|z\|_{\text{TV}} = \int \|T_{\gamma}\|_{\text{TV}} d\sigma(\gamma), \quad \|\text{div}(z)\|_{\text{TV}} = \int \|\text{div}(T_{\gamma})\|_{\text{TV}} d\sigma(\gamma), \qquad (7) $$

and for $\sigma$-a.e. $\gamma$

(i) $\gamma$ is a geodesic curve, i.e. minimizes for $\int g(\gamma(t)) |\gamma'(t)| dt$ for fixed $\gamma(0)$ and $\gamma(1)$,

(ii) $\gamma(0) \in \text{spt } \mu^-$ and $\gamma(1) \in \text{spt } \mu^+$ where $\mu = \mu^+ - \mu^-$ is the Hahn-Jordan decomposition of $\mu = -\text{div } z$.
Proof. We only sketch the proof as it essentially follows the proof of [19, Thm. 1.3.2] with $S = Q$ (that set is not countably 1-rectifiable but this assumption is not needed for our result). First, Smirnov’s theorem [22, Thm. C] implies that $z$ (“totally”) decomposes into two parts, $P$ and $Q$, with $\text{div} P = 0$ and $Q$ is decomposed on curves joining points in the support of $\mu$. The minimality of $z$ for (6) then implies that $P = 0$. Hence $z = Q$, and the same theorem ensures that $z$ can be decomposed as in (7) for some Borel measure $\sigma$. Note that in fact, the three equalities in (7) are in the sense of measures, and integrating the second one against the continuous function $g$, one finds that

$$ \int \int g(\gamma(t))|\gamma'(t))|dt d\sigma(\gamma) = \int gd|z|. $$

By [19, Prop. 2.3.5] (see also [10, Lemma 3.4]), 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$. Then, defining

$$ \tilde{z} = \int T_{\rho_{\gamma(0),\gamma(1)}} d\sigma(\gamma), \quad \text{that is,} \quad \int \varphi \cdot d\tilde{z} = \int \langle T_{\rho_{\gamma(0),\gamma(1)}}, \varphi \rangle d\sigma(\gamma) $$

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

$$ \begin{aligned} \int gd|\tilde{z}| &\le \int \int g(\rho_{\gamma(0),\gamma(1)}(t))|\rho'_{\gamma(0),\gamma(1)}(t))|dt d\sigma(\gamma) \\ &\le \int \int g(\gamma(t))|\gamma'(t))|dt d\sigma(\gamma) = \int gd|z|. \end{aligned} $$

By optimality of $z$, each inequality must be an equality, hence $\sigma$-almost every $\gamma$ is a geodesic curve. Eventually, introduce $S := \int \delta_{\gamma(0)}d\sigma(\gamma)$, $R := \int \delta_{\gamma(1)}d\sigma(\gamma)$. Both are non-negative measures, with

$$ S - R = \int \text{div} T_{\gamma} d\sigma(\gamma) = \text{div} z. $$

Hence, the Hahn-Jordan decomposition theorem ensures that $(\text{div} z)^+ \le S$ and $(\text{div} z)^- \le R$. Then, Equality (7) yields

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

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

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 [16], 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.
## 3 Discrete curves

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).

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

$$ \mathcal{V}_d = \mathbb{R}^{\mathcal{E}} \simeq (\mathbb{R})^{(N-1) \times M} \times (\mathbb{R})^{N \times (M-1)} \quad (8) $$

For a vector field $v \in \mathcal{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}^{\mathcal{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}. \quad (9) $$

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}^{\mathcal{N}}$, given by $\langle D^*v, u \rangle_{\mathbb{R}^{\mathcal{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 on the boundary, $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 (6) is

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

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

Nevertheless, as we illustrate below, 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. One has therefore to turn to other norms, consistent (as $M, N \to \infty$) with the Euclidean norm considered in (6), even if in doing so we lose, at the discrete level, the decomposition theorem of a field into superposition of curves. In the active contour community, that phenomenon is well known, and people prefer to use the Eikonal equation together with fast marching approaches over Dijkstra’s algorithm on the graph of the image [6].
Therefore, we introduce a family of discrete approximations of the convex curve reconstruction problem (6), 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 some image), and the scalar field of curve endpoints $\mu \in \mathbb{R}^N$, the problem is

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

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 (10)) 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 (12) $$

By solving (11) we identify, through the discrete vector field z, dark curves (i.e. where g is close to zero) on a light background (i.e. where g is close to one). An experimental result with (12) 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. All of the experimental results exposed in this document have been implemented with a primal-dual algorithm [3], accelerated by a reprojection on the constraint $D^*z = \mu$ computed with FFTW3 [12], see the description in the supplementary material.

A second instantiation, yielding sharper and more isotropic results, consists in replacing the $\ell^1$-norm in (12) with the Euclidean norm $\|x\| = \sqrt{x_1^2 + x_2^2}$. This corresponds to a consistent discretization of (6) 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 [13], 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 (12) 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 (11). 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 (13) $$
<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Left: 200 × 200 “noisy comma” shape (crop from [16, Fig. 1]), middle: result with weighted $\ell^1$ norm, right: weighted $\ell^2$ norm.</figcaption></figure>

Notice that here, we have defined the coefficients outside of the image with zero padding (i.e. $z_{1/2,j} = z_{N+1/2,j} = z_{i,1/2} = z_{i,M+1/2} = 0$), which is consistent with the vanishing flux condition on the boundary. The results are very isotropic and 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 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 [8].

<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>

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

# 4 An iterative discrete curve reconstruction algorithm

The discrete convex curve reconstruction problem of the previous section takes for granted the fact that the endpoints (i.e. discrete Dirac masses) of the curves
<figure><img src="image_5.png" /><figcaption>Fig. 3. Left: an angiogram (source: Wikipedia), right: geodesic computed by our method (compare for instance with [8, Fig. 4(d)]).</figcaption></figure>

to reconstruct were given. In this section, we propose an iterative algorithm so as to retrieve discrete curves in images with no prior knowledge of their geometry. 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 (g(x) - g_{max}) d|z|(x) : z \text{ solves (6) for } \mu = \sum_{i=1}^n \delta_{t_i} - \sum_{i=1}^n \delta_{s_i} \right\} \quad (14) $$

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 (14). 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 q$ improves the criterion. Also, we sometimes need to decrease the number of sources/sinks.

We now describe the discrete implementation. 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 (11) (e.g. with the primal-dual algorithm [3], see the supplementary material);
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 \times 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 \times 3$ square around a Dirac mass are detailed in Fig. 4. 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.

<figure><img src="image_6.png" /><figcaption>**Fig. 4.** 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>

# 5 Numerical results

We show a few results obtained with our implementation, which is available at https://github.com/majidarth/Discrete-curve-reconstruction. The computation
<figure><img src="image_7.png" /><figcaption>**Fig. 5.** Result of the iterative algorithm, displayed after post-processing, on a 200 × 200 image of 42 bird chromosomes (source: rescaling of an image from [20]). This result is obtained after 100 iterations of the Sec. 4 algorithm, each iteration corresponding to 60 primal-dual steps, for a total of 1 min. 05 s. wall time on an average laptop. Post-processing, which concludes with 5000 primal-dual steps on the resulting set of Dirac masses, adds another 1 min. 12 s. wall time.</figcaption></figure>

of the optimal $z$ given $\mu$ relies on a primal-dual method [3], detailed in the supplementary material. In practice, getting a good estimate for $z$ which solves (11) 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 $\mu$ does not change much from one step to the next). Moreover, we use simple projection steps and an adapted choice of the step sizes. See the supplementary material for details. The implementation is developed in python and does not use parallelization. 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. 5, $\mu$ was initialized as 2000 such pairs of Dirac masses, while there are only 42 chromosomes in the original image.

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, (11) is then computed one last time, and the resulting *z* is the final result. In addition to Fig. 5, we show three results on the “comma” image on Fig. 6: 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.

<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Result of the iterative algorithm on the 200 × 200 “comma” shape. The iteration number *k* for the algorithm of Sec. 4 is indicated. Top: *μ* is initialized as a pair of Dirac masses of opposite intensities randomly chosen in the sublevel set {g ≤ gmax}. Middle: *μ* is initialized as 15 pairs of Dirac masses of opposite intensities randomly chosen in the sublevel set {g ≤ gmax}. Bottom: same as “middle”, but the “comma” shape is blurred with a Gaussian filter of standard deviation 2. At each iteration, 80 primal-dual steps are computed. The middle experiment was the longest, taking 5 min. 07 s. wall time after 390 iterations on an average laptop. Post-processing (with 5000 primal-dual steps on the resulting set of Dirac masses) added 1 min. 03 s. wall time.</figcaption></figure>
# 6 Conclusion

In this paper, we have introduced a variational method for curve extraction which differs from the usual Eikonal equation approach [21], as we represent paths as measure vector fields. It allows to simultaneously compute geodesics between a large numbers of endpoints. We also propose a bi-level method to automatically fit the endpoints, yielding a complete, mostly unsupervised automatic curve extraction algorithm. Further work could focus on more refined averaging to correct anisotropy, analyzing or reproducing recent findings in the context of discrete total variations. We also intend to extend the approach to the more challenging 3D setting.

# Acknowledgements

Antonin Chambolle acknowledges the support of the “France 2030” funding ANR-23-PEIA-0004 (“PDE-AI”).

# References

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

2. 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

3. 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

4. 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

5. 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

6. Cohen, L.D., Kimmel, R.: Global minimum for active contour models: A minimal path approach. International Journal of Computer Vision **24**, 57–78 (1996)

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

8. 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)

9. 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

10. Duval, V., Tovey, R.: Dynamical programming for off-the-grid dynamic inverse problems. ESAIM Control Optim. Calc. Var. **30**, Paper No. 7, 40 (2024). https://doi.org/10.1051/cocv/2023085, https://doi.org/10.1051/cocv/2023085
11. 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

12. 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

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

14. 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

15. 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

16. 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

17. 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

18. 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

19. 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

20. 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

21. 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

22. 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)
# A variational method for curve extraction: supplementary material

Majid Arthaud, Antonin Chambolle, and Vincent Duval

## A Algorithm for the convex optimal path problem

To solve (11) we consider the following primal-dual formulation:

$$ \min_{z \in \mathcal{V}_d} \max_{p \in (\mathbb{R}^N)^2} \langle Az, p \rangle + \chi_{\{D^*z=\mu\}} - \chi_{\{\forall i,j, \|p_{i,j}\|_* \le g_{i,j}\}} $$

where $\| \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 [3] 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:

$$ \begin{cases} 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{cases} $$

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_0 + D(D^*D)^{-1}D^*(z - z_0) $$

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 inverted by means of a DCT (we use FFTW3 [12]). 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 [3] that:

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

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} \quad (2) $$
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 (2) 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. 5, i.e. 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 (1), 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.