Functional Liftings of Vectorial Variational Problems with Laplacian Regularization

TL;DR

This paper introduces a convex relaxation method based on functional lifting for vectorial variational problems with Laplacian regularization, enabling advanced image processing applications like 2D image registration.

Contribution

It develops a novel convex relaxation framework that extends existing multilabeling approaches to higher-order regularization and vectorial data, with theoretical and practical insights.

Findings

Encompasses sublabel-accurate multilabeling as a special case

Provides a mathematical connection between lifted and original functionals

Successfully applied to 2D image registration problems

Abstract

We propose a functional lifting-based convex relaxation of variational problems with Laplacian-based second-order regularization. The approach rests on ideas from the calibration method as well as from sublabel-accurate continuous multilabeling approaches, and makes these approaches amenable for variational problems with vectorial data and higher-order regularization, as is common in image processing applications. We motivate the approach in the function space setting and prove that, in the special case of absolute Laplacian regularization, it encompasses the discretization-first sublabel-accurate continuous multilabeling approach as a special case. We present a mathematical connection between the lifted and original functional and discuss possible interpretations of minimizers in the lifted function space. Finally, we exemplarily apply the proposed approach to 2D image registration problems.

Figures

Equations

$$F(u)={\int }_{\Omega }f(x,u(x),\Delta u(x))dx,$$

$$f(x,z,p)=\rho (x,z)+\eta (p).$$

$$\underset{u}{\min }F(u)={\int }_{\Omega }f(x,u(x),\nabla u(x))dx,$$

$$F(u)={\int }_{\Omega \times \mathbb{R}}\varphi \cdot D{1}_u,$$

$$\mathcal{F}(v):=\underset{\varphi \in X}{\sup }{\int }_{\Omega \times \mathbb{R}}\varphi \cdot Dv,$$

$$\mathcal{C}=\{v:\Omega \times \mathbb{R}\to [0,1]:\underset{z\to -\infty }{\lim }v(x,z)=1,\underset{z\to \infty }{\lim }v(x,z)=0\}.$$

$$X=\{\varphi :\Omega \times \mathbb{R}\to {\mathbb{R}}^{d+1}:$$

$$\text{for every }(x,z)\in \Omega \times \mathbb{R}\},$$

$$F(u)=\mathcal{F}(1_u)$$

$$\mathcal{F}(\mathbf{u}):=\underset{(p,q)\in X}{\sup }{\int }_{\Omega }{\int }_{\Gamma }({\Delta }_{x}p(x,z)+q(x,z))d{\mathbf{u}}_x(z)dx,$$

$$X=\{(p,q):$$

$$z\mapsto p(x,z)\text{ concave }$$

$$\text{and }q(x,z)+f^{\ast }(x,z,{\nabla }_{z}p(x,z))\le 0$$

$$\text{for every }(x,z)\in \Omega \times \Gamma \},$$

$$F(u)\ge \mathcal{F}(\mathbf{u}).$$

$${\Delta }_{x}p(x,u(x))$$

$$-2⟨{\nabla }_x{\nabla }_{z}p(x,u(x)),\nabla u(x)⟩-⟨{\nabla }_{z}p(x,u(x)),\Delta u(x)⟩.$$

$$-{\int }_{\Omega }⟨{\nabla }_x{\nabla }_{z}p(x,u(x)),\nabla u(x)⟩dx$$

$$+{\int }_{\Omega }\sum _{i=1}^d⟨{\partial }_{i}u(x),D_z^{2}p(x,u(x)){\partial }_{i}u(x)⟩dx,$$

$${\int }_{\Omega }{\Delta }_{x}p(x,u(x))dx\le {\int }_{\Omega }⟨{\nabla }_{z}p(x,u(x)),\Delta u(x)⟩dx.$$

$$\mathcal{F}(\mathbf{u})$$

$$={\int }_{\Omega }{\Delta }_{x}p(x,u(x))+q(x,u(x))dx$$

$$\stackrel{\text{eqref}eq:fstar-ineq}{\le }{\int }_{\Omega }{\Delta }_{x}p(x,u(x))-f^{\ast }(x,u(x),{\nabla }_{z}p(x,u(x)))dx$$

$$\stackrel{\text{eqref}eq:lapl-est}{\le }{\int }_{\Omega }⟨{\nabla }_{z}p(x,u(x)),\Delta u(x)⟩-f^{\ast }(x,u(x),{\nabla }_{z}p(x,u(x)))dx$$

$$\le {\int }_{\Omega }f(x,u(x),\Delta u(x))dx,$$

$${\nabla }_{z}p(x,u(x))\in {\partial }_{p}f(x,u(x),\Delta u(x))$$

$$p(x,z)$$

$$q(x,z)$$

$$(\Delta u{)}^{i_0}={\sum }_{l=1}^m(u^{i_l}-u^{i_0})\in {\mathbb{R}}^s.$$

$${\sum }_{l=1}^4(u^{i_l}-u^{i_0})$$

$$\varphi (p):=\{\begin{array}{ll}\mu \Vert {\sum }_{l=1}^m(T_{j_l}{\alpha }^l-T_{j_0}{\alpha }^0)\Vert ,& \text{if }p=\mu {\sum }_{l=1}^m(P^{j_l}{\alpha }^l-P^{j_l}{\alpha }^0),\\ +\infty ,& \text{otherwise,}\end{array}$$

$$\mathcal{K}={\bigcap }_{1\le j_l\le M}\{f\in {\mathbb{R}}^L:$$

$$\text{for any }{\alpha }^l\in {\Delta }_{s+1}^U,l=0,1,\dots ,m\},$$

$$\mathcal{K}=\left\{f\in {\mathbb{R}}^L:f:\Gamma \to \mathbb{R}\text{ is concave and 1-Lipschitz continuous}\right\}.$$

$$\frac{1}{m}\sum _{l=1}^m(f(t^l)-f(t^0))$$

$$\le \Vert \left(\frac{1}{m}\sum _{l=1}^m{t}^l\right)-t^0\Vert =\frac{1}{m}\Vert \sum _{l=1}^m(t^l-t^0)\Vert .$$

$$\frac{f(t^1)+f(t^2)}{2}\le f\left(\frac{t^1+t^2}{2}\right)$$

Full text

1 Introduction

Let $\Omega\subset\mathbb{R}^{d}$ and $\Gamma\subset\mathbb{R}^{s}$ both be bounded sets. In the following, we consider the variational problem of minimizing the functional

[TABLE]

that acts on vector-valued functions $u\in C^{2}(\Omega;\Gamma)$. Convexity of the integrand $f\colon\Omega\times\Gamma\times\mathbb{R}^{s}\to\mathbb{R}$ is only assumed in the last entry, so that $u\mapsto F(u)$ is generally non-convex. The Laplacian $\Delta u$ is understood component-wise and reduces to $u^{\prime\prime}$ if the domain $\Omega$ is one-dimensional.

Variational problems of this form occur in a wide variety of image processing tasks, including image reconstruction, restoration, and interpolation. Commonly, the integrand is split into data term and regularizer:

[TABLE]

As an example, in image registration (sometimes referred to as large-displacement optical flow), the data term $\rho(x,z)=d(R(x),T(x+z))$ encodes the pointwise distance of a reference image $R\colon\mathbb{R}^{d}\to\mathbb{R}^{k}$ to a deformed template image $T\colon\mathbb{R}^{d}\to\mathbb{R}^{k}$ according to a given distance measure $d(\cdot,\cdot)$, such as the squared Euclidean distance $d(a,b)=\frac{1}{2}\|a-b\|_{2}^{2}$. While often a suitable convex regularizer $\eta$ can be found, the highly non-convex nature of $\rho$ renders the search for global minimizers of (1) a difficult problem.

Instead of directly minimizing $F$ using gradient descent or other local solvers, we will aim to replace it by a convex functional $\mathcal{F}$ that acts on a higher-dimensional (lifted) function space. If the lifting is chosen in such a way that we can construct global minimizers of $F$ from global minimizers of $\mathcal{F}$, we can find a global solution of the original problem by applying convex solvers to $\mathcal{F}$. While we cannot claim this property for our choice of lifting, we believe that the mathematical motivation and some of the experimental results show that this approach can be a good basis for future work on global solutions of variational models with higher-order regularization.

Calibrations in variational calculus  The lifted functional $\mathcal{F}$ proposed in this work is motivated by previous lifting approaches for first-order variational problems of the form

[TABLE]

where $F$ acts on functions $u\colon\Omega\to\Gamma$ with $\Omega\subset\mathbb{R}^{d}$ and scalar range $\Gamma\subset\mathbb{R}$.

The calibration method as introduced in [1] gives a globally sufficient optimality condition for functionals of the form (3) with $\Gamma=\mathbb{R}$. Importantly, $f(x,z,p)$ is not required to be convex in $(x,z)$, but only in $p$. The method states that $u$ minimizes $F$ if there exists a divergence-free vector field $\phi\colon\Omega\times\mathbb{R}\to\mathbb{R}^{d+1}$ (a calibration) in a certain admissible set $X$ of vector fields on $\Omega\times\mathbb{R}$ (see below for details), such that

[TABLE]

where $\mathbf{1}_{u}$ is the characteristic function of the subgraph of $u$ in $\Omega\times\mathbb{R}$, $\mathbf{1}_{u}(x,z)=1$ if $u(x)>z$ and [math] otherwise, and $D\mathbf{1}_{u}$ is its distributional derivative. The duality between subgraphs and certain vector fields is also the subject of the broader theory of Cartesian currents [8].

A convex relaxation of the original minimization problem then can be formulated in a higher-dimensional space by considering the functional [5, 19]

[TABLE]

acting on functions $v$ from the convex set

[TABLE]

In both formulations, the set of admissible test functions is

[TABLE]

where $f^{*}(x,z,p):=\sup_{q}\langle p,q\rangle-f(x,z,q)$ is the convex conjugate of $f$ with respect to the last variable. In fact, the equality

[TABLE]

has been argued to hold for $u\in W^{1,1}(\Omega)$ under suitable assumptions on $f$ [19]. A rigorous proof of the case of $u\in BV(\Omega)$ and $f(x,z,p)=f(z,p)$ ($f$ independent of $x$), but not necessarily continuous in $z$, can be found in the recent work [2].

In [17], it is discussed how the choice of discretization influences the results of numerical implementations of this approach. More precisely, motivated by the work [18] from continuous multilabeling techniques, the choice of piecewise linear finite elements on $\Gamma$ was shown to exhibit so-called sublabel-accuracy, which is known to significantly reduce memory requirements.

Vectorial data  The application of the calibration method to vectorial data $\Gamma\subset\mathbb{R}^{s}$, $s>1$, is not straightforward, as the concept of subgraphs, which is central to the idea, does not translate easily to higher-dimensional range. While the original sufficient minimization criterion has been successfully translated [16], functional lifting approaches have not been based on this generalization so far. In [20], this approach is considered to be intractable in terms of memory and computational performance.

There are functional lifting approaches for vectorial data with first-order regularization that consider the subgraphs of the components of $u$ [9, 21]. It is not clear how to generalize this approach to nonlinear data $\Gamma\subset\mathcal{M}$, such as a manifold $\mathcal{M}$, where other functional lifting approaches exist at least for the case of total variation regularization [14].

An approach along the lines of [18] for vectorial data with total variation regularization was proposed in [12]. Even though [17] demonstrated how [18] can be interpreted as a discretized version of the calibration-based lifting, the equivalent approach [12] for vectorial data lacks a fully-continuous formulation as well as a generalization to arbitrary integrands that would demonstrate the exact connection to the calibration method.

Higher-order regularization  Another limitation of the calibration method is its limitation to first-order derivatives of $u$, which leaves out higher-order regularizers such as the Laplacian-based curvature regularizer in image registration [7]. Recently, a functional lifting approach has been successfully applied to second-order regularized image registration problems [15], but the approach was limited to a single regularizer, namely the integral over the $1$-norm of the Laplacian (absolute Laplacian regularization).

Projection of lifted solutions  In the scalar-valued case with first-order regularization, the calibration-based lifting is known to generate minimizers that can be projected to minimizers of the original problem by thresholding [19, Theorem 3.1]. This method is also used for vectorial data with component-wise lifting as in [21]. In the continuous multi-labeling approaches [14, 18, 12], simple averaging is demonstrated to produce useful results even though no theoretical proof is given addressing the accuracy in general. In convex LP relaxation methods, projection (or rounding) strategies with provable optimality bounds exist [11] and can be extended to the continuous setting [13]. We demonstrate that rounding is non-trivial in our case, but will leave a thorough investigation to future work.

Contribution  In Section 2, we propose a calibration method-like functional lifting approach in the fully-continuous vector-valued setting for functionals that depend in a convex way on $\Delta u$. We show that the lifted functional satisfies $\mathcal{F}(\delta_{u})\leq F(u)$, where $\delta_{u}$ is the lifted version of a function $u$ and discuss the question of whether the inequality is actually an equality. For the case of absolute Laplacian regularization, we show that our model is a generalization of [15]. In Section 2.3, we clarify how convex saddle-point solvers can be applied to our discretized model. Section 3 is concerned with experimental results. We discuss the problem of projection and demonstrate that the model can be applied to image registration problems.

2 A calibration method with vectorial second-order terms

2.1 Continuous formulation

We propose the following lifted substitute for $F$:

[TABLE]

acting on functions $\mathbf{u}\colon\Omega\to\mathcal{P}(\Gamma)$ with values in the space $\mathcal{P}(\Gamma)$ of Borel probability measures on $\Gamma$. This means that, for each $x\in\Omega$ and any measurable set $U\subset\Gamma$, the expression $\mathbf{u}_{x}(U)\in\mathbb{R}$ can be interpreted as the “confidence” of an assumed underlying function on $\Omega$ to take a value inside of $U$ at point $x$. A function $u\colon\Omega\to\Gamma$ can be lifted to a function $\mathbf{u}\colon\Omega\to\mathcal{P}(\Gamma)$ by defining $\mathbf{u}_{x}:=\delta_{u(x)}$, the Dirac mass at $u(x)\in\Gamma$, for each $x\in\Omega$.

We propose the following set of test functions in the definition of $\mathcal{F}$:

[TABLE]

where $f^{*}(x,z,q):=\sup_{p\in\mathbb{R}^{s}}\langle q,p\rangle-f(x,z,p)$ is the convex conjugate of $f$ with respect to the last argument.

A thorough analysis of $\mathcal{F}$ requires a careful choice of function spaces in the definition of $X$ as well as a precise definition of the properties of the integrand $f$ and the admissible functions $\mathbf{u}\colon\Omega\to\mathcal{P}(\Gamma)$, which we leave to future work. Here, we present a proof that the lifted functional $\mathcal{F}$ bounds the original functional $F$ from below.

Proposition 1

Let $f\colon\Omega\times\Gamma\times\mathbb{R}^{s}\to\mathbb{R}$ be measurable in the first two, and convex in the third entry, and let $u\in C^{2}(\Omega;\Gamma)$ be given. Then, for $\mathbf{u}\colon\Omega\to\mathcal{P}(\Gamma)$ defined by $\mathbf{u}_{x}:=\delta_{u(x)}$, it holds that

[TABLE]

Proof

Let $p,q$ be any pair of functions satisfying the properties from the definition of $X$. By the chain rule, we compute

[TABLE]

Furthermore, the divergence theorem ensures

[TABLE]

as well as $\int_{\Omega}\Delta\left[p(x,u(x))\right]dx=0$ by the compact support of $p$. As $p\in C_{c}^{2}(\Omega\times\Gamma)$, concavity of $z\mapsto p(x,z)$ implies a negative semi-definite Hessian $D^{2}_{z}p(x,z)$, so that, together with (16)–(17),

[TABLE]

We conclude

[TABLE]

where we used the definition of $f^{*}$ in the last inequality. ∎

By a standard result from convex analysis, $\langle p,g\rangle-f^{*}(x,z,g)=f(x,z,p)$ whenever $g\in\partial_{p}f(x,z,p)$, the subdifferential of $f$ with respect to $p$. Hence, for equality to hold in (15), we would need to find a function $p\in C_{c}^{2}(\Omega\times\Gamma)$ with

[TABLE]

and associated $q(x,z):=-f^{*}(x,z,\Delta u(x))$, such that $(p,q)\in X$ or $(p,q)$ can be approximated by functions from $X$.

Separate data term and regularizer  If the integrand can be decomposed into $f(x,z,p)=\rho(x,z)+\eta(p)$ as in (2), with $\eta\in C^{1}(\mathbb{R}^{s})$ and $u$ sufficiently smooth, the optimal pair $(p,q)$ in the sense of (24) can be explicitly given as

[TABLE]

A rigorous argument that such $p,q$ exist for any given $u$ could be made by approximating them by compactly supported functions from the admissible set $X$ using suitable cut-off functions on $\Omega\times\Gamma$.

2.2 Connection to the discretization-first approach [15]

In [15], data term $\rho$ and regularizer $\eta$ are lifted independently from each other for the case $\eta=\|\cdot\|_{1}$. Following the continuous multilabeling approaches in [6, 18, 12], the setting is fully discretized in $\Omega\times\Gamma$ in a first step. Then the lifted data term and regularizer are defined to be the convex hull of a constraint function, which enforces the lifted terms to agree on the Dirac measures $\delta_{u}$ with the original functional applied to the corresponding function $u$. The data term is taken from [12], while the main contribution concerns the regularizer that now depends on the Laplacian of $u$.

In this section, we show that our fully-continuous lifting is a generalization of the result from [15] after discretization.

Discretization  In order to formulate the discretization-first lifting approach given in [15], we have to clarify the used discretization.

For the image domain $\Omega\subset\mathbb{R}^{d}$, discretized using points $X^{1},\dots,X^{N}\in\Omega$ on a rectangular grid, we employ a finite-differences scheme: We assume that, on each grid point $X^{i_{0}}$, the discrete Laplacian of $u\in\mathbb{R}^{N,s}$, $u^{i}\approx u(X^{i})\in\mathbb{R}^{s}$, is defined using the values of $u$ on $m+1$ grid points $X^{i_{0}},\dots,X^{i_{m}}$ such that

[TABLE]

For example, in the case $d=2$, the popular five-point stencil means $m=4$ and the $X^{i_{l}}$ are the neighboring points of $X^{i_{0}}$ in the rectangular grid. More precisely,

[TABLE]

The range $\Gamma\subset\mathbb{R}^{s}$ is triangulated into simplices $\Delta_{1},\dots,\Delta_{M}$ with altogether $L$ vertices (or labels) $Z^{1},\dots,Z^{L}\in\Gamma$. We write $T:=(Z^{1}|\dots|Z^{L})^{T}\in\mathbb{R}^{L,s},$ and define the sparse indexing matrices $P^{j}\in\mathbb{R}^{s+1,L}$ in such a way that the rows of $T_{j}:=P^{j}T\in\mathbb{R}^{s+1,s}$ are the labels that make up $\Delta_{j}$.

There exist piecewise linear finite elements $\Phi_{k}\colon\Gamma\to\mathbb{R}$, $k=1,\dots,L$ satisfying $\Phi_{k}(t_{l})=1$ if $k=l$, and $\Phi_{k}(t_{l})=0$ otherwise. In particular, the $\Phi_{k}$ form a partition of unity for $\Gamma$, i.e., $\sum_{k}\Phi_{k}(z)=1\text{ for any }z\in\Gamma$. For a function $p\colon\Gamma\to\mathbb{R}$ in the function space spanned by the $\Phi_{k}$, with a slight abuse of notation, we write $p=(p_{1},\dots,p_{L})$, where $p_{k}=p(Z^{k})$ so that $p(z)=\sum_{k}p_{k}\Phi_{k}(z).$

Functional lifting of the discretized absolute Laplacian  Along the lines of classical continuous multilabeling approaches, the absolute Laplacian regularizer is lifted to become the convex hull of the constraint function $\phi:\mathbb{R}^{L}\to\mathbb{R}\cup\{+\infty\}$,

[TABLE]

where $\mu\geq 0$, $\alpha^{l}\in\Delta^{U}_{s+1}$ (for $\Delta^{U}_{s+1}$ the unit simplex) and $1\leq j_{l}\leq M$ for each $l=0,\dots,m$. The parameter $\mu\geq 0$ is enforcing positive homogeneity of $\phi$ which makes sure that the convex conjugate $\phi^{*}$ of $\phi$ is given by the characteristic function $\delta_{\mathcal{K}}$ of a set $\mathcal{K}\subset\mathbb{R}^{L}$. Namely,

[TABLE]

where $t^{l}:=T_{j_{l}}\alpha^{l}$ and $f(t^{l})$ is the evaluation of the piecewise linear function $f$ defined by the coefficients $(f_{1},\dots,f_{L})$ (cf. above). The formulation of $\mathcal{K}$ comes with infinitely many constraints so far.

We now show two propositions which give a meaning to this set of constraints for arbitrary dimensions $s$ of the labeling space and an arbitrary choice of norm in the definition of $\eta=\|\cdot\|$. They extend the component-wise (anisotropic) absolute Laplacian result in [15] to the vector-valued case.

Proposition 2

The set $\mathcal{K}$ can be written as

[TABLE]

Proof

If the piecewise linear function induced by $f\in\mathbb{R}^{L}$ is concave and 1-Lipschitz continuous, then

[TABLE]

Hence, $f\in\mathcal{K}$. On the other hand, if $f\in\mathcal{K}$, then we recover Lipschitz continuity by choosing $t^{l}=t^{1}$, for any $l$ in (30). For concavity, we first prove mid-point concavity. That is, for any $t^{1},t^{2}\in\Gamma$, we have

[TABLE]

or, equivalently, $[f(t^{1})-f(t^{0})]+[f(t^{2})-f(t^{0})]\leq 0$, where $t^{0}=\frac{1}{2}(t^{1}+t^{2})$. This follows from (30) by choosing $t^{0}=\frac{1}{2}(t^{1}+t^{2})$ and $t^{l}=t^{0}$ for $l>2$. With this choice, the right-hand side of the inequality in (30) vanishes and the left-hand side reduces to the desired statement. Now, $f$ is continuous by definition and, for these functions, mid-point concavity is equivalent to concavity. ∎

The following theorem is an extension of [15, Theorem 1] to the vector-valued case and is crucial for numerical performance, as it shows that the constraints in Prop. 2 can be reduced to a finite number:

Proposition 3

The set $\mathcal{K}$ can be expressed using not more than $|\mathcal{E}|$ (nonlinear) constraints, where $\mathcal{E}$ is the set of faces (or edges in the 2D-case) in the triangulation.

Proof

Usually, Lipschitz continuity of a piecewise linear function requires one constraint on each of the simplices in the triangulation, and thus as many constraints as there are gradients. However, together with concavity, it suffices to enforce a gradient constraint on each of the boundary simplices, of which there are fewer than the number of outer faces in the triangulation. This can be seen by considering the one-dimensional case where Lipschitz constraints on the two outermost pieces of a concave function enforce Lipschitz continuity on the whole domain. Concavity of a function $f\colon\Gamma\to\mathbb{R}$ expressed in the basis $(\Phi_{k})$ is equivalent to its gradient being monotonously decreasing across the common boundary between any neighboring simplices. Together, we need one gradient constraint for each inner, and at most one for each outer face in the triangulation. ∎

2.3 Numerical aspects

For the numerical experiments, we restrict to the special case of integrands $f(x,z,p)=\rho(x,z)+\eta(p)$ as motivated in Section 2.1.

Discretization  We base our discretization on the setting in Section 2.2. For a function $p\colon\Gamma\to\mathbb{R}$ in the function space spanned by the $\Phi_{k}$, we note that

[TABLE]

where $A^{j}$ and $b^{j}$ are such that $\alpha=A^{j}z-b^{j}\in\Delta^{U}_{s+1}$ contains the barycentric coordinates of $z$ with respect to $\Delta_{j}$. More precisely, for $\bar{T}^{j}:=(P^{j}T|-e)^{-1}\in\mathbb{R}^{s+1,s+1}$ with $e=(1,\dots,1)\in\mathbb{R}^{s+1}$, we set

[TABLE]

The functions $\mathbf{u}\colon\Omega\to\mathcal{P}(\Gamma)$ are discretized as $u^{ik}:=\int_{\Gamma}\Phi_{k}(z)d\mathbf{u}_{X^{i}}(z)$, hence $u\in\mathbb{R}^{N,L}$. Furthermore, whenever $\mathbf{u}_{x}=\delta_{u(x)}$, the discretization $u^{i}$ contains the barycentric coordinates of $u(X^{i})$ relative to $\Delta_{j}$. In the context of first-order models, this property is described as sublabel-accuracy in [12, 17].

Dual admissibility constraints  The admissible set $X$ of dual variables is realized by discretizing the conditions (12) and (13).

Concavity (12) of a function $p\colon\Gamma\to\mathbb{R}$ expressed in the basis $(\Phi_{k})$ is equivalent to its gradient being monotonously decreasing across the common boundary between any neighboring simplices. This amounts to

[TABLE]

where $g^{j_{1}},g^{j_{2}}$ are the (piecewise constant) gradients $\nabla p(z)$ on two neighboring simplices $\Delta_{j_{1}},\Delta_{j_{2}}$, and $n_{j_{1},j_{2}}\in\mathbb{R}^{s}$ is the normal of their common boundary pointing from $\Delta_{j_{1}}$ to $\Delta_{j_{2}}$.

The inequality (13) is discretized using (35) similar to the one-dimensional setting presented in [17]. We denote the dependence of $p$ and $q$ on $X^{i}\in\Omega$ by a superscript $i$ as in $q^{i}$ and $p^{i}$. Then, for any $j=1,\dots,M$, we require

[TABLE]

which, for $\rho_{j}:=\rho+\delta_{\Delta_{j}}$, can be formulated equivalently as

[TABLE]

The fully discretized problem can be expressed in convex-concave saddle point form to which we apply the primal-dual hybrid gradient (PDHG) algorithm [4] with adaptive step sizes from [10]. The epigraph projections for $\rho_{j}^{*}$ and $\eta$ are implemented along the lines of [18] and [19].

3 Numerical results

We implemented the proposed model in Python 3 with NumPy and PyCUDA. The examples were computed on an Intel Core i7 4.00 GHz with 16 GB of memory and an NVIDIA GeForce GTX 1080 Ti with 12 GB of dedicated video memory. The iteration was stopped when the Euclidean norms of the primal and dual residuals [10] fell below $10^{-6}\cdot\sqrt{n}$ where $n$ is the respective number of variables.

Image registration  We show that the proposed model can be applied to two-dimensional image registration problems (Figure 1 and 2). We used the sum of squared distances (SSD) data term $\rho(x,z):=\frac{1}{2}\|R(x)-T(x+z)\|_{2}^{2}$ and squared Laplacian (curvature) regularization $\eta(p):=\frac{1}{2}\|\cdot\|^{2}$. The image values $T(x+z)$ were calculated using bilinear interpolation with Neumann boundary conditions. After minimizing the lifted functional, we projected the solution by taking averages over $\Gamma$ in each image pixel.

In the first experiment (Figure 1), the reference image $R$ was synthesized by numerically rotating the template $T$ by $40$ degrees. The grid plot of the computed deformation as well as the deformed template are visually very close to the rigid ground-truth deformation (a rotation by 40 degrees). Note that the method obtains almost pixel-accurate results although the range $\Gamma$ of the deformation is discretized on a disk around the origin, triangulated using only $25$ vertices, which is far less than the image resolution.

The second experiment (Figure 2) consists of two coronal slices from a DCE-MRI dataset of a human kidney (data courtesy of Jarle Rørvik, Haukeland University Hospital Bergen, Norway; taken from [3]). The deformation computed using our proposed model is able to significantly reduce the misfit in liver and kidney in the left half while accurately fixing the spinal cord at the right edge.

Projecting the lifted solution  In the scalar-valued case with first-order regularization, the minimizers of the calibration-based lifting can be projected to minimizers of the original problem [19, Theorem 3.1]. In our notation, the thresholding technique used there corresponds to mapping $\mathbf{u}$ to

[TABLE]

which is (provably) a global minimizer of the original problem for any $s\in[0,1)$.

To investigate whether a similar property can hold in our higher-order case, we applied our model with Laplacian regularization $\eta(p)=\frac{1}{2}\|p\|^{2}$ as well as the calibration method approach with total variation regularization to the data term $\rho(x,z)=(|x|-|z|)^{2}$ with one-dimensional domain $\Omega=[-1,1]$ and scalar data $\Gamma=[-1,1]$ using $20$ regularly-spaced discretization points (Figure 3).

The result from the first-order approach is easily interpretable as a composition of two solutions to the original problem, each of which can be obtained by thresholding (40). In contrast, thresholding applied to the result from the second-order approach yields the two hat functions $v_{1}(x)=|x|$ and $v_{2}(x)=-|x|$, neither of which minimizes the original functional. Instead, the solution turns out to be of the form $\mathbf{u}=\frac{1}{2}\delta_{u_{1}}+\frac{1}{2}\delta_{u_{2}}$, where $u_{1}$ and $u_{2}$ are in fact global minimizers of the original problem: namely, the straight lines $u_{1}(x)=x$ and $u_{2}(x)=-x$.

4 Conclusion

In this work we presented a novel fully-continuous functional lifting approach for non-convex variational problems that involve Laplacian second-order terms and vectorial data, with the aim to ultimately provide sufficient optimality conditions and find global solutions despite the non-convexity. First experiments indicate that the method can produce subpixel-accurate solutions for the non-convex image registration problem. We argued that more involved projection strategies than in the classical calibration approach will be needed for obtaining a good (approximate) solution of the original problem from a solution of the lifted problem. Another interesting direction for future work is the generalization to functionals that involve arbitrary second- or higher-order terms.

Acknowledgments  The authors acknowledge support through DFG grant LE 4064/1-1 “Functional Lifting 2.0: Efficient Convexifications for Imaging and Vision” and NVIDIA Corporation.