Globally Optimal Joint Image Segmentation and Shape Matching Based on Wasserstein Modes

Introduction

Object segmentation and matching are fundamental problems in image processing and computer vision as they form the basis for many high-level approaches to understanding an image. They are intimately related: segmentation of the foreground is a prerequisite for matching in a sequential processing pipeline. Whereas, when performed simultaneously, matching with a template of the sought-after object (e.g. starfish, car, etc.) as prior knowledge can help guiding segmentation to become more robust to corruption of local image features through noise, occlusion and other distortions. Naturally the combined problem is more complicated.

Today convex variational methods can solve image labelling and segmentation problems based on local cues exactly or in good approximation. But combining an object segmentation functional with a shape prior entails a delicate trade-off between descriptive power and computational complexity. Sophisticated shape priors are often described by highly non-convex functionals whereas convex shape prior functionals tend to be rather simplistic. Incompatibility between different shape representations within one approach or the requirement of geometric invariance are common causes of difficulty.

In this paper we present a shape prior functional for simultaneous object segmentation and matching which has been designed specifically to address the issues of representation incompatibility and geometric invariance. Using optimal transport and the differential geometric structure of the 2-Wasserstein space for regularization, one can combine appearance modelling, description of statistical shape variations and geometric invariance in a mathematically uniform way. The linear programming formulation of optimal transport due to Kantorovich allows for an adaptive convex relaxation which can be used for a globally optimal branch and bound scheme, thus avoiding the initialization problem from which most non-convex approaches suffer.

2 Related Literature

Variational methods based on convex relaxations have been successfully applied to obtain globally optimal (approximate) solutions to the originally combinatorial image labelling or segmentation problem . The segmentation is usually encoded by (relaxed) indicator functions which allow for simple and convex formulation of local data matching terms and regularizers that encourage local boundary regularity, such as total variation and its generalizations. However, introducing global regularizers, such as shape priors, into such models is difficult. Convex shape priors based on indicator functions are conceivable but tend to be rather simplistic and lack important features such as geometric invariance .

Somewhat complimentary is the representation of shapes by their outline contours. Treated as infinite dimensional manifolds such representations can be used to construct sophisticated shape modelling functionals . But matching contours with local image data usually yields non-convex functionals that can only be optimized locally via gradient descent. Often one has to internally convert the contour to the region representation. Therefore such approaches require a good initialization to yield reasonable results.

Object Registration.

Independent of the segmentation problem, computing meaningful registrations between two fixed objects (e.g. whole images, measures, meshes…) has attracted a lot of attention. Typical applications are shape interpolation, data interpretation and using registrations as a basis for a measure of object similarity. Often one requires invariance of the sought-after registration under isometric transformations of either of the two objects. Major approaches include the framework of diffeomorphic matching and metamorphosis , methods based on physical deformation energies and the metric approach to shape matching . An extension to shapes that in addition to their geometry are equipped with a ‘signal living on the shape’ is presented in .

These methods provide impressive results at the cost of non-convex functionals and high computational complexity. Naïve online combination with object segmentation is thus not possible. In a shape prior based on object matching has been constructed through convex relaxation of the Gromov-Wasserstein distance .

Optimal Transport.

Optimal transport is a popular tool in machine learning and image analysis. It provides a meaningful metric on probability measures by ‘lifting’ a metric from the base space. Thus it is a powerful similarity measure on bag-of-feature representations and other histograms . It is also applied in geometric problems to extract an object registration from the optimal transport plan . However this requires alignment of the objects beforehand. A step towards loosening this constraint is presented in where one optimizes over a suitable class of transformations. The 2-Wasserstein space, induced by optimal transport, exhibits structure akin to a Riemannian manifold . This was exploited in for analysis of spatial variations in observed sets of measures.

3 Contribution and Outline

We present a functional for object segmentation with a shape prior. Motivated by the literature on object registration, we propose to base the prior on matching the foreground proposal to a template object. For this we need to be able to jointly optimize over segmentation and registration. Matching is done via optimal transport and based both on geometry and local appearance information. Foreground and template are represented as metric measure spaces which provides ample flexibility. This encompasses a wide range of spatial data structures (pixels, super-pixels, point clouds, sparse interest points, …) and local appearance features (color, patches, filter responses, …). Inspired by the Riemannian structure of the 2-Wasserstein space is used to model geometric transformations, object-typical deformations and changes in appearance in a uniform way. Hence, the resulting approach is invariant under translation and approximately invariant under rotation and scaling.

It has recently been shown that this way of modelling transformations and deformations is equivalent to modelling based on closed contours but no conversion of shape representation is required during inference. So shape modelling and local appearance matching are performed directly in the same object representation, allowing to combine the local appearance matching of indicator functions with the manifold based shape modelling on contours. Also, explicitly using the conversion during learning greatly simplifies statistical analysis of the training data and avoids difficulties that arise in .

The resulting overall functional is non-convex, but non-convexity is constrained to a low-dimensional variable, making optimization less cumbersome than in typical contour-based approaches or shape matching functionals. Using the linear programming formulation of optimal transport due to Kantorovich, we derive an adaptive convex relaxation and construct a globally optimal branch and bound scheme thereon. Another option is to apply a local alternating optimization scheme. By employing both optimization techniques one after another their respective advantages (no initialization required, speed) can be combined. This allows to construct a ‘coarse’ object localization method and a subsequent more precise segmentation method as different approximate optimization techniques of the very same functional instead of using two different models. Additionally an efficient graph-cut relaxation is discussed.

The paper is organized as follows: In Sect. 2 the mathematical background for the paper is introduced. We touch upon the convex variational framework for image segmentation, optimal transport and its differential geometric aspects and the description of shapes via manifolds of (parametrized) contours. The proposed functional is successively developed throughout Sect. 3. We start in Sect. 3.1 with a basic segmentation functional where optimal transport w.r.t. a reference template is used as a shape prior. This functional has obvious limitations (e.g. lack of geometric invariance). An alleviation is proposed in Sect. 3.2 by introducing additional degrees of freedoms that allow transformation of the template set. These transformations can be used to achieve geometric invariance and to model statistical object variation, learned from training data (Sects. 3.3 and 3.4). In Sect. 4 we discuss two different approaches for optimization: locally, based on alternating descending steps and globally by branch and bound with adaptive convex relaxations (Sects. 4.1 and 4.2). A relaxation that replaces optimal transport by graph cuts for reduced computational cost is derived in 4.3. Numerical experiments are presented in Sect. 5 to illustrate the different features of the approach and to compare the two optimization schemes. A brief conclusion is given at the end.

4 Notation

Mathematical Background

The first term is referred to as data term, the second as regularizer. The data term s\big{(}y,u(y)\big{)} describes how well label $u(y)$ matches pixel $y$, based on local appearance information. The regularizer $R$ introduces prior knowledge to increase robustness to noisy appearance. A common assumption is that boundaries between objects are smooth, a suitable regularizer then is the total variation.

To obtain feasible convex problems the constraint that $u$ must be binary is usually relaxed to the interval $$ and the functional (2.1) is suitably extended onto non-binary functions, such that it is convex. In the case of total variation regularization such an extension may be

where the data term of (2.1) can be equivalently expressed as a linear function in $u$.

Total variation is a local regularizer in the sense that it only depends locally on the (distributional) derivative of its argument. It can thus only account for local noise, i.e. noise that is statistically independent at different points of the image. Although this weakness can be alleviated to some extent by employing non-local total variation , the inherent underlying assumption is often not satisfied: faulty observations caused by illumination changes or occlusion clearly have long range correlations. At the same time, in particular for the problem of object segmentation more detailed prior knowledge might be available that is not exploited by local regularizers: the shape of the sought-after object. A non-local regularizer that encourages the foreground region to have a particular shape is called a shape prior.

In this article we construct a shape prior by regularization of the foreground region with optimal transport. Hence, we interpret $u$ as the density of a measure $\nu$ w.r.t. the Lebesgue measure $\mathcal{L}_{Y}$ on $Y$. The feasible set for $\nu$ will be:

The first constraint ensures that $\nu\in\textnormal{SegMeas}(Y,M)$ has a density which is a relaxed indicator function. The second constraint fixes the overall mass of $\nu$ to $M$. This is necessary to make it comparable by optimal transport.

2 Optimal Transport

is referred to as the set of couplings between $\mu$ and $\nu$. It is the set of non-negative measures on $X\times Y$ with marginals $\mu$ and $\nu$ respectively.

and the Riemannian inner product for two tangent vectors is given by the $L^{2}$ inner product w.r.t. $\mu$:

Analogous to (2.7) first order variations of a measure $\mu$ along a given tangent vector $t$ are described by

The Jacobian determinant of $T_{\lambda}=\operatorname{id}+\lambda\cdot t$ is

And by the change of variables formula the density of ${T_{\lambda}}_{\sharp}\mu$ is given by

Clearly the concept of optimal transport generalizes to non-negative measures of any (finite) mass, as long as the mass of all involved measures is fixed to be identical. An extensive introduction to optimal transport and the structure of Wasserstein spaces is given in . A nice review of the Riemannian viewpoint can be found in and is further investigated in for sufficiently regular measures.

3 Contour Manifolds and Shape Measures

The tangent space $T_{c}\textnormal{Emb}$ at a given curve $c\in\textnormal{Emb}$ is represented by smooth vector fields on $S^{1}$, indicating first order deformation:

This linear structure is a useful basis for analysis of shapes, represented by closed simple contours, and construction of shape priors thereon (see Sect. 1.2).

Let Diff denote the set of smooth automorphisms on $S^{1}$. In shape analysis one naturally wants to identify different parametrizations of the same curve. This can be achieved by resorting to the quotient manifold $B=\textnormal{Emb}/\textnormal{Diff}$ of equivalence classes of curves, equivalence $c_{1}\sim c_{2}$ between $c_{1},c_{2}\in\textnormal{Emb}$ given if there exists a $\varphi\in\textnormal{Diff}$ such that $c_{1}=c_{2}\circ\varphi$. We write $[c]$ for the class of all curves equivalent to $c$.

One finds that for some $a\in T_{c}\textnormal{Emb}$ the component which is locally tangent to the contour corresponds to a first order change in parametrization of $c$. ‘Actual’ changes of the shape can always be represented by scalar functions on $S^{1}$ that describe deformations which are locally normal to the contour:

where $H$ indicates that this belongs to the horizontal bundle on Emb w.r.t. the quotient $B$. For smooth paths in Emb one can always find an equivalent path such that the tangents lie in $H_{c}\textnormal{Emb}$. While splitting off reparametrization is very elegant from a mathematical perspective, it remains a computational challenge when handling parametrized curves numerically (see for example ).

This maps $a$ to a uniquely determined $t$. We denote this map by $f_{c}$ (depending on the basis contour $c$) and write $t=f_{c}(a)$.

Note that $t=f_{c}(a)$ has constant divergence on $\Omega(c)=\operatorname{spt}F(c)$. Hence by virtue of (2.12) one finds to first order of $\lambda$ that $\mu(\lambda)=(\operatorname{id}+\lambda\cdot t)_{\sharp}F(c)$ has constant density on its support and is therefore itself a shape measure.

This means that describing shapes via shape measures and appropriate tangent vectors thereon is mathematically equivalent to describing shapes by contours modulo parametrization and deformations. Thus we can construct shape priors for regularization with optimal transport, based on measures, without any representation conversion during inference and without having to handle parametrization ambiguities numerically.

Regularization with Optimal Transport

Note that this is conceptually different from matching approaches where a certain local image feature (usually intensity or gray-level) is directly converted into a density. The limitations of this are discussed in in the context of ‘colored currents’. In brief, one problem is, for example, that only one dimensional features can be described. Another is, that, by converting features to density, different, a priori equally important image regions, are assigned different densities and thus have a different influence on the optimizer.

We use the measure to indicate the location of the sought-after object. Local image data is handled in a different fashion: for this we introduce a suitable feature space $\mathcal{F}$. Depending on the image this may be the corresponding color space. It may however also be a more elaborate space spanned by small image patches or local filter responses. We then assume that any point $y\in Y$ is equipped with some $f_{y}\in\mathcal{F}$ which we refer to as the observed feature. We can thus consider every pixel to be a point in the enhanced space $Y\times\mathcal{F}$ with coordinates $(y,f_{y})$.

For regularization with optimal transport we need to provide a prototype, referred to as template. Let $X$ be a set whose geometry will model the shape of the object of interest. It will be equipped with a measure $\mu$ which should usually be the Lebesgue measure on $X$, having density $1$ everywhere, to indicate that ‘all of $X$ is part of the object’. The constant $M$ specifying the total mass for feasible segmentations $\nu$ will be the mass of $\mu$:

Additionally, we describe the appearance of the template by associating to all elements $x\in X$ corresponding $f_{x}\in\mathcal{F}$, the expected features.

Combining this, we can construct a functional for rating the plausibility of a segmentation proposal $\nu\in\textnormal{SegMeas}(Y,M)$:

The first term is the minimal matching cost between the segmentation region and the template via optimal transport with a cost function that combines the geometry and appearance. The second term can contain other typical components of a segmentation functional, for example a local boundary regularizer (cf. Sect. 2.1). The functional is illustrated in Fig. 1a.

Any non-isometric deformation between template foreground and the object will be uniformly penalized by the geometric part of the corresponding optimal transport cost. No information on more or less common deformations (learned from a set of training samples) can be encoded.

Since the mass $M$ of $\mu$, related to the size of the template $X$, equals the mass of $\nu$, this determines the size of the foreground object in $Y$. Hence, the presented functionals imply that one must know the scale of the sought-after object beforehand. This is not possible in all applications.

2 Wasserstein Modes

The function $F$ can be used to introduce statistical knowledge on the distribution of the coefficients $\lambda$. The enhanced functional is illustrated in Fig. 1b.

Functional (3.5) is generally non-convex. For fixed $\lambda$ it is convex in $\nu$. For fixed $\nu$ and a fixed coupling $\pi$ in the optimal transport term it is convex in $\lambda$ if transformations are of the form (3.4) and $F$ is convex. Joint non-convexity does not come as a surprise. It is in fact easy to see that a meaningful isometry invariant segmentation functional with explicitly modelled transformations is bound to be non-convex (Fig. 2).

For optimization of (3.5) assume we first eliminate the high-dimensional variable $\nu$ through minimization (which is a convex problem). One is then left with:

This is in general non-convex, but the dimensionality of $\lambda$ is typically very low (of the order of 10). We can thus still hope to find globally optimal solutions by means of non-convex optimization. We will present a corresponding branch and bound scheme in Sect. 4.2.

When the feature space $\mathcal{F}$ has an appropriate linear structure a natural generalization of (3.4) is to not only model geometric transformations of $X$ but also of the expected features $f_{x}$. In analogy to (3.4) consider

This will be useful when the appearance of the object is known to be subject to variations or when a feature is affected by geometric transformations: for example the expected response to an oriented local filter will need to be changed when the object is rotated. The corresponding generalized functional is

We will further study this generalization in Sect. 5. Meanwhile, for the sake of simplicity we constrain ourselves to purely geometric modes.

3 Geometric Invariance

The framework provided by transformations (3.4) and functional (3.5) allows to introduce geometric invariance into the segmentation / matching approach. In this section we will consider translations, (approximate) rotations and scale transformations. Scale transformations will play a special role as they change the mass of the template.

the corresponding coefficients $\lambda_{\textnormal{t}1},\lambda_{\textnormal{t}2}$ parametrize translations of the template. Further, let $R(\phi)$ be the 2-dimensional rotation matrix by angle $\phi$. Then the mode

Note that $t_{\textnormal{t}1},t_{\textnormal{t}2}$ and $t_{\textnormal{r}}$ have zero divergence. Hence, to first order the implied transformations do not alter the density of $\mu$. For explicit invariance under translations and rotations the modelling function $F$ in (3.5) should be constant w.r.t. the coefficients $\lambda_{\textnormal{t}1},\lambda_{\textnormal{t}2}$ and $\lambda_{\textnormal{r}}$.

Scale.

The size of $X$ and $\mu$ determines the size of the object within the image. In many applications the scale is not known beforehand, thus dynamical resizing of the template during the search is desirable. With slight extensions the framework of transformations can be employed to introduce as a scale-mode into the approach. Let

By the change of variable formula (cf. (2.11,2.12)) the density of ${T_{\lambda}}_{\sharp}\,\mu$ is given by

Thus, introducing a scale mode into (3.5) yields

where we have scaled $\mu$ by the appropriate factor in the feasible set for $\pi$ and we have normalized the first term by a factor of $(1+\lambda_{\textnormal{s}})^{-2}$ to make the term scale invariant. Depending on whether scale invariance is desired the terms $F(\lambda)$ and $G(\nu)$ may need to be rescaled appropriately, too. The feasible set for $\nu$ in $E_{\textnormal{s}}$ is \textnormal{SegMeas}\big{(}Y,(1+\lambda_{\textnormal{s}})^{2}\cdot M\big{)}.

While the modes for translation and rotation leave the area of the template unaltered, statistical deformation modes that we learn from sample data will in general have non-zero divergence. Handling changes in mass will require some extra care during optimization. Therefore we will decompose such modes into a divergence-free part and a contribution of the scale-component.

4 Statistical Variation

One of the limitations of (3.3) discussed in Sec. 3.1 is that non-isometric variations of the template object are uniformly penalized by the geometric component of the corresponding optimal transport cost. However, not all deformations with the same optimal transport cost are equally likely. It may be necessary to reweigh the distance to more accurately model common and less common deformations.

For contour based shape priors a model of statistical object variations is typically learned from samples in a tangent space approximation of the contour manifold. In the tangent space approximation to the Wasserstein space ${\mathcal{W}_{2}}$ was used to analyze typical deformations in a dataset of densities. But mimicking the learning procedure on the contour manifold with optimal transport involves some unsolved problems.

The first problem is to find an appropriate footpoint for the tangent space approximation, i.e. a point by the associated tangent space of which we want to approximate the manifold to first order. One should pick a point which is close to all training samples. Typically one chooses a suitable mean, in a more general metric setting the natural generalization is the Karcher mean. Computation of the barycenter on ${\mathcal{W}_{2}}$ is a non-trivial problem , which has recently been made more accessible through Entropy smoothing . However it becomes more involved when one wants to take geometric invariances into account and impose the constraint of constant density on the support. In the $L^{2}$-mean of the density functions was picked as footpoint after aligning the centers of mass and the principal axes of the samples. Though this is not necessarily an ideal choice (the $L^{2}$-mean of the densities can be very far from some of the samples) it seems to work for smooth densities with limited variations. It will not extend to the binary densities that we consider in this paper since their $L^{2}$-mean need not be binary. In the problem was tentatively solved by manually picking a ‘typical’ sample from the training set as the footpoint.

The second problem is how one maps the samples into the tangent space of the footpoint. A natural choice is the logarithmic map, or some approximation thereof. Recall from Sect. 2.2 that tangent vectors on the manifold of measures are curl-free vector fields and that the logarithmic map is basically obtained by taking the relative transport map. There are some issues with the application to object segmentation: The vector fields computed by the logarithmic map need not have constant divergence, although fluctuations are typically small enough to be ignored for practical purposes. A second issue is that the vector fields are in general not smooth between measures with non-smooth densities, as in our case. This leads to unreasonable interpolations and unwanted artifacts during statistical analysis of the vector fields representing the sample set.

In this paper we circumvent both problems by employing the diffeomorphism between the manifold of contours and the manifold of shape measures (see Sect. 2.3). This allows us to outsource the shape learning problem to the contour representation where established methods for finding a good mean and tangent vectors are available (for example ).

Concretely we used the contour metric and the corresponding approximate algorithmic framework based on gradient descent and dynamic programming presented in for computing the Karcher mean of a set of training shapes and for mapping the training-samples onto the tangent space at the mean via the logarithmic map. We then performed a principal component analysis w.r.t. the Riemannian inner product to extract the dominating modes of shape variation within the training set, together with their observed standard deviation $\{(t_{i},\sigma_{i})\}$. The results we obtained were stable under choosing different initializations. Learning of the class ‘starfish’ is illustrated in Fig. 3. The standard deviations $\sigma_{i}$ were then used to define $F(\lambda)$ to model a Gaussian distribution on the statistical mode parameters:

where $\gamma$ is a parameter determining the weight of $F$ w.r.t. the other functional components.

5 Background Modelling

The previous sections describe how to model the sought-after object via a template, i.e. they focus on the image foreground. Let us now briefly comment on the background.

Sometimes information on the expected appearance of the background is available. This can be incorporated by a linear contribution to $G$ (3.5):

where a positive (negative) coefficient $g(y)$ indicates that a given point is likely to be part of the background (foreground) (cf. Sect. 2.1). Such linear terms can be absorbed into the optimal transport term:

That is, the background appearance model leads to an effective shift of the foreground assignment costs: $c(x,y)\rightarrow c(x,y)+g(y)$.

In other situations it may be desirable to impose that the region directly around the foreground object does not look like foreground itself. An example for such a situation and the corresponding solution are discussed with numerical examples in Sect. 5, see Fig. 5.

Optimization

Functional (3.5) is generally non-convex. It is convex in $\nu$ for fixed $\lambda$ and it is convex in $\lambda$ under suitable conditions (see Sect. 3.2). Based on this, an alternating optimization scheme is conceivable for divergence-free modes. This has also been proposed in [11, Sect. 3.2.1]. We require the following reformulation of (3.6):

Computing (3.6) involves a nested optimization problem over $\nu\in\textnormal{SegMeas}(Y,M)$ and then $\pi\in\Pi(\mu,\nu)$. Given a coupling $\pi\in\Pi(\mu,\nu)$ the marginal $\nu$ can be reconstructed via projection: $\nu={\textnormal{Proj}_{Y}}_{\sharp}\pi$. This allows to reformulate the optimization of (3.6) directly in terms of couplings. Let

and let the feasible set for $\pi$ in $\hat{E}$ be

Then for fixed $\lambda$ one has by construction

and for any optimizer $\pi^{\ast}$ of $\hat{E}$ the marginal ${\textnormal{Proj}_{Y}}_{\sharp}\pi^{\ast}$ is an optimizer of $E$.

Functional $\hat{E}(\lambda,\pi)$ is separately convex in $\lambda$ and $\pi$ for transformations of the form (3.4) and convex $F$. For some initial $\lambda^{1}$ consider the following sequence for $k=1,2,\ldots$:

The sequence of energies $\hat{E}(\lambda^{1},\pi^{1})\rightarrow\hat{E}(\lambda^{2},\pi^{1})\rightarrow\hat{E}(\lambda^{2},\pi^{2})\rightarrow\ldots$ is non-increasing and converges.

Since $\lambda^{k}$ is feasible when determining $\lambda^{k+1}$, one has $\hat{E}(\lambda^{k+1},\pi^{k})\leq\hat{E}(\lambda^{k},\pi^{k})$. Likewise $\pi^{k}$ is a feasible point for computing $\pi^{k+1}$ so $\hat{E}(\lambda^{k+1},\pi^{k+1})\leq\hat{E}(\lambda^{k+1},\pi^{k})$. Hence, the sequence of energies is non-increasing. As $\hat{E}$ is bounded from below, the sequence of energies must converge. ∎

Unfortunately this cannot be extended to modes with non-zero divergence, as changing $\lambda_{\textnormal{s}}$ changes the feasible set \textnormal{SegCoupl}\big{(}Y,(1+\lambda_{\textnormal{s}})^{2}\cdot\mu\big{)} for $\pi$. Thus $\pi^{k}$ need not be feasible for the problem that determines $\pi^{k+1}$ and the sequence of energies created may be increasing. We will provide a workaround for this in the next section (Remark 4.7).

The alternating scheme (4.4) is fast and tends to converge after few iterations. But obviously it need not converge to a global optimum and the result depends on the initialization $\lambda^{1}$. Therefore, similar to contour based segmentation functionals it must be applied with care. In practice application to ‘large’ transformations, e.g. translations and rotations, works only if a good initial guess is available (see Fig. 9). On the other hand it achieves decent results on smaller transformations, as most statistically learned deformations are.

2 Globally Optimal Branch and Bound

For handling large displacement transformations, one needs a global optimization scheme. As discussed in Remark 3.3, for fixed $\lambda$ we can eliminate $\nu$ by a separate convex optimization. One obtains (3.6):

This function is in general non-convex but low dimensional. We thus strive for a non-convex global optimization scheme.

Given Remark 4.1 $E_{1}(\lambda)$ can be written as

If $G$ is zero, then by inserting suitable dummy nodes, computing $E_{1}(\lambda)$ can be written as an optimal transport problem for which efficient solvers are available.

where we have again merged the nested optimizations as above. All occurrences of $\lambda$ are optimized separately and independently over $\Lambda$. By introducing a nested sequence of feasible sets

we obtain an adaptive convex relaxation of $E_{1}(\lambda)$ over $\Lambda$. The relaxation becomes tighter as the set becomes smaller. For application in a branch and bound scheme the following properties are required:

thm:E2 The functional $E_{2}$ has the following properties:

$E_{2}(\Lambda)\leq E_{1}(\lambda)\,\forall\,\lambda\in\Lambda$,

$\lim_{\Lambda\rightarrow\{\lambda_{0}\}}E_{2}(\Lambda)=E_{1}(\lambda_{0})$,

$\Lambda_{1}\subset\Lambda_{2}\Rightarrow E_{2}(\Lambda_{1})\geq E_{2}(\Lambda_{2})$.

Property (i): For any $\lambda\in\Lambda$ obviously

So for any fixed $\pi\in\textnormal{SegCoupl}(Y,\mu)$ (overriding the minimization in (4.6,4.7)) have $E_{2}(\Lambda)\leq E_{1}(\lambda)$. Consequently this inequality will also hold after minimization w.r.t. $\pi$.

For the limit property (ii) note that the functions {c_{\textnormal{geo}}}\big{(}T_{\lambda}(x),y\big{)} and $F(\lambda)$ are continuous functions of $\lambda$. Hence, when $\Lambda\rightarrow\{\lambda_{0}\}$ all involved minimizations will converge towards the respective function values at $\lambda_{0}$ and $E_{2}$ converges as desired.

For the hierarchical bound property (iii) note that for fixed $\pi$ in (4.7) minimization over the larger set $\Lambda_{2}$ will never yield the larger result for all occurrences of $\lambda$. This relation will then also hold after minimization. ∎

With the aid of $E_{2}$ one can then construct a branch and bound scheme for optimization of $E_{1}$. Let

be a finite list of $\lambda$-parameter sets $\Lambda_{i}$ and lower bounds $b_{i}$ on $E_{1}$ on these respective sets. For such a list consider the following refinement procedure: refine(L):

Find the element $(\Lambda_{i^{\ast}},b_{i^{\ast}})\in L$ with the smallest lower bound $b_{i^{\ast}}$.

Let $\textnormal{{subdiv}}(\Lambda_{i^{\ast}})=\{\Lambda_{i^{\ast},j}\}_{j}$ be a subdivision of the set $\Lambda_{i^{\ast}}$ into smaller sets.

Compute $b_{i^{\ast},j}=E_{2}(\Lambda_{i^{\ast},j})$ for all $\Lambda_{i^{\ast},j}\in\textnormal{{subdiv}}(\Lambda_{i^{\ast}})$.

Remove $(\Lambda_{i^{\ast}},b_{i^{\ast}})$ from $L$ and add $\{(\Lambda_{i^{\ast},j},b_{i^{\ast},j})\}_{j}$ for $\Lambda_{i^{\ast},j}\in\textnormal{{subdiv}}(\Lambda_{i^{\ast}})$.

thm:refinement Let $L$ be a list of finite length. Let the subdivision in refine be such that any set will be split into a finite number of smaller sets, and that any two distinct points will eventually be separated by successive subdivision. Set $\textnormal{{subdiv}}(\{\lambda_{0}\})=\{\{\lambda_{0}\}\}$. Then repeated application of refine to the list $L$ will generate an adaptive piecewise constant underestimator of $E_{1}$ throughout the union of the sets $\Lambda$ appearing in $L$. The sequence of smallest lower bounds will converge to the global minimum of $E_{1}$.

Obviously the sequence of smallest lower bounds is non-decreasing and never greater than the minimum of $E_{1}$ throughout the considered region (see \threfthm:E2 (iii) and (i)). So it must converge to a value which is at most this minimum. Assume that $\{\Lambda_{i}\}_{i}$ is a sequence with $\Lambda_{i+1}\in\textnormal{{subdiv}}(\Lambda_{i})$ such that $E_{2}(\Lambda_{i})$ is a subsequence of the smallest lowest bounds of $L$ (there must be such a sequence since $L$ is finite). Since subdiv will eventually separate any two distinct points, this sequence must converge to a singleton $\{\lambda_{0}\}$ and the corresponding subsequence of smallest lowest bounds converges to $E_{2}(\{\lambda_{0}\})=E_{1}(\lambda_{0})$. Since the sequence of smallest lowest bounds converges, and the limit is at most the minimum of $E_{1}$, $E_{1}(\lambda_{0})$ must be the minimum. ∎

When the global optimum is unique, one can see that there also is a subsequence of $\lambda$-sets, converging to the global optimum.

In practice we start with a coarse grid of hypercubes covering the space of reasonable $\lambda$-parameters (e.g. translation throughout the image, rotation within bounds where the approximation is valid and the deformation-coefficients in ranges according to the statistical model) and the respective $E_{2}$-bounds. Any hypercube with the smallest bound will then be subdivided into equally sized smaller hypercubes, leading to an adaptive $2^{n}$-tree cover on the considered parameter range.

The refinement is stopped, when the interval with the lowest bound has edge lengths that correspond to an uncertainty in $T_{\lambda}(x)$ which is in the range of the discretization of $X$ and $Y$. Further refinement would only reveal structure determined by rasterization effects.

The optimum of $E_{1}$ w.r.t. modes that have large displacements (such as translation and rotation) tends to be rather distinct, i.e. there is a small, steep basin around the optimal position. The hierarchical optimization scheme then works rather efficiently.

On the other hand, modes that model smaller, local displacements (e.g. those learned from training samples), often have broad, shallow basins around the optimal value. The branch and bound scheme can then take longer to converge.

Therefore it suggests itself to combine the two optimization schemes: the hierarchical approach is used to determine a good initial guess for translation, rotation and a coarse estimate for the smaller modes. For this the alternating scheme is not applicable due to the non-convexity. But once the broad basin around the global optimum is located, the branch and bound scheme may become inefficient. Conversely, using the estimate of the hierarchical scheme as initialization, we can then expect that the alternating method will give reasonable results.

In the presence of a scale mode one can define $E_{\textnormal{s},1}$ and $E_{\textnormal{s},2}$ equivalent to $E_{1}$ and $E_{2}$ with slight adaptations.

where in the second line we have merged the nested optimization over $\nu$ and $\pi$, see Remark 4.1. To obtain $E_{\textnormal{s},2}(\Lambda)$ all occurrences of $\lambda$ will again be replaced by independent separate optimizations over $\Lambda$. To handle the dependency of the feasible set on $\lambda_{\textnormal{s}}$ consider the following set:

Obviously \textnormal{SegCoupl}\big{(}Y,(1+\lambda_{\textnormal{s}})^{2}\cdot\mu\big{)}\subset\textnormal{SegCoupl}\big{(}Y,(1+\lambda_{\textnormal{s,l}})^{2}\cdot\mu,(1+\lambda_{\textnormal{s,u}})^{2}\cdot\mu\big{)} as long as $\lambda_{\textnormal{s,l}}\leq\lambda_{\textnormal{s}}\leq\lambda_{\textnormal{s,u}}$. Then a possible definition of $E_{\textnormal{s},2}$ equivalent to (4.7) is

where $\lambda_{\textnormal{s,l}}$ and $\lambda_{\textnormal{s,u}}$ are the infimum and supremum of $\lambda_{\textnormal{s}}$ in $\Lambda$. It is easy to see that $E_{\textnormal{s},2\textnormal{a}}$ satisfies \threfthm:E2 w.r.t. $E_{\textnormal{s},1}$. The proof is analogous.

If $G$ is zero the definition of $E_{\textnormal{s},2\textnormal{a}}$ can be improved upon. Consider the following lemma:

thm:Superlinearity For some cost function $c$ and $m>0$ let

Then $f(m_{2})/m_{2}\geq f(m_{1})/m_{1}$ for $m_{2}>m_{1}$.

Assume $f(m_{2})<(m_{2}/m_{1})\cdot f(m_{1})$ for $m_{2}>m_{1}$ and let $\pi^{\ast}_{2}$ be an optimizer for $f(m_{2})$. Then $(m_{1}/m_{2})\cdot\pi^{\ast}_{2}$ is feasible for computation of $f(m_{1})$ and one has

With the aid of \threfthm:Superlinearity one then finds that the following is a suitable variant of $E_{\textnormal{s},2\textnormal{a}}$:

The advantages over $E_{\textnormal{s},2\textnormal{a}}$ are a tighter scaling factor and a simpler feasible set for the optimal transport term.

The alternating optimization scheme presented in Sect. 4.1 only works with zero-divergence modes. The hierarchical optimization scheme can be used to extend this to the scale mode. The non-scale coefficients are determined by separate optimization as before, see (4.4b). The new coefficient $\lambda_{\textnormal{s}}^{k+1}$ and $\pi^{k+1}$ are jointly determined by global hierarchical optimization, while keeping the other mode coefficients fixed (this replaces (4.4a)). This hierarchical scheme will only go over one degree of freedom and thus be very quick. Again one finds a non-increasing sequence that must eventually converge.

3 Graph Cut Relaxation

Both alternating and hierarchical optimization require solving a lot of optimal transport problems. Even with efficient solvers this will quickly become computationally expensive as the size of $X$ and $Y$ or the number of modes increases. If $G$ is non-zero then usually even more so because dedicated optimal transport solvers can no longer be applied directly to compute $E_{1}(\lambda)$. Therefore, in this section we present a mass-constraint relaxation that, for suitable choice of $G$, turns computation of $E_{1}(\lambda)$ into a min-cut problem. This can be solved very fast with dedicated algorithms and therefore the relaxation yields a huge speed-up.

Throughout this section let $X$ and $Y$ be discrete sets, e.g. pixels or super-pixels. The Lebesgue measure on $Y$ is approximated by

for subsets $\sigma\subset Y$, where $m_{y}$ is the area of super-pixel $y$. Any $\nu\in\textnormal{SegMeas}(Y,M)$ can then be expressed as

for all $\sigma\subset Y$ with some function $u_{\nu}\,\colon\,Y\rightarrow$. Let $G$ be a total-variation-like local boundary regularizer of $\nu$, expressed in terms of $u_{\nu}$:

where $\mathcal{G}$ is the set of super-pixel neighbours and $a_{y,y^{\prime}}$ is a weight that models the likelihood of a boundary between neighbours $y$ and $y^{\prime}$. Such weights can be constructed from feature dissimilarity in $y,y^{\prime}$, from the response of edge detectors and from the length of the boundary.

We now relax the template-marginal constraint from the coupling set $\Pi(\mu,\nu)$ and allow $\nu$ to have arbitrary mass. So the feasible set of $\nu$ will be

This is (2.3) without the mass constraint. The ‘couplings’ $\pi$ will be taken from the set

Merging optimizations (see Remark 4.1) yields the feasible set

The relaxed equivalent of $E_{1}$ (4.6) that we consider in this section is

Let $\pi^{\ast}$ be an optimizer of $E_{\textnormal{r},1}(\lambda)$ for some configuration $\lambda$. If $({\textnormal{Proj}_{Y}}_{\sharp}\pi^{\ast})(y)>0$ for some $y\in Y$, this mass will come from the cheapest $x\in X$ for this $y$, since there is no longer any constraint on the mass on $X$. The linear matching in the first term simplifies to a nearest neighbour matching for each $y\in Y$. This implies that the minimization in (4.23) over $\pi\in\textnormal{SegCoupl}(Y)$ can be simplified to a minimization over $\nu\in\textnormal{SegMeas}(Y)$. Therefore (4.23) is equivalent to

We express now $\nu$ in terms of $u_{\nu}$, see (4.18), and plug in the form of the regularizer $G$ (4.19). This yields

For fixed $\lambda$ this is a convex formulation of the max-flow / min-cut problem with nodes $Y$ and edges $\mathcal{G}$. The edge-weight between $y\in Y$ and the sink is given by ${c_{\textnormal{min}}}(y,\lambda)\cdot m_{y}$ and the weights of the edges between $y,y^{\prime}\in Y$ by $a_{y,y^{\prime}}$. This problem can be solved very efficiently by dedicated algorithms, see for example .

Both the alternating method and the hierarchical scheme, Sects. 4.1 and 4.2, can be applied directly to the optimization of $E_{\textnormal{r},1}$. The sequence equivalent to (4.4) will provide a non-increasing converging sequence of energies. Since the dependence of the feasible set on the mass of $\mu$ has disappeared, it can also be extended to the scale mode. Also, handling the scale mode in the hierarchical scheme is simplified.

Functional (4.26) can be interpreted as a binary Markov random field (MRF) with labels fore- and background ($u\in\{1,0\}$) and a latent object configuration variable $\lambda$. Such enhanced MRFs have been used in with the latent variables describing layered pictorial structures and in with graph-based shape models. Optimization of a general class of such models via branch and bound has been discussed in . A main difference of the approach presented here and is that the shape variations are not captured implicitly in the hierarchical cluster of sample shapes but explicitly and smoothly in the set of learned Wasserstein modes.

Numerical Examples

We will now present some numerical examples for joint image segmentation and shape matching with Wasserstein modes. The scope of these examples is to transparently show the key properties of the functional (geometric invariance, response to noisy data etc.) and to demonstrate its applicability to different types of geometric data and features.

As discussed in Sect. 3.3 the functional component $F$, modelling the distribution of the deformation parameter $\lambda$ (c.f. (3.5)), was not depending on the $\lambda$-entries that describe translation, rotation and scale. For the statistical modes we modelled a simple Gaussian as given by (3.16). The weight $\gamma$ was set to a small value, i.e. we ‘trusted’ the data for small deformations and mainly wanted to keep the deformations from becoming too large, where the linear deformation model does no longer work very well.

The number of used modes ranged between 3 and 8 for branch and bound, up to about 14 for the alternating scheme. As discussed in Sect. 4.2, for the initial covering $L$ of the parameter space for $\lambda$, we used a grid of $n$-dimensional hypercubes: for the translation components ranging over the area of the image, for rotation and scale within the limits where the numerical approximation is valid and for the statistical modes depending on the observed standard deviations during learning.

A very important parameter in the functional is the relative weight between the geometric and the appearance cost function, ${c_{\textnormal{geo}}}$ and $c_{\mathcal{F}}$. When the appearance features are very noisy, we tend to put more trust on the geometric component and thus the predefined deformation modes. For very reliable data we may accept a previously unknown deformation to better match the observed features. Some intuition on how to choose this relative weight may be gained from Fig. 8.

Optimization Algorithms.

In most experiments numerical optimization was carried out in two steps, starting with branch and bound over the modes with largest deformations, followed by alternating optimization over all modes (see Remark 4.5). As pointed out in Remark 3.2 continuous solvers cannot be applied since the marginal $\nu$ is unknown. Therefore we rely on discrete algorithms. For numerical optimization of $E_{2}(\Lambda)$ we implemented two different methods:

For $G=0$, i.e. in the absence of an additional segmentation term on the marginal $\nu$ (c.f. (3.5)), the functional $E_{2}(\Lambda)$ (4.7) can be evaluated by using a dedicated optimal transport solver. For this we wrote a c++ implementation of the Hungarian method .

When $G$ is a discrete total-variation-like local regularity prior (see Sect. 4.3, eq. (4.19)) evaluation of $E_{2}(\Lambda)$ can be written as a linear program, which we solved with CPLEX.

The top level, i.e. everything except for the calls to optimize $E_{2}(\Lambda)$ was implemented in Mathematica.

In practice we used the first variant for the branch and bound stage and the second variant for the subsequent alternating optimization stage. The reasoning behind this is that the total variation of a segmentation depends mostly on its local properties and can be vary significantly without altering its global configuration, which is what we look for during the branch and bound optimization. TV is then added during the ‘fine-tuning’ in the alternating stage.

Reducing Complexity in Practice.

To reduce computational complexity, we sampled the cost function ${c_{\textnormal{geo}}}(x,y)+c_{\mathcal{F}}(f_{x},f_{y})$ for fixed $x$ only at positions $y$ close to $x$. When $y$ is very far from $x$ the high geometric cost will make the assignment very unlikely. The cut-off radius around $x$ is chosen according to the range of $c_{\mathcal{F}}$ and the size of the mode parameter set $\Lambda$ during branch and bound. Global optimality of the sub-sampled cost-function w.r.t. the dense model can be checked by introducing ‘overflow’ variables with suitable assignment costs for each $x$: as long as no mass is put onto these overflow variables, the optimizer of the reduced model is also globally optimal in the dense model.

Computational Complexity and Runtime.

Although we only used experimental code, which was far from being optimized for performance we briefly comment on the observed running-times to give the reader a general idea of the applicability. Experiments were performed on a standard desktop computer with an Intel Core i7 processor at 3.4 GHz and 16 GB RAM. The branch and bound scheme, which is the computationally most demanding part, was parallelized over the processor cores. The alternating optimization is much less demanding and consequently converges much faster. The discrete templates had several 100 points, the discrete images, super-pixel segmentations, etc. several 1000 points.

For the branch and bound scheme the running time is determined by how many of the tree of bounds at different scales have to be explored until a minimizer is found. This number is sensitive to several factors: it grows exponentially with the number of degrees of freedom. Also, it depends on the specific problem instance and how well the global optimum is pronounced. In the presence of strong noise or multiple similarly good minima the scheme will naturally take longer as in a problem with only one distinct solution. Consequently it is not really possible to accurately estimate the number of required bounds beforehand, i.e. to give an overall expected complexity estimate of the branch and bound scheme.

During our experiments we observed running times from under a minute for 3-4 modes on ‘easy problems’ up to about a day for 7-8 modes on very noisy and large instances. Instances shown in this section were mostly set up such that branch and bound would take 10 minutes at most.

Of course the running time also depends strongly on the problem dimensions. Fortunately, the flexible mathematical framework provides means for reducing the problem dimensions easily by working for example on an over-segmentation with super-pixels instead of on the full pixel grid. The loss of resolution can often be compensated for by adding a local regularizer.

2 Numerical Results

We start with some synthetic experiments to transparently illustrate different properties of the functional. For these experiments the feature cost function $c_{\mathcal{F}}(f_{x},f_{y})$ was chosen to be constant w.r.t. $x$, i.e. every template point expects the same features and the template has a homogeneous appearance. This corresponds to a classifier that tries to locally asses for each pixel whether it is part of the fore- or background.

A shape model of a bunny is learned from several different views. The subsequent task is then to find a novel view (within the range of the training views) among a collection of different shapes. Branch and bound was used to optimize over translations, rotation and scale of the object. On these degrees of freedom the alternating scheme is prone to getting stuck in a poor local minimum, if initialized on the wrong shape. Afterwards the alternating scheme was applied to account for non-isometric variations due to perspective. Additionally it is shown how the rotation invariance can be extended to large angles. The results of this experiment are illustrated in Fig. 4.

Background Modelling.

Note that in order to locate the bunny correctly, we sometimes also need to model the image background in some way. This can be done implicitly by ensuring that the boundary of the foreground is aligned with detected contours in the image via a weighted TV-like term through $G(\nu)$ in (3.5). A more explicit approach is to extend the template to include a small region ‘looking like background’ around the foreground (see Sect. 3.5). This is demonstrated in Fig. 5.

More examples on detecting objects in a noisy environment and on restoring shapes from distorted detections are given in Figs. 6 and 7.

Interaction of Regularizers.

Now let us study the interaction between the different components of the functional. Let $G$ be the discrete total variation of $\nu$ (4.19). For now we ignore deformations and simply take a fixed template. That is we consider the following functional:

where we have introduced weights $\tau$ and $\sigma$. In Fig. 8 it is illustrated how the optimal segmentations depend on $\tau$ and $\sigma$ in the presence of different types of noise.

Alternating Optimization.

In Fig. 9 the behaviour of the alternating optimization scheme is elucidated. In particular it becomes apparent how in noisy problems the scheme easily gets stuck in poor local minima. This is a general problem of local optimization methods and proves the importance of the globally optimal branch and bound scheme to provide a proper initial starting point.

Super-pixels.

An important feature of functional (3.5) is that its discrete version readily encompasses a wide range of data structures. As the computational complexity strongly depends on the size of the discretizations of $X$ and $Y$ it may be reasonable to apply the functional not directly to the pixel level but to a coarser over-segmentation as for example provided by super-pixels. Some examples with the class ‘starfish’ are given in Fig. 10. Fig. 11 shows some of the involved non-isometric deformations to illustrate the range of the linear modes model and also one example where the limit of the linear expansion has been reached.

In Fig. 12 the scale invariance of the approach is demonstrated by actually deliberately breaking it. The same functional is optimized twice, but with a different prior on the allowed object scale. Depending on the admissible scale, once the large and once the small clownfish is segmented. Such a task can only be solved with global optimization techniques.

So far we have only considered the case where $c_{\mathcal{F}}(f_{x},f_{y})$ was constant w.r.t. $x$. However, computationally there is no increase in complexity if we pick a more general feature cost. The potential of this additional freedom is now demonstrated on an example with the UIUC database (see for example ). This is a set of gray level side views of parking cars. Locating these cars cannot be approached with a homogeneous foreground / background detector, as no consistent separation based on local appearance features seems to be possible.

Therefore we now learn local detectors for each point of the template $X$ separately and based on these compute an inhomogeneous $c_{\mathcal{F}}$. As features we use local histograms of the image color and its gradient. We compute assignments between the learned template and the training cars (both shapes fixed, only geometric, no appearance cost). Based on these assignments we extract for each template point $x$ the collection of expected features $f_{x}$. Then, on a test image $Y$ we compare for each super-pixel $y\in Y$ its histogram of features $f_{y}$ with the distribution of expected features $f_{x}$ on each template point via an optimal transport based histogram distance (see e.g. ). These comparison costs were used as costs $c_{\mathcal{F}}(f_{x},f_{y})$.

We want to emphasize at this point that we do in no way champion this particular choice of features and this choice does not constitute a part of our presented framework. We merely seek to provide a transparent set-up to demonstrate the benefit of locally adaptive template appearance without obstruction through more complicated feature acquisition and processing.

Fig. 13 gives an impression of the functions $c_{\mathcal{F}}(f_{x},f_{y})$ obtained in this way. Obviously, for a single template point $x\in X$ the associated cost is very noisy and not very informative. We can thus only hope that through the combination of all template pixels and the knowledge about their relative spatial arrangement we can identify the positions of the cars.

Since the variation of the shapes of the cars is small we only consider translations during branch and bound for locating the cars. Geometric flexibility beyond that is provided by the optimal transport matching. In this way on 10 out of 15 test images the global optimum correctly corresponded to a car (some images show multiple cars). As baseline we performed a simple Hough transform which failed to correctly locate any car. Fig. 14 gives some example cases and also illustrates a failed case.

A similar experiment was performed in . There the main focus was on modelling the boundary of the cars whereas here we concentrate on its region. Both approaches can incorporate both cues from the object interior as well as its boundary. In Sect. 4.3 it was discussed how is closely related to the graph-cut relaxation of our functional. The most significant difference is how in our approach the geometric variability is explicitly modelled by a linear space of modes whereas in it is implicitly encoded in a hierarchical clustering.

We have already mentioned in Remark 3.4 that the Wasserstein modes can also be extended beyond geometric variations to the feature component. This is of particular use when an expected feature is known to change under a certain geometric transformation. For example the orientation of an expected gradient changes with rotation. More generally, a vector valued feature $f_{x}$ will have to be transformed by

the Jacobian of the applied transformation, to preserve it’s ‘relative orientation’ within the template, and we see that this yields a linear deformation on the feature space.

Here we provide a simple example to point out the potential of this flexibility. We now assume that both location and expected feature of a template point vary with the transformations. We model this by linearly expanding $c_{\mathcal{F}}$ in $\lambda$ around the origin. That is we choose (c.f. (3.7-3.9)):

where $c_{\mathcal{F},i}(f_{x},f_{y})$ is the partial derivative of the feature component of \hat{c}\big{(}\hat{T}_{\lambda}(x),(y,f_{y})\big{)} w.r.t. $\lambda_{i}$ evaluated at zero (thus giving the first order change along the feature component of $\hat{t}_{i}$). Both discussed optimization schemes can easily be adapted to this extension.

As a toy example we will be looking for apples. Unripe, small apples are assumed to be green, ripe, large apples should have a reddish color. That is, the expected color varies with size (Naturally the apparent size of an apple on the image depends strongly on the distance from the camera. But we will generously overlook this for the sake of the demonstration.) The results of our search for fruit are illustrated in Fig. 15.

Point Clouds.

Last but not least we want to further illustrate the flexibility of the numerical framework by applying it to a scenario with point clouds. This is relevant when one does not deal with dense images but only with sparse interest points. We give a transparent, synthetic example in Fig. 16.

Conclusion

We have presented a functional for simultaneous image segmentation and shape matching to correctly locate and segment objects within images under noisy conditions.

Matching is based on optimal transport with a cost function that combines geometric plausibility with consistency of appearance features. Through the convex Kantorovich formulation in terms of coupling measures the functional can naturally be combined with other segmentation terms known from convex variational image segmentation. To implement geometric invariances and to account for non-isometric shape variations we introduced additional degrees of freedom, drawing from the Riemannian structure of the 2-Wasserstein space. Through an equivalence relation of the class of shape measures with closed contours this enabled us to introduce well established shape analysis tools from the contour regime into the segmentation approach while remaining in the measure representation.

While the resulting functional is non-convex, this non-convexity is constrained to a low dimensional variable which allowed us to devise an adaptive convex relaxation on which a globally optimal branch & bound optimization scheme could be constructed. Alternatively, a faster but only locally optimal alternating optimization scheme was discussed. While it seems impractical to run the branch and bound scheme on a high number of deformation modes, it still provides a consistent way to find good initializations for the alternating scheme, thus overcoming a severe problem in many other segmentation / matching approaches. Determining a good initial guess and the subsequent ‘fine tuning’ are based on the very same model and only differ in the application of the optimization scheme. To reduce numerical complexity, a graph-cut relaxation was discussed.

In Sect. 5 we presented a series of numerical examples to demonstrate various aspects of the approach. The basic behaviour of the branch and bound scheme was illustrated as well as the limitations of the alternating scheme. It was shown how the location and shape of the optimal segmentations depend on noise and how different kinds of noise can at least partially be handled by properly choosing the weights between the different terms of the functional. We put a particular focus on illustrating the flexibility in both spatial data structure (pixels, super-pixels, point clouds) as well as in incorporating different types of knowledge on the object appearance (spatially varying, adaptive to deformations).

In the presented state a major limitation of the functional is the linearity of the modes: this makes it difficult to handle large deformations. In this respect other approaches such as the LDDMM framework are already much further developed, yet focus on smooth registration mappings without addressing variational segmentation simultaneously and explicitly. On the other hand we notice that in terms of handling local feature data this approach is similarly flexible (compare for example with ). Also, we consider the branch and bound scheme as an important step towards coherently solving the initialization problem.

Future work should therefore focus on making the deformations more flexible and powerful while trying to retain the ability to obtain robust initializations.

Acknowledgement. This work was supported by the DFG, grant GRK 1653.

References