[1]Fabian Parzer

1]Faculty of Mathematics, University of Vienna, Oskar-Morgenstern-Platz 1,
Vienna, 1090, Austria

2]Johann Radon Institute for Computational and Applied Mathematics (RICAM),
Altenbergerstrasse 69, Linz, 4040, Austria

3]Christian Doppler Laboratory for Mathematical Modeling and Simulation of
Next Generations of Ultrasound Devices (MaMSi), Oskar-Morgenstern-Platz 1,
Vienna, 1090, Austria

# Uncertainty Quantification for Scale-Space Blob Detection

<EMAIL_ADDRESS>Clemens Kirisits<EMAIL_ADDRESS>Otmar Scherzer<EMAIL_ADDRESS>[  [  [

###### Abstract

We consider the problem of blob detection for uncertain images, such as images
that have to be inferred from noisy measurements. Extending recent work
motivated by astronomical applications, we propose an approach that represents
the uncertainty in the position and size of a blob by a region in a three-
dimensional scale space. Motivated by classic tube methods such as the taut-
string algorithm, these regions are obtained from level sets of the minimizer
of a total variation functional within a high-dimensional tube. The resulting
non-smooth optimization problem is challenging to solve, and we compare
various numerical approaches for its solution and relate them to the
literature on constrained total variation denoising. Finally, the proposed
methodology is illustrated on numerical experiments for deconvolution and
models related to astrophysics, where it is demonstrated that it allows to
represent the uncertainty in the detected blobs in a precise and physically
interpretable way.

###### keywords:

uncertainty quantification, blob detection, scale space, total variation
regularization

##  1 Introduction

Blob detection is a classic task in computer vision. Here, we mean by a _blob_
a round structure with a roughly Gaussian intensity profile. In order to
simultaneously estimate the position and size of such blobs, detection methods
often rely on the scale-space representation of an image, which represents the
image at different levels of smoothing, allowing to distinguish low-scale from
high-scale structures. This approach is commonly refered to as scale-space
blob detection. The most well-known example of this is the Laplacian-of-
Gaussians (LoG) method Lindeberg (1998), which is based on the premise that,
for the Gaussian scale-space representation, the local extrema of the
normalized Laplacian are good indicators for the position and size of a blob.

In many applications – in particular in astronomy – the image of interest is
not known a-priori but has to be reconstructed from noisy measurements. This
means that the image comes with significant uncertainties, and it is important
to take these uncertainties into account when performing blob detection.

A particular example from astronomy is integrated-light stellar population
recovery Boecker et al. (2020), where the problem is to detect stellar
populations as blobs in a two-dimensional mass density that is reconstructed
from an observed spectrum. For this problem, the present authors with co-
authors have developed an uncertainty-aware version of the Laplacian-of-
Gaussians blob detection method, ULoG Parzer et al. (2023). ULoG was
formulated as a tube method that detects significant blobs by computing a
tube-shaped confidence region for the uncertain scale-space representation,
and then solving a minimization problem designed to obtain a representative
that exhibits the "least amount of blobs".

While the results of the ULoG method were satisfying for this particular
application, it only yielded a very rudimentary representation of the
uncertainty with respect to the position and size of a blob.

In this paper, we propose an improved method that aims to resolve this issue.
The basic idea is to represent the uncertainty in a blob by a region in scale
space which represents the possible variation in position and size.

To obtain these regions, we formulate an optimization problem that enforces
solutions with piecewise constant normalized Laplacian, from which the desired
blob regions can be extracted as level sets. The formulation of the
optimization problem uses ideas from total variation (TV) regularization,
which is why we refer to the novel method as TV-ULoG.

###  1.1 Contributions

  * •

We introduce the TV-ULoG method for blob detection with uncertainties. The
proposed approach is flexible and can be adapted to a wide range of
applications, in particular Bayesian imaging. We also discuss connections to
the taut-string algorithm and constrained total-variation denoising.

  * •

We extensively discuss the numerical treatment of the resulting non-smooth,
bound-constrained convex optimization problem. We compare approaches based on
smoothing the dual or the primal problem, and an interior-point approach based
on reformulating the optimization problem as SOCP.

  * •

Finally, we illustrate the TV-ULoG method on numerical test cases for
astronomical imaging and deconvolution.

###  1.2 Organization of the paper

The paper is organized as follows:

  * •

In the remainder of this section, we review related work (Section 1.3) and
introduce notation that is used throughout the paper (Section 1.4).

  * •

In Section 2, we recall the necessary prerequisites on scale-space blob
detection. We focus on the Gaussian scale-space representation and the
Laplacian-of-Gaussians blob detection method, which are fundamental for the
rest of the paper.

  * •

In Section 3, we describe our tube-based approach for uncertainty-aware blob
detection. After discussing scale-space aspects of uncertainty quantification
and the ULoG method from our previous work, we derive the novel TV-ULoG
method.

  * •

In Section 4, we discuss in detail the numerical implementation of TV-ULoG.
The majority of the section is devoted to the solution of the resulting
optimization problem.

  * •

In Section 5, we demonstrate the method on two numerical test cases. We also
use these test cases to evaluate the performance of the proposed optimization
methods.

  * •

The paper ends with a conclusion in Section 6.

###  1.3 Related work

We have based our work on the Laplacian-of-Gaussians method for scale-space
blob detection Lindeberg (1998) since it is well-known and its mathematical
formulation is simple, making it easier to extend it to the case of uncertain
images. Alternative methods for blob detection are discussed e.g. in the
reviews Mikolajczyk et al. (2005); Lindeberg (2013); Li et al. (2015). Some
general references on scale-space methods for image processing and computer
vision are Lindeberg (1994); Sporring et al. (1997); Lindeberg (1994);
Weickert (1998).

Our work can be seen as an instance of a statistical scale space method
Holmström and Pasanen (2017), but this is not a focus of this paper. In
particular the works Godtliebsen et al. (2004); Godtliebsen and Øigård (2005);
Thon et al. (2012); Pasanen et al. (2013) are similar since they also study
uncertain signals in scale space. However, our approach differs through its
formulation as a tube method and the specific focus on blob detection.

Another related line of work formulates significance tests for image
structures as convex optimization problems Repetti et al. (2019); Price et al.
(2021, 2020). This methodology relies on concentration inequalities Pereyra
(2017) and is computationally very efficient, but does not automatically
detect the position and scale of structures since the presence of a structure
must be formulated as user-specified hypothesis.

To our knowledge, before Parzer et al. (2023), the specific problem of
uncertainty-aware blob detection has not been addressed previously.

###  1.4 Notation

  * •

For $n\in\mathbb{N}$, we define the discrete range
$[n]:=\left\\{1,\ldots,n\right\\}$.

  * •

$\mathbb{R}_{+}:=[0,\infty)$ denotes the non-negative real numbers

  * •

For a vector ${\bm{x}}\in\mathbb{R}^{n}$, its Euclidean norm is denoted by
$\left\|{\bm{x}}\right\|:=\sqrt{\sum_{i=1}^{n}x_{i}^{2}}$.

  * •

The convolution of two functions is denoted by $f*g({\bm{x}}):=\int
f({\bm{y}})g({\bm{x}}-{\bm{y}})\,\mathrm{d}{\bm{y}}$.

  * •

The spatial Laplace operator is denoted by
$\Delta:=\partial_{x_{1}}^{2}+\partial_{x_{2}}^{2}$.

  * •

The probability distribution of a random element $U$ is denoted by
${\mathbb{P}}_{U}$ Kallenberg (2002). Its corresponding density function is
denoted by $p_{U}$, if it exists. Given another random element $V$, the
conditional distribution of $U$, given $V=v$, is denoted by
${\mathbb{P}}_{U|V}(\cdot|v)$, with corresponding conditional density
$p_{U|V}(\cdot|v)$ (see also Schervish (1995)).

  * •

Given functions $u^{\mathrm{low}},u^{\mathrm{upp}}:D\to\mathbb{R}$ on a domain
$D\subset\mathbb{R}^{n}$, with $u^{\mathrm{low}}\leq u^{\mathrm{upp}}$, we
call the set of functions

$\displaystyle:=\\{$ $\displaystyle u:D\to\mathbb{R}~{}:$ (1) $\displaystyle
u^{\mathrm{low}}({\bm{x}})\leq u({\bm{x}})\leq
u^{\mathrm{upp}}({\bm{x}})~{}\forall{\bm{x}}\in D\\}$

a _tube_. Similarly, given two vectors
${\bm{u}}^{\mathrm{low}},{\bm{u}}^{\mathrm{upp}}\in\mathbb{R}^{n}$, we call
the set of vectors

$\displaystyle:=\\{{\bm{u}}\in\mathbb{R}^{n}~{}:~{}u^{\mathrm{low}}_{i}\leq
u_{i}\leq u^{\mathrm{upp}}_{i}$ (2) $\displaystyle\forall i\in[n]\\}$

a _discrete tube_. This definition is straightforward to extend to higher-
dimensional objects, such as discrete images.

  * •

We denote the characteristic function of a set $C$ by

$\displaystyle\chi_{C}({\bm{v}}):=\begin{cases}0,&\text{if }{\bm{v}}\in C,\\\
\infty,&\text{otherwise}.\end{cases}$

  * •

We let $\bm{0}_{n}\in\mathbb{R}^{n}$ denote the zero-vector in
$\mathbb{R}^{n}$ and $\bm{1}_{n}\in\mathbb{R}^{n}$ denote the vector with all
entries equal to 1. Similarly, $\bm{0}_{m\times n}\in\mathbb{R}^{m\times n}$
denotes the zero matrix and $\bm{1}_{m\times n}\in\mathbb{R}^{m\times n}$
denotes the matrix with all entries equal to 1. Also,
$\bm{e}_{i}^{n}\in\mathbb{R}^{n}$ denotes the $i$-th basis vector in
$\mathbb{R}^{n}$, with entries $(\bm{e}_{i}^{n})_{j}=\delta_{ij}$.

  * •

Given a nonempty closed convex set $C\subset\mathbb{R}^{d}$,
$P_{C}(x):=\operatorname*{argmin}_{{\bm{y}}\in
C}\left\|{\bm{y}}-{\bm{x}}\right\|$ denotes the projection on $C$.

  * •

For a set $S$, $2^{S}$ denotes its power set.

##  2 Scale-space blob detection

In this section, we review the scale-space approach to blob detection that
underlies the rest of this paper. We focus on the Gaussian scale-space
representation, which we introduce in Section 2.1. Then, we review the classic
Laplacian-of-Gaussians method for blob detection in Section 2.2.

###  2.1 Scale-space representations

In the computer vision literature, the mathematical theory of describing
images at different scales is known as _scale space theory_ Lindeberg (1994).
A scale space representation of an image $f:\mathbb{R}^{2}\to\mathbb{R}$ is a
function $u:\mathbb{R}^{2}\times\mathbb{R}_{+}\to\mathbb{R}$ which depends on
an additional third parameter $t$ that represents physical scale. The
_Gaussian scale-space representation_ is the most studied example, due to its
simple formulation and the fact that it is the unique linear scale space
representation that satisfies a series of intuitive axioms that formalize the
notion of scale Alvarez et al. (1993); Florack et al. (1992).

The Gaussian scale-space representation of an image $f:D\to\mathbb{R}$ on a
domain $D\subset\mathbb{R}^{2}$ is defined as the solution
$u:D\times[0,\infty)\to\mathbb{R}$ of the diffusion equation

$\displaystyle\partial_{t}u({\bm{x}},t)$ $\displaystyle=\frac{1}{2}\Delta
u({\bm{x}},t),\quad$ $\displaystyle({\bm{x}},t)\in D\times(0,\infty),$ (3)
$\displaystyle u({\bm{x}},0)$ $\displaystyle=f({\bm{x}}),\quad$
$\displaystyle{\bm{x}}\in D,$ $\displaystyle\partial_{\bm{\nu}}u({\bm{x}},t)$
$\displaystyle=0,\quad$ $\displaystyle({\bm{x}},t)\in\partial
D\times(0,\infty).$

Here, we imposed Neumann boundary conditions, following Duits et al. (2003).
In the following, we will denote with $\Phi$ the solution operator of (3) that
maps an image $f$ to its scale-space representation $u$. One can show that
$\Phi$ is well-defined under suitable assumptions (see e.g. Brezis (2011)).

For the rest of this paper, we will not consider other scale-space
representations. In particular, we will always mean the representation $u$
defined in (3) when we write "the scale-space representation" of an image $f$.

###  2.2 Blob detection

The scale-space representation of an image is often used to detect the
position and size of features of interest. A particular example of this is the
well-known Laplacian-of-Gaussians method for blob detection Lindeberg (1998).
It is a special case of the differential-invariants approach for feature
detection which we introduce next.

####  2.2.1 Feature detection from differential invariants

Image features often correspond well to local extrema of _scale-normalized
derivatives_ of the image’s scale-space representation $u$, i.e. combinations
of the scaled partial derivatives

$\displaystyle\tilde{\partial}_{x_{i}}u({\bm{x}},t):=\sqrt{t}\partial_{x_{i}}u({\bm{x}},t),\quad
i=1,2$ (4)

(cf. Lindeberg (2013)). The scale-normalization is necessary to achieve an
intuitive scale-invariance property: A feature of scale $t$ in the image $f$
should correspond to a feature of scale $s\cdot t$ in the rescaled image
$f_{s}({\bm{x}}):=f({\bm{x}}/\sqrt{s})$. E.g. zooming in by a factor 2
increases the scale of all features by a factor 4.

####  2.2.2 The Laplacians-of-Gaussians method

The Laplacian-of-Gaussians method is a special case of the methodology
described in Section 2.2.1 for blob detection. It uses the local minimizers of
the scale-normalized Laplacian of $u$ as indicators for blobs in an image,
where the scale-normalized Laplacian is given by

$\displaystyle\tilde{\Delta}u({\bm{x}},t):=(\tilde{\partial}_{x_{1}}^{2}+\tilde{\partial}_{x_{2}}^{2})u({\bm{x}},t)=t\Delta
u({\bm{x}},t),$ (5)

That is, a local minimizer or maximizer of $\tilde{\Delta}u$ at
$({\bm{x}},t)\in D\times\mathbb{R}_{+}$ indicates the presence of a blob-like
shape with center ${\bm{x}}$ and scale $t$.

######  Example 1.

To understand what is meant by "blob-like shape", let us define the
prototypical blob with center $\bm{m}\in\mathbb{R}^{2}$ and size $s$ as the
symmetrical Gaussian

$\displaystyle f({\bm{x}}):=\frac{1}{2\pi
s}\exp\left(-\frac{\left\|{\bm{x}}-\bm{m}\right\|^{2}}{2s}\right),\quad{\bm{x}}\in\mathbb{R}^{2}.$

Its scale-space representation is

$\displaystyle
u({\bm{x}},t)=\frac{1}{2\pi(s+t)}\exp\left(-\frac{\left\|{\bm{x}}-\bm{m}\right\|^{2}}{2(s+t)}\right),$

and the scale-normalized Laplacian is

$\displaystyle\tilde{\Delta}u({\bm{x}},t)=t\frac{\left\|{\bm{x}}-\bm{m}\right\|^{2}-2(s+t)}{2(s+t)^{2}}u({\bm{x}},t),$

which has a unique local minimum at $(\bm{m},s)$. This means that the position
and size of the prototypical Gaussian blob $f$ are exactly recovered from the
local minimizer of $\tilde{\Delta}u$. Note that the normalization in (5) is
important for detecting the scale, since the un-normalized Laplacian $\Delta
u$ has a unique local minimum at $(\bm{m},0)$, which does not indicate the
blob size. The relation of the Laplace operator to blob detection is also
discussed in Marr and Hildreth (1980); Voorhees and Poggio (1988).

The prototypical Gaussian blob also motivates the common visualization of the
results of the Laplacian-of-Gaussians method, where a detected blob
$(\bm{x},t)$ is visualized by a circle with center ${\bm{x}}$ and radius
proportional to $\sqrt{t}$. Usually, the radius $\sqrt{t}$ is used for one-
dimensional signals and the radius $\sqrt{2t}$ is used for two-dimensional
signals (images).

##  3 Tube-based uncertainty quantification for blob detection

The Laplacian-of-Gaussians method described in Section 2.2 detects blobs from
local minimizers of $\tilde{\Delta}u$, where $u$ is the Gaussian scale-space
representation of $f$ given by (3). The purpose of this paper is to extend
this methodology to the case where the image $f$ is uncertain, for example
when it has to be estimated from noisy measurements.

###  3.1 Incorporating uncertainty into scale-space methods

Consider the problem of recovering an image of interest $f^{*}:D\to\mathbb{R}$
(called the _ground truth_) from a noisy measurement $y$ given by

$\displaystyle y=\mathcal{G}(f^{*})+w,$ (6)

where $\mathcal{G}$ is an operator that represents the measurement process,
and $w$ is noise. The presence of the noise implies that any estimate of the
image $f^{*}$ comes with uncertainties. In the Bayesian approach Geman and
Geman (1986); Hanson (1993) to imaging, these uncertainties are taken into
account by modelling $f^{*}$, $y$ and $w$ as realizations of random quantities
$F$, $Y$ and $W$ that are related by

$\displaystyle Y=\mathcal{G}(F)+W.$ (7)

The assumed distribution ${\mathbb{P}}_{F}$ of the random image $F$ is called
the _prior_ , since it encodes a-priori assumptions on the unknown image.
Using (7) and statistical assumptions on $W$, one can define a likelihood in
the form of a conditional probability distribution ${\mathbb{P}}_{Y|F}$.
Recall that, for given $f$, ${\mathbb{P}}_{Y|F}(\cdot|f)$ is a probability
measure that represents the distribution of $Y$ given $F=f$. For the
construction of conditional probability distributions on abstract spaces (e.g.
infinite-dimensional function spaces), see for example (Kallenberg, 2002,
chapter 6) or (Klenke, 2014, chapter 8.3).

Together, the prior and the likelihood determine the so-called _posterior_
distribution ${\mathbb{P}}_{F|Y}$, which quantifies the uncertainty with
respect to $F$, conditional on observing $Y$, through Bayes’ rule (see e.g.
Kaipio and Somersalo (2005); Dashti and Stuart (2017) for further reference).

Treating the image of interest as random means that its scale-space
representation also needs to be modelled as random object. Let $\Phi$ denote
the solution operator of the diffusion equation (3) that maps an image $f$ to
its corresponding scale-space representation $u=\Phi f$. Then the scale-space
representation of the random image $F$ is the random function $U=\Phi F$. The
posterior distribution ${\mathbb{P}}_{U|Y}$ of $U$ is then determined by the
posterior distribution ${\mathbb{P}}_{F|Y}$ of $F$ through the usual
transformation rules: Given an observation $Y=y$, the posterior probability
that $U$ lies in a set of functions $A$ is given by

$\displaystyle{\mathbb{P}}_{U|Y}(A{~{}|~{}}y)$
$\displaystyle={\mathbb{P}}_{\Phi F|Y}(A{~{}|~{}}y)$
$\displaystyle={\mathbb{P}}_{F|Y}(\Phi^{-1}(A){~{}|~{}}y).$ (8)

The problem of uncertainty-aware blob detection can then be rephrased as
finding local minimizers of $\tilde{\Delta}U$ for uncertain $U$. Developing a
practical method to solve this problem requires a suitable representation of
the uncertainty encoded in the abstract posterior distribution
${\mathbb{P}}_{U|Y}$. In this paper, we assume that the uncertainty with
respect to $U$ is represented by a _credible scale-space tube_. That is, we
assume knowledge of functions
$u^{\mathrm{low}},u^{\mathrm{upp}}:D\times\mathbb{R}_{+}\to\mathbb{R}$ such
that

${\mathbb{P}}_{U|Y}([u^{\mathrm{low}},u^{\mathrm{upp}}]{~{}|~{}}y)\geq
1-\alpha,$ (9)

for a small credibility parameter $\alpha\in(0,1)$. We restrict ourselves to
tube-shaped regions for three main reasons: First, such tubes can be obtained
in many applications (see Remark 1). Second, this choice leads to a relatively
simple formulation of our method as bound-constrained convex optimization
problem (see Section 3.4.2 below). Finally, it is motivated by existing tube
methods from density estimation, which we quickly review in Section 3.2 below.

######  Remark 1.

Since $U$ depends linearly on $F$, it is often straightforward to compute
tubes that satisfy (9) using sampling-based inference (e.g. Markov chain Monte
Carlo (MCMC) Brooks et al. (2011)) or approaches based on analytic
approximations (e.g. variational inference Blei et al. (2017)). We illustrate
this in Appendix A, where we present a method for estimating the scale-space
tube using MCMC samples. This method is used in the numerical examples of
Section 5.

#### The scale-uncertainty tradeoff

The structure of the tube $[u^{\mathrm{low}},u^{\mathrm{upp}}]$ will in
general represent a tradeoff between scale and uncertainty in the image $F$:
At lower scales, the uncertainty is higher since small-scale features are more
difficult to detect from noisy observations. At higher scales, the uncertainty
decreases as this local variability is filtered out and the ability to resolve
finer details is lost. This phenomenon is a special case of the fundamental
bias-variance tradeoff that applies to uncertain signals in general Wasserman
(2006). In the language of image processing, it corresponds to the simple
intuition that coarse structures are easier to identify than fine details,
given limited data. In particular, the bounds $u^{\mathrm{low}}$ and
$u^{\mathrm{upp}}$ in general are not scale-space representations of
corresponding images $f^{\mathrm{low}},f^{\mathrm{upp}}$. This was also
demonstrated in the previous work Parzer et al. (2023).

###  3.2 Tube methods in density estimation

The task of estimating the local extrema of an uncertain signal has previously
been studied in density estimation, where it is adressed using tube methods
Mammen (1995); Mammen and van de Geer (1997); Davies and Kovac (2004);
Obereder et al. (2007). For an unknown density $g:D\to\mathbb{R}$, one
considers a tube $[g^{\mathrm{low}},g^{\mathrm{upp}}]$ that is typically
constructed as a corridor of fixed width around a noisy measurement of the
density.

The minimizer $\bar{g}\in[g^{\mathrm{low}},g^{\mathrm{upp}}]$ of a tube-
constrained optimization problem

$\displaystyle\min_{g}\quad$ $\displaystyle\mathcal{J}(g)$ s. t.
$\displaystyle g^{\mathrm{low}}\leq g\leq g^{\mathrm{upp}}$

then serves as representative with the "smallest amount of features" among all
signals in the tube, where the "amount of features" is encoded in the choice
of the cost function $\mathcal{J}$. A prominent example is the one-dimensional
taut-string algorithm Hartigan and Hartigan (1985); Davies (1995). It uses the
choice $\mathcal{J}(g)=\int\sqrt{1+\left|g^{\prime}\right|^{2}}$ which is
known to yield an estimate $\bar{g}$ such that its derivative
$\bar{g}^{\prime}$ is piecewise-constant and minimizes the number of local
extrema amongst all functions with anti-derivative in the given tube Davies
and Kovac (2001).

###  3.3 Review of the ULoG method

Motivated by these classic tube methods for density estimation, we define a
suitable cost function $\mathcal{J}(u)$ that serves as a proxy for the number
of local extrema of $\tilde{\Delta}u$, and then obtain a desired
representative $\bar{u}$ as the minimizer of the constrained optimization
problem

$\displaystyle\min_{u}$ $\displaystyle\mathcal{J}(u)$ (10) s. t.
$\displaystyle u^{\mathrm{low}}\leq u\leq u^{\mathrm{upp}}\text{ on
}D\times\mathbb{R}_{+}.$

In previous work Parzer et al. (2023) we used the cost function

$\displaystyle\mathcal{J}(u)=\int_{D\times\mathbb{R}_{+}}(\tilde{\Delta}u)^{2}.$
(11)

The particular difference to our new method is that we minimize a different
energy function, which takes into account derivatives with respect to $t$ (see
(14) below), while the optimization problem (10) can be decoupled in time,
leading to a group of bounded linear least-squares problems that can be solved
independently.

In analogy to the Laplacians-of-Gaussians method, we then identified
"significant blobs" as the local minimizers of the solution $\bar{u}$ of (10)
(see also Example 1). The resulting method was correspondingly called the
"Uncertainty-aware Laplacian-of-Gaussians" (ULoG) method.

The ULoG method showed satisfying results in simulations, where minimizers of
(10) with $\mathcal{J}$ given by (11) typically exhibited a scale-normalized
Laplacian $\tilde{\Delta}\bar{u}$ that attained its local minima at distinct
points which corresponded very well to significant blobs in the ground truth
image (cf. Parzer et al. (2023)).

However, the choice (11) was motivated more by computational simplicity than
theoretical considerations. Furthermore, a main limitation of this approach is
that it only provides a limited view of the uncertainty of scale-space blobs,
since the uncertainty with respect to a blob is represented by a single point
$({\bm{x}},t)\in D\times\mathbb{R}_{+}$. The scale parameter $t$ is difficult
to interpret, since it is both related to the expected physical size of the
blob and to the uncertainty with respect to its center. That is, a large value
of $t$ can correspond both to the presence of a large, certain blob in $f$, or
a small one with uncertain position.

###  3.4 A total variation-based approach

A remedy is to instead represent the uncertainty with respect to a blob by a
three-dimensional connected _region_ in $D\times\mathbb{R}_{+}$ that contains
the uncertain blob with high probability. The geometry of this region then
provides a more nuanced view of the possible variation in position and size.

We determine the desired regions using again a tube-based approach. To this
end, we consider a cost function that leads to minimizers of (10) with
_piecewise-constant_ normalized Laplacian. In that case, the local minima of
$\tilde{\Delta}\bar{u}$ are attained at flat connected regions which can
easily be extracted using thresholding.

As cost function, we propose to use the total variation of $\tilde{\Delta}u$
Caselles et al. (2015). This can be motivated by the fact that in one
dimension, minimizing the total variation of $u^{\prime}$ within the tube
$[u^{\mathrm{low}},u^{\mathrm{upp}}]$ leads to an estimate with very similar
properties as the mode-minimizing taut-string estimate Davies et al. (2009).
On the other hand, this approach is consistent with insights from image
processing, where total variation regularization is a common tool to achieve
piecewise-constant reconstructions Rudin et al. (1992).

####  3.4.1 Outline of the method

We start with a high-level description of the proposed method (see also the
schematic overview in Figure 1).

  1. 1.

We assume that we are given a prior distribution ${\mathbb{P}}_{F}$ for the
image of interest and a likelihood ${\mathbb{P}}_{Y|F}$ (determined by the
forward model (7) and statistical assumptions on the noise). For a concrete
observation $Y=y$, these determine a posterior distribution
${\mathbb{P}}_{F|Y}(\cdot|y)$ through Bayes’ rule (see Section 3.1), which in
turn determines a posterior distribution ${\mathbb{P}}_{U|Y}(\cdot|y)$ for the
uncertain scale-space representation $U$ (see (8)).

  2. 2.

For a given credibility level $\alpha\in(0,1)$, we compute a credible scale-
space tube $[u^{\mathrm{low}},u^{\mathrm{upp}}]$ that satisfies (9), e.g using
the MCMC-based method described in Appendix A.

  3. 3.

Next, we solve the tube-constrained total variation-based optimization problem
formulated in Section 3.4.2 below. This yields a minimizer
$\bar{u}\in[u^{\mathrm{low}},u^{\mathrm{upp}}]$.

  4. 4.

From the computed minimizer $\bar{u}$, we extract the regions in
$D\times\mathbb{R}_{+}$ where $\tilde{\Delta}\bar{u}$ attains its local
minima. In practice this is done using the thresholding procedure described in
Section 4.3.

  5. 5.

The extracted blob regions can now be visualized using the method outlined in
Section 4.4.

Since this approach can be seen as a modification of the previous ULoG method
(Section 3.3), we will refer to it as TV-ULoG.

ptforward model $Y=\mathcal{G}(F)+W$noise modellikelihood
${\mathbb{P}}_{Y|F}$prior ${\mathbb{P}}_{F}$data $y$posterior
${\mathbb{P}}_{F|Y}(\cdot|y)$${\mathbb{P}}_{U|Y}(\cdot|y)$$[u^{\mathrm{low}},u^{\mathrm{upp}}]$$\bar{u}$blob
regionsBayes’ rulepushforwardptcomputation of credible tubes, see e.g.
Appendix Aptsolution of the optimization problem (14)extraction procedure
Figure 1: Overview of the proposed approach

####  3.4.2 Formulation of the optimization problem

To prepare the precise mathematical formulation of the resulting optimization
problem, we define the scale-normalized total variation of a function
$u:D\times\mathbb{R}_{+}\to\mathbb{R}$ by

$\displaystyle\tilde{\mathrm{TV}}(u):=\int_{D\times\mathbb{R}_{+}}\left\|\tilde{\nabla}_{{\bm{x}},t}u\right\|,$
(12)

where

$\displaystyle\tilde{\nabla}_{{\bm{x}},t}:=\begin{bmatrix}\tilde{\partial}_{x_{1}}\\\
\tilde{\partial}_{x_{2}}\\\ \tilde{\partial}_{t}\end{bmatrix}$ (13)

is the scale-normalized gradient operator. Here, $\tilde{\partial}_{x_{1}}$
and $\tilde{\partial}_{x_{2}}$ are the scale-normalized spatial derivatives
defined in (4), while $\tilde{\partial}_{t}$ is defined as

$\displaystyle\tilde{\partial}_{t}u({\bm{x}},t):=t\partial_{t}u({\bm{x}},t).$

Following the discussion at the start of this section, we suggest to use the
choice $\mathcal{J}(u)=\tilde{\mathrm{TV}}(\tilde{\Delta}u)$ in (10). To allow
for a finite-difference discretization of the involved differential operators,
we have to assume suitable boundary conditions. Our particular choice of
Neumann boundary conditions on $u$ and $\tilde{\Delta}u$ was mostly motivated
by ease of implementation (see Remark 3). In summary, we arrive at the
following formulation:

$\displaystyle\min_{u}$ $\displaystyle\tilde{\mathrm{TV}}(\tilde{\Delta}u)$
(14) s. t. $\displaystyle u^{\mathrm{low}}\leq u\leq u^{\mathrm{upp}}\text{ on
}D\times\mathbb{R}_{+},$ $\displaystyle\partial_{\bm{\nu}}u(\bm{x},t)=0\text{
on }\partial D\times(0,\infty),$
$\displaystyle\partial_{\bm{\nu}}\tilde{\Delta}u(\bm{x},t)=0\text{ on
}\partial D\times(0,\infty).$

In Section 4.2, we discuss the discretization of (14) and consider various
optimization algorithms for finding a minimizer of the resulting non-smooth
convex optimization problem. Extraction and visualization of the desired
regions from the computed minimizer are described in Section 4.3 and Section
4.4.

######  Remark 2.

The scale-normalization in (13) is necessary to achieve a scale-invariance
property analogous to the one described in Section 2.2.1: If $\bar{u}_{s}$ is
a minimizer of the scaled problem

$\displaystyle\min_{u_{s}}\quad$
$\displaystyle\tilde{\mathrm{TV}}(\tilde{\Delta}u_{s})$ s. t. $\displaystyle
u_{s}^{\mathrm{low}}\leq u_{s}\leq u_{s}^{\mathrm{upp}}\text{ on
}D\times\mathbb{R}_{+},$
$\displaystyle\partial_{\bm{\nu}}u_{s}(\bm{x},t)=0\text{ on }\partial
D\times\mathbb{R}_{+},$
$\displaystyle\partial_{\bm{\nu}}\tilde{\Delta}u_{s}(\bm{x},t)=0\text{ on
}\partial D\times\mathbb{R}_{+},$

where

$\displaystyle
u_{s}^{\mathrm{low}}({\bm{x}},t)=u^{\mathrm{low}}(\sqrt{s}{\bm{x}},st),$
$\displaystyle
u_{s}^{\mathrm{upp}}({\bm{x}},t)=u^{\mathrm{upp}}(\sqrt{s}{\bm{x}},st),$

then $\bar{u}$ given by

$\displaystyle\bar{u}({\bm{x}},t)=\bar{u}_{s}({\bm{x}}/\sqrt{s},t/s)\qquad\text{for
all }({\bm{x}},t)\in D\times\mathbb{R}_{+}$

is a minimizer of the original problem (14) and vice versa. See also Lindeberg
(1998) for more motivation behind such scaling properties.

######  Remark 3.

The Neumann boundary condition on $u$ in (14) is motivated by the definition
of the scale-space representation (3). The motivation for the second Neumann
boundary condition on $\tilde{\Delta}u$ is that it leads to the usual formula
for the discrete total variation when combined with a forward-difference
approximation (see Section 4.2).

##  4 Numerical implementation

In this section we discuss how the TV-ULoG method is implemented in practice.
In Section 4.1, we discuss the discretization of the optimization problem
(14). Then, we present methods to solve the resulting discrete optimization
problem in Section 4.2. Section 4.3 describes how the desired blob regions can
be extracted from the computed minimizer, while Section 4.4 outlines possible
visualization of these regions.

###  4.1 Discrete scale-space representation

To discuss the numerical implementation of the TV-ULoG method, we assume that
the image domain $D\subset\mathbb{R}^{2}$ is rectangular and has been
discretized into a grid $(\bm{x}_{ij})_{i\in[N_{1}],j\in[N_{2}]}$, with
uniform grid size $h_{1}>0$ in $x_{1}$-direction and uniform grid size
$h_{2}>0$ in $x_{2}$-direction, such that the discrete image of interest is
given by a matrix ${\bm{f}}\in\mathbb{R}^{N_{1}\times N_{2}}$.

For the scale discretization, it was suggested in Lindeberg (1998) to use
discrete scales $0<t_{1}<\ldots<t_{K}$ that are exponentially increasing, i.e.

$\displaystyle t_{k+1}=bt_{k},\quad k\in[K-1]$ (15)

for some $b>1$.

The discrete scale-space representation is then defined through a suitable
discretization $\bm{\Phi}:\mathbb{R}^{N_{1}\times
N_{2}}\to\mathbb{R}^{N_{1}\times N_{2}\times K}$ of the solution operator
$\Phi$ of the diffusion equation (3). Then, the discrete scale-space
representation of ${\bm{f}}$ is given by the three-dimensional array

$\displaystyle\bm{u}:=\bm{\Phi}{\bm{f}}\in\mathbb{R}^{N_{1}\times N_{2}\times
K}.$ (16)

In practice, the solution operator $\bm{\Phi}$ is often implemented as a
convolution with a suitable discrete convolution kernel (see e.g. Lindeberg
(1994)).

Following the discussion in Section 3.1, we then assume access to a credible
scale-space tube $[{\bm{u}}^{\mathrm{low}},{\bm{u}}^{\mathrm{upp}}]$ for the
discrete scale-space representation $\bm{U}=\bm{\Phi}{\bm{F}}$, such that,
analogously to (9), there holds

$\displaystyle{\mathbb{P}}_{{\bm{U}}|{\bm{Y}}}([{\bm{u}}^{\mathrm{low}},{\bm{u}}^{\mathrm{upp}}]{~{}|~{}}{\bm{y}})\geq
1-\alpha,$

for given $\alpha\in(0,1)$. How such a tube can be estimated in practice from
MCMC samples is described in Appendix A.

In Section 4.2, we present how, given
$[{\bm{u}}^{\mathrm{low}},{\bm{u}}^{\mathrm{upp}}]$, the non-smooth
optimization problem (14) can be solved numerically after discretization.
Then, in Section 4.3, we describe a procedure that extracts the desired blob
regions in scale space from the computed minimizer in a way that is robust
against numerical errors. In Section 4.4, we discuss visualizations of these
regions that meaningfully represent the uncertainty of the scale-space blobs.

###  4.2 Numerical treatment of the optimization problem

In order to solve the optimization problem (14) numerically, we have to define
suitable discretizations for the differential operators in the objective.

To discretize the scale-normalized Laplacian (5), we use the common central-
difference scheme, i.e. we define

$\displaystyle(\tilde{\bm{\Delta}}\bm{u})_{i,j,k}:=t_{k}$
$\displaystyle\left(\frac{u_{i+1,j,k}-2u_{i,j,k}+u_{i-1,j,k}}{h_{1}}\right.$
(17)
$\displaystyle\left.+\frac{u_{i,j+1,k}-2u_{i,j,k}+u_{i,j-1,k}}{h_{2}}\right),$

where we mirror $u$ at the boundaries of the index range, i.e. we set

$\displaystyle u_{0,j,k}:=u_{2,j,k},\quad u_{N_{1}+1,j,k}:=u_{N_{1}-1,j,k},$
(18) $\displaystyle u_{i,0,k}:=u_{j,2,k},\quad
u_{i,N_{2}+1,k}:=u_{j,N_{2}-1,k}$

for all $(i,j,k)\in[N_{1}]\times[N_{2}]\times[K]$. This choice implements the
Neumann boundary condition in (3) (see example 5.49 in Bredies and Lorenz
(2018)).

For the total variation (12), we use an isotropic discretization. To this end,
we define the forward differences in $x_{1}$\- and $x_{2}$-direction by

$\displaystyle(\tilde{\bm{\nabla}}_{{\bm{x}},t}\bm{u})_{i,j,k}:=\begin{bmatrix}(\tilde{\bm{\nabla}}_{x_{1}}\bm{u})_{i,j,k}\\\
(\tilde{\bm{\nabla}}_{x_{2}}\bm{u})_{i,j,k}\\\
(\tilde{\bm{\nabla}}_{t}\bm{u})_{i,j,k}\end{bmatrix}\in\mathbb{R}^{3},$ where
$\displaystyle(\tilde{\bm{\nabla}}_{x_{1}}\bm{u})_{i,j,k}:=t_{k}^{1/2}\frac{u_{i+1,j,k}-u_{i,j,k}}{h_{1}},$
$\displaystyle(\bm{\tilde{\bm{\nabla}}}_{x_{2}}\bm{u})_{i,j,k}:=t_{k}^{1/2}\frac{u_{i,j+1,k}-u_{i,j,k}}{h_{2}},$

where we formally define

$\displaystyle u_{N_{1}+1,j,k}:=u_{N_{1},j,k},\quad
u_{i,N_{2}+1,k}:=u_{i,N_{2},k},$ $\displaystyle u_{i,j,K+1}:=u_{i,j,K},$

for all $(i,j,k)\in[N_{1}]\times[N_{2}]\times[K]$. Similar to (18), this
choice implements the Neumann boundary condition on $\tilde{\Delta}u$.

For the non-uniform scale grid (15), the forward difference approximation of
the scale derivative $\partial_{t}u({\bm{x}},t)$ is

$\displaystyle(\tilde{\bm{\nabla}}_{t}\bm{u})_{i,j,k}$
$\displaystyle=t_{k}\frac{u_{i,j,k+1}-u_{i,j,k}}{t_{k+1}-t_{k}}$
$\displaystyle=\frac{1}{b-1}(u_{i,j,k+1}-u_{i,j,k}).$

We combine the forward difference approximations in a discrete scale-space
gradient $\tilde{\bm{\nabla}}_{{\bm{x}},t}:\mathbb{R}^{N_{1}\times N_{2}\times
K}\to\mathbb{R}^{N_{1}\times N_{2}\times K\times 3}$, given by

$(\tilde{\bm{\nabla}}_{{\bm{x}},t}\bm{u})_{i,j,k}:=\begin{bmatrix}(\tilde{\bm{\nabla}}_{x_{1}}\bm{u})_{i,j,k}\\\
(\tilde{\bm{\nabla}}_{x_{2}}\bm{u})_{i,j,k}\\\
(\tilde{\bm{\nabla}}_{t}\bm{u})_{i,j,k}\end{bmatrix}\in\mathbb{R}^{3}.$ (19)

With definition (19), we define the isotropic scale-normalized total variation
of $\bm{u}\in\mathbb{R}^{N_{1}\times N_{2}\times K}$ as

$\displaystyle\tilde{\mathbf{TV}}(\bm{u})=\sum_{k=1}^{K}\sum_{i=1}^{N_{1}}\sum_{j=1}^{N_{2}}\left\|(\tilde{\bm{\nabla}}_{{\bm{x}},t}\bm{u})_{i,j,k}\right\|.$
(20)

In summary, our discretization of (14) reads as

$\displaystyle\min_{\bm{u}\in\mathbb{R}^{N_{1}\times N_{2}\times K}}$
$\displaystyle\tilde{\mathbf{TV}}(\tilde{\bm{\Delta}}\bm{u})$ (21) s.t.
$\displaystyle{\bm{u}}^{\mathrm{low}}\leq\bm{u}\leq{\bm{u}}^{\mathrm{upp}}.$

This is a bound-constrained, non-smooth convex optimization problem. For ease
of reference, we will also refer to (21) as the _TV-ULoG optimization
problem_.

To our knowledge, this particular problem has not been previously considered
in the image processing literature. The closest candidate is the constrained
total variation (CTV) denoising problem: Given a noisy image
$\bm{f}^{\delta}\in\mathbb{R}^{N_{1}\times N_{2}}$ and bounds
${\bm{f}}^{\mathrm{low}},{\bm{f}}^{\mathrm{upp}}\in\mathbb{R}^{N_{1}\times
N_{2}}$ with ${\bm{f}}^{\mathrm{low}}\leq{\bm{f}}^{\mathrm{upp}}$, the
discrete CTV-denoising problem is

$\displaystyle\min_{\bm{f}\in\mathbb{R}^{N_{1}\times N_{2}}}$
$\displaystyle\left\|\bm{f}-\bm{f}^{\delta}\right\|^{2}+\lambda\mathbf{TV}(\bm{f})$
(22) s. t.
$\displaystyle{\bm{f}}^{\mathrm{low}}\leq\bm{f}\leq{\bm{f}}^{\mathrm{upp}}.$

Here, $\lambda>0$ is a tunable regularization parameter and
$\mathbf{TV}(\cdot)$ is the discrete total variation. The problem (22) has
been considered e.g. by Beck and Teboulle in their seminal work Beck and
Teboulle (2009), where they propose the fast gradient projection (FGP) method
for its numerical solution. The FGP method can be summarized as applying
Nesterov-accelerated projected gradient descent to the dual of (22). In
Section 4.2.1, we show that the FGP method can be applied to a Nesterov-
smoothed dual of (21). In Section 4.2.2, we describe an analogous method that
instead applies FGP to the Nesterov-smoothed primal problem.

However, as is illustrated further in Section 5, both methods converge very
slowly for our particular problem. For this reason, we present a third
alternative approach for the solution of (21) in Section 4.2.3. It formulates
(21) as an equivalent second-order cone program (SOCP) and solves it with an
interior-point method. This method is also motivated by analogous approaches
for TV-minimization Goldfarb and Yin (2005); Yin et al. (2007).

####  4.2.1 Dual smoothing

To prepare the discussion of the dual-smoothing approach, we recall the notion
of Fenchel duality Bot et al. (2009): Let $\mathbb{U},\mathbb{V}$ be finite-
dimensional Euclidean spaces, $\bm{A}:\mathbb{U}\to\mathbb{V}$ a linear
function, and $\phi:\mathbb{U}\to[-\infty,\infty]$,
$\psi:\mathbb{V}\to[-\infty,\infty]$ be proper convex functions. Then the
Fenchel dual of the optimization problem

$\displaystyle\min_{{\bm{u}}\in\mathbb{U}}\left\\{\phi({\bm{u}})+\psi(\bm{A}{\bm{u}})\right\\}$
(23)

is

$\displaystyle\max_{{\bm{v}}\in\mathbb{V}}\left\\{-\phi^{*}(-\bm{A}^{*}{\bm{v}})-\psi^{*}({\bm{v}})\right\\},$
(24)

where $\phi^{*}$ and $\psi^{*}$ denote the convex conjugates of $\phi$ and
$\psi$, respectively.

Let now $\mathbb{U}=\mathbb{R}^{N_{1}\times N_{2}\times K}$ and
$\mathbb{V}=\mathbb{R}^{N_{1}\times N_{2}\times K\times 3}$, such that (21) is
of the form (23) with

$\displaystyle\phi({\bm{u}}):=\chi_{[{\bm{u}}^{\mathrm{low}},{\bm{u}}^{\mathrm{upp}}]}({\bm{u}}),\qquad{\bm{u}}\in\mathbb{R}^{N_{1}\times
N_{2}\times K},$ (25)
$\displaystyle\psi({\bm{v}}):=\sum_{i,j,k}\left\|{\bm{v}}_{i,j,k}\right\|,\qquad{\bm{v}}\in\mathbb{R}^{N_{1}\times
N_{2}\times K\times 3},$ (26)
$\displaystyle{\bm{A}}:=\tilde{\bm{\nabla}}_{{\bm{x}},t}\tilde{\bm{\Delta}}:\mathbb{R}^{N_{1}\times
N_{2}\times K}\to\mathbb{R}^{N_{1}\times N_{2}\times K\times 3}.$ (27)

The next proposition gives an explicit expression for the dual of the TV-ULoG
optimization problem (21). It is defined with reference to the convex set

$\displaystyle S:=\left\\{{\bm{v}}\in\mathbb{R}^{N_{1}\times N_{2}\times
K\times 3}~{}:~{}\left\|{\bm{v}}_{i,j,k}\right\|\leq 1\text{ for all
}i,j,k\right\\}.$ (28)

######  Proposition 1.

The Fenchel dual of (21) is

$\max_{{\bm{v}}\in
S}\min_{{\bm{w}}\in[{\bm{u}}^{\mathrm{low}},{\bm{u}}^{\mathrm{upp}}]}\left<{\bm{A}}^{\top}{\bm{v}},{\bm{w}}\right>,$
(29)

where ${\bm{A}}$ is as in (27).

###### Proof.

Let $\phi$ and $\psi$ be given by (25) and (26), respectively. By example 4.2
in Beck (2017), we have

$\displaystyle\phi^{*}(\bm{v})=\max_{{\bm{w}}\in[{\bm{u}}^{\mathrm{low}},{\bm{u}}^{\mathrm{upp}}]}\left<{\bm{w}},{\bm{v}}\right>.$
(30)

Furthermore, by theorem 4.12 and example 4.4.12 in Beck (2017),

$\displaystyle\psi^{*}({\bm{v}})$
$\displaystyle=\sum_{i,j,k}\chi_{B_{1}(0)}({\bm{v}}_{i,j,k})$
$\displaystyle=\chi_{S}({\bm{v}}),$ (31)

where the set $S\subset\mathbb{R}^{N_{1}\times N_{2}\times K\times 3}$ is
given by (28). The proof follows if we plug (30) and (31) into (24). ∎

The dual problem (29) could be solved with projected subgradient methods Beck
(2017). However, a more efficient approach is to consider a smoothed
approximation instead, since this allows to apply accelerated methods such as
FGP: Given a convex optimization problem of the form

$\min_{{\bm{u}}\in C_{1}}~{}\phi({\bm{u}}),\text{ where
}\phi({\bm{u}})=\max_{{\bm{w}}\in C_{2}}\left<\bm{B}{\bm{u}},{\bm{w}}\right>,$
(32)

where $C_{1}$ and $C_{2}$ are bounded closed convex sets, Nesterov Nesterov
(2005) proposed to approximate the objective functional $\phi$ by

$\displaystyle\phi_{\mu}(\bm{u}):=\max_{{\bm{w}}\in
C_{2}}\left\\{\left<\bm{B}{\bm{u}},{\bm{w}}\right>-\frac{\mu}{2}\left\|{\bm{w}}\right\|^{2}\right\\}.$
(33)

The associated optimization problem

$\displaystyle\min_{{\bm{u}}\in C_{1}}\phi_{\mu}({\bm{u}})$

is called the Nesterov smoothing of (32). It is a convex problem with smooth
objective, and can hence be solved fast using accelerated first-order methods.

The next proposition derives the Nesterov smoothing of the dual problem (29).
The derivation is analogous to the proof of proposition 4.1 in Beck and
Teboulle (2009). For completeness, we have provided the proof below.

######  Proposition 2.

The Nesterov smoothing corresponding to (29) is given by

$\min_{{\bm{v}}\in
S}\quad\left\\{\left\|\frac{1}{\mu}{\bm{A}}^{\top}{\bm{v}}\right\|^{2}-\left\|(I-P_{[{\bm{u}}^{\mathrm{low}},{\bm{u}}^{\mathrm{upp}}]})(\frac{1}{\mu}{\bm{A}}^{\top}{\bm{v}})\right\|^{2}\right\\}.$
(34)

###### Proof.

First, (29) is equivalent to

$\displaystyle\min_{{\bm{v}}\in
S}\max_{{\bm{w}}\in[{\bm{u}}^{\mathrm{low}},{\bm{u}}^{\mathrm{upp}}]}\left<-{\bm{A}}^{\top}{\bm{v}},{\bm{w}}\right>.$

This is of the form (32) with $C_{1}=S$,
$C_{2}=[{\bm{u}}^{\mathrm{low}},{\bm{u}}^{\mathrm{upp}}]$ and
$\bm{B}=-{\bm{A}}^{\top}$. Hence, its Nesterov smoothing is

$\displaystyle\min_{{\bm{v}}\in K}\quad\phi_{\mu}({\bm{v}}),$
$\displaystyle\phi_{\mu}({\bm{v}}):=\max_{{\bm{w}}\in[{\bm{u}}^{\mathrm{low}},{\bm{u}}^{\mathrm{upp}}]}\left\\{\left<-{\bm{A}}^{\top}{\bm{v}},{\bm{w}}\right>-\frac{\mu}{2}\left\|{\bm{w}}\right\|^{2}\right\\}.$

By completing the square, one can show

$\displaystyle\phi_{\mu}({\bm{v}})=\max_{{\bm{w}}\in[{\bm{u}}^{\mathrm{low}},{\bm{u}}^{\mathrm{upp}}]}$
$\displaystyle\left\\{\frac{\mu}{2}\left\|\frac{1}{\mu}{\bm{A}}^{\top}{\bm{v}}\right\|^{2}\right.$
(35)
$\displaystyle\quad\left.-\frac{\mu}{2}\left\|{\bm{w}}+\frac{1}{\mu}{\bm{A}}^{\top}{\bm{v}}\right\|^{2}\right\\}.$

Recall the distance minimizing property of the orthogonal projection: For a
convex set $C$, we have

$\displaystyle\min_{{\bm{w}}\in
C}\left\|\bm{q}-{\bm{w}}\right\|^{2}=\left\|(I-P_{C})(\bm{q})\right\|^{2}.$

In particular, this means

$\displaystyle\max_{{\bm{w}}\in[{\bm{u}}^{\mathrm{low}},{\bm{u}}^{\mathrm{upp}}]}\left\|{\bm{w}}+\frac{1}{\mu}{\bm{A}}^{\top}{\bm{v}}\right\|^{2}$
$\displaystyle\qquad=\left\|(I-P_{[{\bm{u}}^{\mathrm{low}},{\bm{u}}^{\mathrm{upp}}]})(-\frac{1}{\mu}{\bm{A}}^{\top}{\bm{v}}))\right\|^{2}.$

Inserting (4.2.1) in (35) then yields

$\displaystyle\phi_{\mu}({\bm{v}})=$
$\displaystyle\frac{\mu}{2}\left(-\left\|(I-P_{[{\bm{u}}^{\mathrm{low}},{\bm{u}}^{\mathrm{upp}}]})(-\frac{1}{\mu}{\bm{A}}^{\top}{\bm{v}})\right\|^{2}\right.$
$\displaystyle\quad\left.+\frac{1}{\mu}\left\|{\bm{A}}^{\top}{\bm{v}}\right\|^{2}\right).$

This finishes the proof. ∎

The smoothed problem (34) is now amenable to minimization with the fast
gradient projection method given by Algorithm 1 (cf. Beck (2017)).

Algorithm 1 FGP

1:A convex function $\phi:\mathbb{U}\to\mathbb{R}$ on a Euclidean space
$\mathbb{U}$, a convex set $C\subset\mathbb{U}$, a steplength $\beta>0$ and an
initial guess ${\bm{u}}_{0}\in\mathbb{U}$.

2:${\bm{v}}_{0}={\bm{u}}_{0}$;

3:$t_{0}=1$;

4:for k = 0,…,K do

5: ${\bm{u}}_{k+1}=P_{C}({\bm{v}}_{k}-\beta\nabla\phi({\bm{v}}_{k}))$;

6: $t_{k+1}=\frac{1}{2}(1+\sqrt{1+4t_{k}^{2}})$;

7:
${\bm{v}}_{k+1}={\bm{u}}_{k+1}+\frac{t_{k}-1}{t_{k}+1}({\bm{u}}_{k+1}-{\bm{u}}_{k})$;

8:end for

######  Remark 4.

Problem (34) is identical to the dual problem of CTV-denoising (see
proposition 4.1 in Beck and Teboulle (2009)), except for the definition of the
operator ${\bm{A}}$. Indeed, problem (34) is the Fenchel dual of the following
regularized version of problem (21),

$\displaystyle\min_{\bm{u}\in\mathbb{R}^{N_{1}\times N_{2}\times K}}\quad$
$\displaystyle\tilde{\mathbf{TV}}(\tilde{\bm{\Delta}}\bm{u})+\frac{1}{2\mu}\left\|\bm{u}\right\|^{2}$
s.t.
$\displaystyle{\bm{u}}^{\mathrm{low}}\leq\bm{u}\leq{\bm{u}}^{\mathrm{upp}}.$

######  Remark 5.

Finding a good value for the smoothing parameter $\mu$ in (33) is difficult in
practice. In Nesterov (2005), it is shown that, for all
${\bm{u}}\in\mathbb{U}$,

$\displaystyle\phi_{\mu}({\bm{u}})\leq\phi({\bm{u}})\leq\phi_{\mu}({\bm{u}})+\mu
D_{2},$ (36)

where $D_{2}=\max_{{\bm{v}}\in C_{2}}\frac{1}{2}\left\|{\bm{v}}\right\|^{2}$.
This means that an accuracy of $\epsilon$ in the objective is ensured if
$\mu\leq\epsilon/D_{2}$. However, the bound (36) is often too conservative and
it can be advantageous to choose $\mu>\epsilon/D_{2}$, since larger values of
$\mu$ allow for larger stepsize and faster convergence.

####  4.2.2 Primal smoothing

As an alternative to the dual approach described in Section 4.2.1, we can also
apply Nesterov-smoothing directly to the primal problem (21).

######  Proposition 3.

Let ${\bm{A}}$ be given by (27) and define the Huber loss function Huber
(1981)

$h_{\mu}:\mathbb{R}^{3}\to\mathbb{R},\quad
h_{\mu}(\bm{v}):=\begin{cases}\left\|\bm{v}\right\|+\frac{\mu}{2},&\text{if
}\left\|\bm{v}\right\|\geq\mu,\\\
\frac{\left\|\bm{v}\right\|^{2}}{2\mu},&\text{otherwise}.\end{cases}$ (37)

Then the Nesterov smoothing of the TV-ULoG optimization problem (21) is

$\displaystyle\min_{\bm{u}\in\mathbb{R}^{N_{1}\times N_{2}\times K}}$
$\displaystyle\sum_{i,j,k}h_{\mu}(({\bm{A}}\bm{u})_{i,j,k}),$ (38) s.t.
$\displaystyle{\bm{u}}^{\mathrm{low}}\leq\bm{u}\leq{\bm{u}}^{\mathrm{upp}}.$

###### Proof.

Let

$\displaystyle\phi(\bm{u}):=\sum_{i,j,k}\left\|(\bm{A}\bm{u})_{i,j,k}\right\|,$
(39)

such that we can write the TV-ULoG optimization problem (21) as

$\displaystyle\min_{\bm{u}\in[{\bm{u}}^{\mathrm{low}},{\bm{u}}^{\mathrm{upp}}]}\phi(\bm{u}).$
(40)

Using the identity

$\displaystyle\left\|{\bm{v}}\right\|=\max_{{\bm{w}}\in
B_{1}(0)}\left<{\bm{w}},{\bm{v}}\right>,$

in (39) yields

$\displaystyle\phi({\bm{u}})$
$\displaystyle=\sum_{i,j,k}\max_{{\bm{w}}_{i,j,k}\in
B_{1}(0)}\left<{\bm{w}}_{i,j,k},({\bm{A}}{\bm{u}})_{i,j,k}\right>$
$\displaystyle=\max_{{\bm{w}}\in S}\left<{\bm{w}},\bm{A}\bm{u}\right>,$

where $S$ is again given by (28). Hence, (40) is of the form (32) with
$C_{1}=[{\bm{u}}^{\mathrm{low}},{\bm{u}}^{\mathrm{upp}}]$ and $C_{2}=S$, which
means that the Nesterov-smoothed objective is given by (cf. (33))

$\displaystyle\phi_{\mu}(\bm{u})$ $\displaystyle=\max_{{\bm{w}}\in
S}\left(\left<{\bm{w}},{\bm{u}}\right>-\frac{\mu}{2}\left\|{\bm{w}}\right\|^{2}\right)$
$\displaystyle=\sum_{i,j,k}\max_{{\bm{w}}_{i,j,k}\in
B_{1}(0)}\left(\left<{\bm{w}}_{i,j,k},(\bm{A}{\bm{u}})_{i,j,k}\right>-\frac{\mu}{2}\left\|{\bm{w}}_{i,j,k}\right\|^{2}\right).$
(41)

It is straightforward to show (e.g. using Lagrange multipliers) that, for any
vector ${\bm{v}}$,

$\displaystyle\max_{{\bm{w}}\in
B_{1}(0)}\left(\left<{\bm{v}},{\bm{w}}\right>-\frac{\mu}{2}\left\|{\bm{w}}\right\|^{2}\right)=h_{\mu}({\bm{v}}),$
(42)

where $h_{\mu}$ is defined in (37). Using (42) in (41) proves the proposition.
∎

Problem (38) is a smooth optimization problem with bound constraints for which
there exists a huge number of optimization methods, such as the previously
introduced FGP method or the popular L-BFGS-B method Byrd et al. (1995); Zhu
et al. (1997); Morales and Nocedal (2011).

####  4.2.3 Interior-point method

Many problems involving TV-regularization, in particular $L^{1}$-TV denoising
and CTV-denoising, can be formulated equivalently as second-order cone
programs (SOCP) and then be solved with interior-point methods Yin et al.
(2007). This strategy has the advantage that it does not require additional
smoothing. Interior-point methods are also known to be much more robust
against ill-conditioning of the KKT system. They require the solution of a
large linear system in every step, which is why they sometimes do not scale
well to larger problems. However, if the special structure of the linear
system can be exploited, interior-point methods can yield state-of-the-art
performance Goldfarb and Yin (2005).

We call an optimization problem a SOCP if it can be written in the following
form Boyd and Vandenberghe (2004):

$\displaystyle\min_{{\bm{v}}\in\mathbb{R}^{n}}$
$\displaystyle\left<\bm{\xi},{\bm{v}}\right>$ (43) s. t.
$\displaystyle\left\|\bm{B}_{i}{\bm{v}}+\bm{c}_{i}\right\|\leq\left<\bm{d}_{i},{\bm{v}}\right>+\eta_{i},\quad
i\in[m],$ $\displaystyle\bm{H}{\bm{v}}=\bm{h},$

where $\bm{\xi},\bm{d}_{i}\in\mathbb{R}^{n}$,
$\bm{B}_{i}\in\mathbb{R}^{n_{i}\times n}$, $\bm{c}_{i}\in\mathbb{R}^{n_{i}}$,
$\eta_{i}\in\mathbb{R}$, $\bm{H}\in\mathbb{R}^{p\times n}$, and
$\bm{h}\in\mathbb{R}^{p}$.

The next proposition states how the TV-ULoG optimization problem (21) can be
brought into the form (43). To prepare the proof, we bring the primal problem
in flattened form so that we work with vectors instead of three-dimensional
arrays. To this end, let $N:=N_{1}\cdot N_{2}\cdot K$ and define the
flattening operator

$\displaystyle\mathrm{flat}:\mathbb{R}^{N_{1}\times N_{2}\times
K}\to\mathbb{R}^{N},$
$\displaystyle\mathrm{flat}(\bm{u})_{\sigma(i,j,k)}=u_{i,j,k},$

where

$\displaystyle\sigma:[N_{1}]\times[N_{2}]\times[K]\to[N],$
$\displaystyle\sigma(i,j,k)=iN_{2}K+jK+k$

is an enumeration of $[N_{1}]\times[N_{2}]\times[K]$. Let

$\displaystyle\bm{z}^{\mathrm{low}}:=\mathrm{flat}({\bm{u}}^{\mathrm{low}}),\quad\bm{z}^{\mathrm{upp}}:=\mathrm{flat}({\bm{u}}^{\mathrm{upp}}).$
(44)

Then, one can find matrices
$\tilde{\bm{A}}_{1},\ldots,\tilde{\bm{A}}_{N}\in\mathbb{R}^{3\times N}$ such
that

$\displaystyle(\tilde{\bm{A}}_{\sigma(i,j,k)}\mathrm{flat}({\bm{u}}))_{r}=({\bm{A}}{\bm{u}})_{i,j,k,r}$
(45)

for all ${\bm{u}}\in\mathbb{R}^{N_{1}\times N_{2}\times K}$, $r\in[3]$ and
$(i,j,k)\in[N_{1}]\times[N_{2}]\times[K]$. Under these definitions, it is easy
to see that the TV-ULoG optimization problem (21) is equivalent to

$\displaystyle\min_{\bm{z}\in\mathbb{R}^{N}}$
$\displaystyle\sum_{\ell=1}^{N}\left\|\tilde{\bm{A}}_{\ell}\bm{z}\right\|$
(46) s. t.
$\displaystyle\bm{z}^{\mathrm{low}}\leq\bm{z}\leq\bm{z}^{\mathrm{upp}},$

where $\bm{z}=\mathrm{flat}({\bm{u}})$.

######  Proposition 4.

The optimization problem (46) is equivalent to the SOCP (43), under the
identifications

$\displaystyle n=2N,\quad m=3N,\quad\bm{\xi}=\begin{bmatrix}\bm{0}_{N}\\\
\bm{1}_{N}\end{bmatrix},$ (47) $\displaystyle\bm{B}_{\ell}=\bm{0}_{1\times
2N},\quad\bm{B}_{N+\ell}=\bm{0}_{1\times
n},\quad\bm{B}_{2N+\ell}=\begin{bmatrix}\tilde{\bm{A}}_{\ell}&\bm{0}_{3\times
n}\end{bmatrix},$
$\displaystyle\bm{c}_{\ell}=0,\quad\bm{c}_{N+\ell}=0,\quad\bm{c}_{2N+\ell}=\bm{0}_{3},$
$\displaystyle\bm{d}_{\ell}=\bm{e}_{\ell}^{n},\quad\bm{d}_{N+\ell}=-\bm{e}_{\ell}^{n},\quad\bm{d}_{2N+i}=\bm{e}_{N+\ell},$
$\displaystyle\eta_{\ell}=-z^{\mathrm{low}}_{\ell},\quad\eta_{N+\ell}=z^{\mathrm{upp}}_{\ell},\quad\eta_{2N+\ell}=0,\quad\ell\in[N],$
$\displaystyle\bm{H}=\bm{0}_{1\times n},\quad\bm{h}=0.$

In particular, if $\bar{{\bm{v}}}\in\mathbb{R}^{2N}$ is a minimizer of (43),
then $\bar{{\bm{u}}}$ given by

$\displaystyle\bar{u}_{i,j,k}=\bar{v}_{\sigma(i,j,k)}$

is a minimizer of the TV-ULoG optimization problem (21).

###### Proof.

Note that (46) is equivalent to

$\displaystyle\min_{\bm{z}\in\mathbb{R}^{N},\bm{q}\in\mathbb{R}^{N}}$
$\displaystyle\sum_{\ell=1}^{N}q_{i}$ (48) s. t.
$\displaystyle\bm{z}^{\mathrm{low}}\leq\bm{z}\leq\bm{z}^{\mathrm{upp}},$
$\displaystyle\left\|\tilde{\bm{A}}_{\ell}\bm{z}\right\|\leq
q_{\ell},\quad\ell\in[N].$

Let us combine the optimization variables $\bm{z}$ and $\bm{q}$ in a single
vector ${\bm{v}}=[\bm{z},\bm{q}]\in\mathbb{R}^{2N}$, such that (48) becomes

$\displaystyle\min_{{\bm{v}}\in\mathbb{R}^{2N}}$
$\displaystyle\sum_{\ell=1}^{N}v_{N+\ell},$ (49) s. t. $\displaystyle
z^{\mathrm{low}}_{\ell}\leq v_{\ell}\leq z^{\mathrm{upp}}_{\ell},$
$\displaystyle\left\|\begin{bmatrix}\tilde{\bm{A}}_{\ell}&\bm{0}_{3\times
N}\end{bmatrix}{\bm{v}}\right\|\leq v_{N+\ell},\quad\ell\in[N].$

It is now straightforward to check that if we plug the identifications (47)
into the SOCP standard form (43), we obtain precisely (49). ∎

The resulting SOCP can then be solved efficiently with existing interior-point
solvers (see also Section 5.3), exploiting the sparse structure of the
discretized differential operator ${\bm{A}}$.

###  4.3 Extraction of blob regions

Let $\bar{\bm{u}}$ be a numerical solution of the TV-ULoG optimization problem
(21). Ideally, $\tilde{\bm{\Delta}}\bar{\bm{u}}$ is piecewise constant such
that it attains its local minima on index sets
$\mathcal{M}_{1},\ldots,\mathcal{M}_{S}\subset[N_{1}]\times[N_{2}]\times[K]$
which quantify our uncertainty with respect to the blobs of the uncertain
image (see Section 3).

However, due to numerical errors, $\tilde{\bm{\Delta}}\bar{\bm{u}}$ will
typically only be approximately piecewise constant or exhibit artifacts such
as staircasing. For this reason, we use the following thresholding procedure
to extract the desired regions.

  1. 1.

Let $\bm{a}:=\tilde{\bm{\Delta}}\bar{\bm{u}}\in\mathbb{R}^{N_{1}\times
N_{2}\times K}$.

  2. 2.

Detect the local minimizers
$\bm{m}_{s},\ldots,\bm{m}_{S}\in[N_{1}]\times[N_{2}]\times[K]$ of $\bm{a}$.

  3. 3.

For each local minimizer $m_{s}=(i_{s},j_{s},k_{s})$, detect the corresponding
plateau, i.e. the largest connected component
$\mathcal{M}_{s}\in[N_{1}]\times[N_{2}]\times[K]$ with
$m_{s}\in\mathcal{M}_{s}$ such that

$\displaystyle a_{i,j,k}\leq ra_{i_{s},j_{s},k_{s}}\quad\text{for all
}(i,j,k)\in\mathcal{M}_{s}.$

Here, $r\in(0,1)$ is a relative threshold that determines the size of the
resulting regions. If we choose $r$ closer to 1, the resulting regions will be
tighter but the results of the method will be less robust against numerical
errors.

###  4.4 Visualization

The result of the extraction step are regions
$\mathcal{M}_{1},\ldots,\mathcal{M}_{S}\subset[N_{1}]\times[N_{2}]\times[K]$.
Since these are sets in a discrete three-dimensional space, they are difficult
to visualize directly. The extent of the regions along the first two axes (the
spatial domain) expresses uncertainty with respect to position of the
corresponding blob, while the extent along the third axis (the scale)
corresponds to uncertainty in the size of the corresponding blob (see the
middle panel in Figure 7 below).

We suggest to visualize the uncertainty in scale and position by introducing
two projections on the pixel grid $[N_{1}]\times[N_{2}]$ which are easy to
visualize as images.

The first projection is the direct projection on the image domain,

$\displaystyle\Pi_{1}:2^{[N_{1}]\times[N_{2}]\times[K]}\to
2^{[N_{1}]\times[N_{2}]},$ (50)
$\displaystyle\Pi_{1}(\mathcal{M}):=\left\\{(i,j)~{}:~{}\exists
k\in[K]:(i,j,k)\in\mathcal{M}\right\\}.$

The set $\Pi_{1}(\mathcal{M})$ is a two-dimensional region that contains the
centers of all blobs in $\mathcal{M}$, and can thus be used to visualize the
uncertainty in the center position of the uncertain blob.

The second projection is motivated by the visualization for the discrete
Laplacians-of-Gaussians method, where a point
$(i,j,k)\in[N_{1}]\times[N_{2}]\times[K]$ is visualized by a (discretized)
circle $B_{\sqrt{2t_{k}}}(i,j)$. Taking the possibly different grid sizes
$h_{1}$ and $h_{2}$ into account, the set
$B_{r}(i,j)\subset[N_{1}]\times[N_{2}]\times[K]$ is defined as

$\displaystyle
B_{r}(i,j)=\left\\{(i^{\prime},j^{\prime})~{}:~{}\sqrt{h_{1}^{2}(i^{\prime}-i)^{2}+h_{2}^{2}(j^{\prime}-j)^{2}}\leq
r\right\\}.$

We then define the projection of a set
$\mathcal{M}\subset[N_{1}]\times[N_{2}]\times[K]$ onto the union over all
circles corresponding to points in $\mathcal{M}$ by

$\displaystyle\Pi_{2}:2^{[N_{1}]\times[N_{2}]\times[K]}\to
2^{[N_{1}]\times[N_{2}]},$ (51)
$\displaystyle\Pi_{2}(\mathcal{M}):=\bigcup_{(i,j,k)\in\mathcal{M}}B_{\sqrt{2t_{k}}}(i,j).$

Together, $\Pi_{1}(\mathcal{M})$ and $\Pi_{2}(\mathcal{M})$ allow to visualize
the uncertainty in the blob center and the blob scale within a set
$\mathcal{M}$. An example of this visualization is provided in Figure 7
(middle panel). Figure 6 gives an example how this visualization can be used
on one-dimensional signals.

###  4.5 Further remarks

As mentioned in the beginning of Section 4.2, the TV-ULoG optimization problem
(21) has many parallels to CTV-denoising. For example, Weiss et al. (2009)
considers constrained-TV regularization and proposes a Nesterov smoothing
approach similar to the one we describe in Section 4.2.2. The TV-ULoG
optimization problem could also be attacked with ADMM, similar to the method
presented in Chan et al. (2013).

However, there are two significant differences between the TV-ULoG
optimization problem and the CTV-denoising problem (22), apart from the higher
dimensionality: First, the CTV-denoising problem depends through the TV-term
on the first derivative of the image, while in the TV-ULoG optimization
problem we take the total variation of the normalized Laplacian, which means
that we depend on the third-order derivatives. The resulting discrete
differential operator ${\bm{A}}$ has a high condition number, which makes the
TV-ULoG optimization problem more difficult to solve numerically. This problem
is amplified by the fact that there is no denoising term in the TV-ULoG
optimization problem (21), since this term otherwise has a regularizing
effect. These differences could explain the problems of the proposed first-
order methods we observed in our numerical experiments (see Section 5.3).

Alternatively, problem (21) could be solved with a semi-smooth Newton method
Hintermüller et al. (2002), which is likely more accurate than the first-order
methods. Such an approach would scale similar as the interior-point strategy.

A final option that we did not investigate further are graph-cut methods
Chambolle (2005); Darbon and Sigelle (2006), which could potentially be very
efficient for our approach, since the result of the TV-ULoG method should be
robust against the quantization error. However, such a method would be more
difficult to implement.

##  5 Numerical experiments

In this section, we illustrate the tube-based blob detection methods
introduced in this paper on two Bayesian inverse problems.

We start with a one-dimensional problem which mostly serves didactical
purposes, since the scale-space representation of a one-dimensional signal is
two-dimensional and can therefore be visualized as an image. This allows us to
further illustrate and discuss the ideas underlying our approach.

As second example we consider a more challenging two-dimensional imaging
problem from stellar dynamics since this application was a main motivation for
the present work.

The numerical experiments were implemented in Python 3.10. The source code is
available from the GitHub repository github.com/FabianKP/tvulog, which also
provides an exact list of the used packages. The reported computation times
were measured on a PC with 64 GB RAM and 8 3.6-GHz Intel i7-7700 CPUs.

###  5.1 One-dimensional deconvolution

####  5.1.1 Problem setup

Figure 2: Setup of the deconvolution problem. Left panel: The ground truth
${\bm{f}}^{*}$. Middle panel: The convolved ground truth $\bm{G}{\bm{f}}^{*}$.
Right panel: The noisy data $\bm{y}=\bm{G}{\bm{f}}^{*}+\bm{w}^{*}$.

We tested the proposed method on the problem of blob detection in one-
dimensional Bayesian deconvolution. That is, we considered the task of
identifying blobs in a one-dimensional discrete signal, modelled as
realization of a random vector ${\bm{F}}$, using data ${\bm{Y}}$ which is
given by the noisy convolution

$\displaystyle{\bm{Y}}=\bm{G}{\bm{F}}+\bm{W},$

where $\bm{G}:\mathbb{R}^{N}\to\mathbb{R}^{N}$ is a convolution operator, and
$\bm{W}$ is zero-mean Gaussian noise. While we described the proposed
methodology only for the case of two-dimensional signals (images), the
adaptation to one-dimensional signals is completely straightforward, and we
skip it to avoid repetition.

To define a full Bayesian model for the deconvolution problem, we assign a
zero-mean Gaussian prior on ${\bm{F}}$,

$\displaystyle p_{\bm{F}}=\mathcal{N}(\bm{0},\bm{\Sigma}),$ (52)

where $\bm{\Sigma}\in\mathbb{R}^{N\times N}$ is the prior covariance matrix.
In our experiments, we used a prior covariance $\bm{\Sigma}$ corresponding to
a Gaussian random Markov field prior.

Assuming that the noise $\bm{W}$ is uncorrelated with constant standard
deviation $\gamma$, we arrive at the likelihood

$\displaystyle
p_{{\bm{Y}}|{\bm{F}}}(\cdot{~{}|~{}}{\bm{f}})=\mathcal{N}(\bm{G}{\bm{f}},\gamma^{2}{\bm{I}}),\quad{\bm{f}}\in\mathbb{R}^{N}.$
(53)

Combining (52) and (53) via Bayes theorem yields the posterior density

$\displaystyle
p_{{\bm{F}}|{\bm{Y}}}({\bm{f}}{~{}|~{}}{\bm{y}})\propto\exp\left(-\frac{1}{2\gamma^{2}}\left\|\bm{G}{\bm{f}}-{\bm{y}}\right\|^{2}-\frac{1}{2}\left\|\bm{\Sigma}^{-1/2}{\bm{f}}\right\|^{2}\right).$

####  5.1.2 Simulation

Figure 3: The lower and upper bound of the scale-space tube
$[{\bm{u}}^{\mathrm{low}},{\bm{u}}^{\mathrm{upp}}]$ for the one-dimensional
deconvolution problem. The scale-space representation of the posterior mean is
plotted in the right panel for comparison. Figure 4: Slice through the scale-
space tube $[{\bm{u}}^{\mathrm{low}},{\bm{u}}^{\mathrm{upp}}]$ for a fixed
scale.

For our numerical experiment, we used a sinusoidal ground truth
${\bm{f}}^{*}\in\mathbb{R}^{N\times N}$ (see Figure 2) with $N=200$, from
which we simulated noisy data $\bm{y}=\bm{G}{\bm{f}}^{*}+\bm{w}$, where
$\bm{w}$ was generated from $\mathcal{N}(\bm{0},\gamma^{2}\bm{I})$ with
$\gamma=0.03$.

We then generated 10 000 MCMC samples using the Linear-RTO sampler provided by
the CUQIpy Python package Riis et al. (2023). Using these MCMC samples, we
computed a scale space tube
$[{\bm{u}}^{\mathrm{low}},{\bm{u}}^{\mathrm{upp}}]$ for the uncertain scale-
space representation ${\bm{U}}=\bm{\Phi}{\bm{F}}$ using the heuristic method
described in Appendix A for the crediblity parameter $\alpha=0.05$
(corresponding to 95% credibility). We used exponentially increasing scales
(cf. (15))

$\displaystyle t_{k}$ $\displaystyle=b^{k-1}t_{\text{min}},\qquad k\in[K],$
(54) $\displaystyle b$
$\displaystyle=\left(\frac{t_{\text{max}}}{t_{\text{min}}}\right)^{\frac{1}{K-1}},$

with $K=30$, $t_{\text{min}}=1$ and $t_{\text{max}}=70^{2}$. The two-
dimensional objects ${\bm{u}}^{\mathrm{low}}$ and ${\bm{u}}^{\mathrm{upp}}$
are visualized in Figure 3. Since it is hard to see the difference between the
lower and uppper bound with the naked eye, we have plotted a horizontal slice
(that is, a slice for a fixed scale) in Figure 4. For reference, we also
computed a point estimate for the signal of interest ${\bm{f}}$ in form of the
posterior mean, given by

$\displaystyle{\bm{f}}^{\text{mean}}=\frac{1}{S}\sum_{s=1}^{S}{\bm{f}}^{(s)},$

where $({\bm{f}}^{(s)})_{s=1}^{S}$ are the computed MCMC samples. We denote
the scale-space representation of ${\bm{f}}^{\text{mean}}$ by
$\bm{u}^{\text{mean}}$.

####  5.1.3 Results

Figure 5: The result of the TV-approach for the one-dimensional deconvolution
problem. First row, left panel: the scale-normalized Laplacian of the
posterior mean; middle panel: the scale-normalized Laplacian of the minimizer
of (21); right panel: the blob regions in scale space extracted from
$\tilde{\Delta}\bar{u}$. Second row: Corresponding horizontal slice for fixed
scale $t_{k}$, $k=15$. Figure 6: Plot of the posterior mean
${\bm{f}}^{\text{mean}}$ for the one-dimensional deconvolution problem
together with the uncertain blobs. The projected blob sets are visualized by
blue horizontal bars, where the solid bar indicates the center projection and
the dotted bar indicates the scale projection. The ground truth ${\bm{f}}^{*}$
(dotted line) is also plotted for reference.

To solve the optimization problem (21), we used the interior-point approach
since it was by far the most efficient (see Section 5.3 below). The normalized
Laplacian of the computed minimizer is plotted in Figure 5 (middle panel).
Compared to the normalized Laplacian of the posterior mean (left panel), the
scale-normalized Laplacian is approximately piecewise constant and attains
local minima on four clearly separated regions, which were extracted using the
thresholding procedure described in Section 4.3 (right panel). Since these
regions are difficult to make out with the bare eye, a horizontal slice (that
is, for fixed scale) is plotted in the second row of Figure 5.

In Figure 6, we plot the extracted blob regions using the procedure described
in Section 4.4, together with the posterior mean and the ground truth for
comparison. The horizontal solid bars are obtained from the projection
$\Pi_{1}$ (50). That is, they indicate the intervals in which we expect the
blob _centers_ to lie with 95%-confidence. Similarly, the dotted bars are
obtained from the projection $\Pi_{2}$ (51) and indicate the maximal extent of
the uncertain blob.

###  5.2 Integrated-light stellar population recovery

####  5.2.1 Problem setup

Next, we revisit the problem of integrated-light stellar population recovery.
This is a constrained linear imaging inverse problem of the form

$\displaystyle{\bm{Y}}=\bm{G}{\bm{F}}+\bm{W},\qquad\bm{F}\geq 0,$

where ${\bm{F}}$ is a two-dimensional non-negative density function (modelled
as random image), ${\bm{Y}}$ is a measured light spectrum and $\bm{W}$ is
zero-mean uncorrelated Gaussian noise. The observation operator
$\bm{G}:\mathbb{R}^{N_{1}\times N_{2}}\to\mathbb{R}^{M}$ is the discretization
of an integral operator that models how the density $\bm{F}$ influences the
spectrum. We do not discuss the details of this problem and the Bayesian
modelling and instead refer the reader to the previous work Parzer et al.
(2023).

####  5.2.2 Simulation

For our numerical experiment, we simulated a realization ${\bm{y}}$ of the
noisy spectrum ${\bm{Y}}$ from a ground truth ${\bm{F}}={\bm{f}}^{*}$ (Figure
7, top panel) and generated 10 000 posterior samples (after 5000 burn-in
iterations) using the SVD-MCMC method described in (Parzer et al., 2023,
section 4). As in the one-dimensional deconvolution example, we computed a
scale-space tube $[{\bm{u}}^{\mathrm{low}},{\bm{u}}^{\mathrm{upp}}]$ using the
method described in Appendix A for the credibility parameter $\alpha=0.05$ and
discrete scales given by (54) with $K=16$, $t_{\text{min}}=1$ and
$t_{\text{max}}=30^{2}$.

For reference, we also computed a point estimate for the signal of interest
${\bm{F}}$, this time in form of the maximum-a-posteriori estimate
${\bm{f}}^{\text{MAP}}$.

####  5.2.3 Results

Figure 7: Results of ULoG and TV-ULoG for integrated-light stellar population
recovery. Top panel: Ground truth ${\bm{f}}^{*}$ from which the mock data was
generated. Middle panel: MAP estimate ${\bm{f}}^{\text{MAP}}$ together with
projection of blob sets. The two projections $\Pi_{1}$ and $\Pi_{2}$ are
indicated by solid and dashed blue lines, respectively. Bottom panel: Results
of ULoG superimposed on the MAP estimate.

We computed a minimizer of the associated optimization problem (21) using the
interior-point method, since it was again the most efficient (see also Section
5.3). Two blob regions were extracted and visualized using the procedures
described in Section 4.3 and Section 4.4. In Figure 7 (middle panel), the
projected blob regions are plotted together with the MAP estimate
${\bm{f}}^{\text{MAP}}$.

####  5.2.4 Comparison with ULoG

The bottom panel of Figure 7 shows the result of the ULoG method (see Section
3.3). Recall that in the ULoG method, we compute a minimizer $\bar{\bm{u}}$ of

$\displaystyle\min_{\bm{u}\in\mathbb{R}^{N_{1}\times N_{2}\times K}}$
$\displaystyle\left\|\tilde{\bm{\Delta}}\bm{u}\right\|^{2}$ s.t.
$\displaystyle{\bm{u}}^{\mathrm{low}}\leq\bm{u}\leq{\bm{u}}^{\mathrm{upp}},$

and then determine the local minimum points of $\tilde{\Delta}\bar{\bm{u}}$. A
local minimum point $(i,j,k)$ of $\bar{\bm{u}}$ is then visualized by a dashed
blue circle with center $(i,j)$ and radius $\sqrt{2t_{k}}$ (see also Example
1). In contrast, the proposed TV-ULoG method yields a representative
$\bar{\bm{u}}$ such that $\tilde{\Delta}\bar{\bm{u}}$ attains its local minima
on connected regions (see also Section 3.4), which are visualized in the
middle panel of Figure 7 using the method described in Section 4.4.

We see that both methods detect two blobs with confidence. However, the
projected blob regions obtained from TV-ULoG allow for a better localization
of the blob center. More importantly, the projected regions have a clear
interpretation: The inner regions contain the center of the uncertain blob
with 95%-probability. The outer regions contain the corresponding blob
circles. In contrast, the dashed circles provided by ULoG have a less clear
interpretation since they mix scale and uncertainty information (see the
discussion after (11)).

###  5.3 Comparison of optimization methods

Both for the one-dimensional deconvolution example and the stellar recovery
example, we compared the performance of the optimization strategies proposed
in Section 4.2 for the solution of the optimization problem (21), namely the
dual-smoothing approach (Section 4.2.1), the primal-smoothing approach
(Section 4.2.2) and the interior-point method (Section 4.2.3).

Both for the dual- and primal-smoothing approach we used the FGP method (see
Algorithm 1) to solve the smoothed optimization problem, while we used ECOS
Domahidi et al. (2013) for the interior-point approach.

The results of our comparison are plotted in Figure 8 and Figure 9.

Figure 8: Comparison of the computation time for the different optimization
methods for the deconvolution problem (see Section 5.1). The objective is
normalized so that the minimum is at 0 and the initial value is at 1. Figure
9: Comparison of the computation time for the different optimization methods
for the stellar recovery problem (see Section 5.2). The objective is
normalized so that the minimum is 0.

In both cases, the interior-point approach is able to achieve much higher
accuracy. For the deconvolution problem, the first-order methods take
considerably longer and are not able to achieve a high precision. For the
stellar recovery problem, the primal smoothing method converges faster to a
low-accuracy solution, but is not able to further improve.

Also, both the dual- and primal-smoothing approach have the additional
disadvantage that they depend on the choice of a smoothing parameter $\mu$,
where a larger value of $\mu$ corresponds to more smoothing but higher
approximation error (see Remark 5). In Figure 10, we plotted the performance
of the dual and primal smoothing approach on the stellar recovery problem for
different choices of $\mu$. The trade-off between speed of convergence and
achieveable accuracy is clearly visible. If the smoothing parameter is chosen
too small, the first-order methods do not converge in a practical amount of
time. This was in particular the case for the choice of smoothing parameter
suggested by the bound (36). We tested many different choices, but did not
find a configuration for which the performance of the first-order methods was
comparable to the interior-point method. Furthermore, we also tested other
solvers for the smoothed optimization problem, such as the L-BFGS-B method,
but the performance was similar to the FGP method.

Figure 10: The performance of the dual- and primal-smoothing approach for the
stellar recovery problem, for various choices of smoothing parameter $\mu$.
For each choice the plot shows the performance of 200 000 FGP iterations.

Since it is necessary to achieve very high precision in the objective function
of (21) to obtain the desired piecewise-constant solutions, the interior-point
method should be the method of choice, since it has the additional advantage
that it does not require hand-tuning of a smoothing parameter.

##  6 Conclusion

In this work, we have developed a novel approach for blob detection in
uncertain images that represents the uncertainty in the detected blobs by
regions in scale space. These regions are obtained from the minimizer of a
non-smooth optimization problem. Using similarities to CTV-denoising, we
proposed three approaches for the numerical solution of the discretized
problem. We also described how the scale space regions can be visualized on
the image domain in an interpretable way.

The proposed method was illustrated on two numerical examples – one-
dimensional deconvolution and integrated-light stellar population recovery –
where it yielded clear results that were consistent with the ground truth. We
also evaluated the performance of the different optimization methods and
observed that the interior-point method outperformed the other two approaches,
assumedly because it is more robust against the ill-conditioning of the
problem.

Our proposed method is flexible since it only requires access to a tube in
which the scale-space representation of the uncertain image lies with high
probability. Such a tube can be computed for many applications, for example in
the important special case of Bayesian imaging.

Finally, the proposed methods are not specific to astronomical applications,
although they were originally developed in that context. The methodology could
be applied in any setting were blob detection in uncertain signals is
relevant, for example in medical and geophysical imaging.

Acknowledgments

FP and OS were supported by the Austrian Science Fund (FWF), with SFB F68
“Tomography Across the Scales”, project F6807-N36 (Tomography with
Uncertainties). The financial support by the Austrian Federal Ministry for
Digital and Economic Affairs, the National Foundation for Research, Technology
and Development and the Christian Doppler Research Association is gratefully
acknowledged. The authors are grateful to Prashin Jethwa and Glenn van de Ven
for the productive cooperation within the SFB. The authors would also like to
thank Noemi Naujoks for feedback on the manuscript.

##  Appendix A Estimation of scale-space tubes from MCMC samples

In this appendix, we describe how, given an observation ${\bm{Y}}={\bm{y}}$,
we can estimate a credible scale-space tube
$[{\bm{u}}^{\mathrm{low}},{\bm{u}}^{\mathrm{upp}}]$ using MCMC samples
${\bm{f}}^{(1)},\ldots,{\bm{f}}^{(S)}$ from the posterior distribution
${\mathbb{P}}_{{\bm{F}}|{\bm{Y}}}(\cdot{~{}|~{}}{\bm{y}})$. Recall that MCMC
is a class of methods that produce approximate samples from a probability
distribution using its density function Brooks et al. (2011). Given MCMC
samples, the standard way to estimate the posterior probability that
${\bm{F}}$ is contained in a given set $A\subset\mathbb{R}^{N_{1}\times
N_{2}}$ is Meyn and Tweedie (2009)

$\displaystyle{\mathbb{P}}_{{\bm{F}}|{\bm{Y}}}(A{~{}|~{}}{\bm{y}})\approx\frac{1}{S}\sum_{s=1}^{S}\mathbbm{1}_{A}({\bm{f}}^{(s)}),$
(55)

where the indicator function $\mathbbm{1}_{A}$ is defined as

$\displaystyle\mathbbm{1}_{A}({\bm{f}}):=\begin{cases}1,&\text{if }{\bm{f}}\in
A,\\\ 0,&\text{otherwise.}\end{cases}$

For our purposes, we want to find the minimal-volume tube
$[{\bm{u}}^{\mathrm{low}},{\bm{u}}^{\mathrm{upp}}]$ such that

$\displaystyle{\mathbb{P}}_{\bm{U}|{\bm{Y}}}([{\bm{u}}^{\mathrm{low}},{\bm{u}}^{\mathrm{upp}}]{~{}|~{}}{\bm{y}})\geq
1-\alpha.$ (56)

Recall that for a set $A\subset\mathbb{R}^{N_{1}\times N_{2}\times K}$,
${\mathbb{P}}_{\bm{U}|{\bm{Y}}}(A{~{}|~{}}{\bm{y}})$ is given by the
pushforward

$\displaystyle{\mathbb{P}}_{{\bm{U}}|{\bm{Y}}}(A{~{}|~{}}{\bm{y}})={\mathbb{P}}_{{\bm{F}}|{\bm{Y}}}(\bm{\Phi}^{-1}(A){~{}|~{}}{\bm{y}}),$
(57)

where $\bm{\Phi}$ denotes the operator that maps an image ${\bm{f}}$ to its
discrete scale-space representation (cf. Section 4.1). Inserting (57) in (56)
and using the approximation (55), we obtain

$\displaystyle 1-\alpha$
$\displaystyle\leq{\mathbb{P}}_{{\bm{F}}|{\bm{Y}}}(\bm{\Phi}^{-1}([{\bm{u}}^{\mathrm{low}},{\bm{u}}^{\mathrm{upp}}]){~{}|~{}}{\bm{y}})$
$\displaystyle\approx\frac{1}{S}\sum_{s=1}^{S}\mathbbm{1}_{\bm{\Phi}^{-1}([{\bm{u}}^{\mathrm{low}},{\bm{u}}^{\mathrm{upp}}])}({\bm{f}}^{(s)})$
$\displaystyle=\frac{1}{S}\left|\left\\{s\in[S]~{}:~{}\bm{\Phi}({\bm{f}}^{(s)})\in[{\bm{u}}^{\mathrm{low}},{\bm{u}}^{\mathrm{upp}}]\right\\}\right|,$

where for a set $A$, $\left|A\right|$ denotes its cardinality. This leads to a
condition on $[{\bm{u}}^{\mathrm{low}},{\bm{u}}^{\mathrm{upp}}]$ in terms of
samples:

$\displaystyle\left|\left\\{s\in[S]~{}:~{}\bm{\Phi}({\bm{f}}^{(s)})\in[{\bm{u}}^{\mathrm{low}},{\bm{u}}^{\mathrm{upp}}]\right\\}\right|\geq(1-\alpha)\cdot
S.$ (58)

Given samples ${\bm{f}}^{(1)},\ldots,{\bm{f}}^{(S)}$ it thus remains to find
the smallest-volume tube $[{\bm{u}}^{\mathrm{low}},{\bm{u}}^{\mathrm{upp}}]$
that satisfies (58). However, as was already discussed in (Parzer et al.,
2023, section 5.3), solving this problem exactly is computationally infeasible
even for moderately sized images. Instead, one has to use heuristic approaches
that determine a small but not minimal tube that satisfies (58).

Our method is based on the idea of sorting the samples in descending order
according to their probability, and then using bisection to find the smallest-
volume tube that satisfies (58) among the increasing sequence of tubes spanned
by the ordered samples. In detail: Let $p_{{\bm{F}}|{\bm{Y}}}$ denote the
density function of ${\mathbb{P}}_{{\bm{F}}|{\bm{Y}}}$. Using a sorting
algorithm, we find a reordering $(i_{s})_{s=1}^{S}$ such that

$\displaystyle
p_{\bm{F}|{\bm{Y}}}(\bm{f}^{(i_{1})}{~{}|~{}}{\bm{y}})\geq\ldots\geq
p_{\bm{F}|{\bm{Y}}}(\bm{f}^{(i_{S})}{~{}|~{}}{\bm{y}}).$ (59)

Then, we compute for each sample ${\bm{f}}^{(i)}$ its scale-space
representation $\bm{u}^{(i)}=\bm{\Phi}({\bm{f}}^{(i)})$. For each $s\in[S]$,
the discrete scale-space representations
$\bm{u}^{(i_{1})},\ldots,\bm{u}^{(i_{s})}$ span a tube
$[\bm{u}^{\mathrm{low},s},\bm{u}^{\mathrm{upp},s}]$ given by

$\displaystyle\bm{u}^{\mathrm{low},s}_{i,j,k}$
$\displaystyle:=\min\left\\{u_{ijk}^{(i_{r})}~{}:~{}r\in[s]\right\\},$ (60)
$\displaystyle\bm{u}^{\mathrm{upp},s}_{i,j,k}$
$\displaystyle:=\max\left\\{u_{ijk}^{(i_{r})}~{}:~{}r\in[s]\right\\}.$

(The tube $[\bm{u}^{\mathrm{low},s},\bm{u}^{\mathrm{upp},s}]$ is simply the
smallest-volume tube that contains
$\bm{u}^{(i_{1})},\ldots,\bm{u}^{(i_{s})}$). Since these tubes are
monotonically increasing, the desired tube can be found very efficiently using
bisection. The detailed pseudocode for this method is given in Algorithm 2.

Note that this algorithm is different from the method used in our previous
work (Parzer et al., 2023, section 5.3), which in particular did not use
evaluations of the posterior density function. In our numerical experiments,
the new method performed better, in the sense that it yielded tubes that
contained the same number of samples but had lower volume.

Algorithm 2 Credible scale-space tubes from samples

1:Samples ${\bm{f}}^{(1)},\ldots,{\bm{f}}^{(S)}$ from
$p_{{\bm{F}}|{\bm{Y}}}(\cdot{~{}|~{}}{\bm{y}})$, a credibility parameter
$\alpha\in(0,1)$, and a maximum number of bisection steps $M\in\mathbb{N}$.

2:Find a reordering $i_{1},\ldots,i_{S}$ such that (59) holds;

3:Compute ${\bm{u}}^{(i)}=\bm{\Phi}{\bm{f}}^{(i)}$ for all $i\in[S]$;

4:$S_{\alpha}=\lceil(1-\alpha)S\rceil$;

5:Set $T=[\bm{u}^{\mathrm{low},S_{\alpha}},\bm{u}^{\mathrm{upp},S_{\alpha}}]$
(cf. Equation 60);

6:$K=\left|\left\\{s\in[S]~{}:~{}\bm{u}^{(i_{s})}\in T\right\\}\right|$;

7:if $K=S_{\alpha}$ then return $T$

8:else

9: Set $K_{\text{min}}=1$ and $K_{\text{max}}=S_{\alpha}$;

10: for $m=1,\ldots,M$ do

11: if $K>S_{\alpha}$ then

12: $K_{\text{max}}=K$;

13: $K=\frac{1}{2}(K_{\text{min}}+K)$;

14: else if $K<S_{\alpha}$ then

15: $K_{\text{min}}=K$;

16: $K=\frac{1}{2}(K+K_{\text{max}})$;

17: else

18: break;

19: end if

20: Set $T=[{\bm{u}}^{\mathrm{low},K},{\bm{u}}^{\mathrm{upp},K}]$;

21: $K=\left|\left\\{s\in[S]~{}:~{}\bm{u}^{(i_{s})}\in T\right\\}\right|$;

22: end for

23: return $T$;

24:end if

## References

  *   * Lindeberg (1998) Lindeberg, T.: Feature Detection with Automatic Scale Selection. Int. J. Comput. Vis. 30(2), 79–116 (1998) https://doi.org/10.1023/A:1008045108935
  * Boecker et al. (2020) Boecker, A., Alfaro-Cuello, M., Neumayer, N., Martín-Navarro, I., Leaman, R.: Recovering age-metallicity distributions from integrated spectra: Validation with MUSE data of a nearby nuclear star cluster. Astrophys. J. 896(1), 13 (2020) https://doi.org/10.3847/1538-4357/ab919d
  * Parzer et al. (2023) Parzer, F., Jethwa, P., Boecker, A., Alfaro-Cuello, M., Scherzer, O., van de Ven, G.: Uncertainty-aware blob detection with an application to integrated-light stellar population recoveries. Astron. Astrophys. 674 (2023) https://doi.org/10.1051/0004-6361/202244739
  * Mikolajczyk et al. (2005) Mikolajczyk, K., Tuytelaars, T., Schmid, C., Zisserman, A., Matas, J., Schaffalitzky, F., Kadir, T., Gool, L.V.: A Comparison of Affine Region Detectors. Int. J. Comput. Vision. 65(1-2), 43–72 (2005) https://doi.org/10.1007/s11263-005-3848-x
  * Lindeberg (2013) Lindeberg, T.: Scale Selection Properties of Generalized Scale-Space Interest Point Detectors. J. Math. Imaging. Vis. 46(2), 177–210 (2013) https://doi.org/10.1007/s10851-012-0378-3
  * Li et al. (2015) Li, Y., Wang, S., Tian, Q., Ding, X.: A survey of recent advances in visual feature detection. Neurocomputing 149, 736–751 (2015) https://doi.org/10.1016/j.neucom.2014.08.003
  * Lindeberg (1994) Lindeberg, T.: Scale-Space Theory in Computer Vision. Springer, Boston (1994)
  * Sporring et al. (1997) Sporring, J., Nielsen, M., Florack, L., Johansen, P. (eds.): Gaussian Scale-Space Theory. Springer Science+Business Media, Dordrecht (1997)
  * Lindeberg (1994) Lindeberg, T.: Scale-space theory: A basic tool for analysing structures at different scales. J. Appl. Stat. 21(2), 225–270 (1994)
  * Weickert (1998) Weickert, J.: Anisotropic Diffusion in Image Processing. B. G. Teubner, Stuttgart (1998)
  * Holmström and Pasanen (2017) Holmström, L., Pasanen, L.: Statistical Scale Space Methods. Int. Stat. Rev. 85(1), 1–30 (2017) https://doi.org/10.1111/insr.12155
  * Godtliebsen et al. (2004) Godtliebsen, F., Marron, J.S., Chaudhuri, P.: Statistical significance of features in digital images. Image Vis. Comput. 22(13), 1093–1104 (2004) https://doi.org/10.1016/j.imavis.2004.05.002
  * Godtliebsen and Øigård (2005) Godtliebsen, F., Øigård, T.A.: A visual display device for significant features in complicated signals. Comput. Stat. Data Anal. 48(2), 317–343 (2005) https://doi.org/10.1016/j.csda.2004.01.010
  * Thon et al. (2012) Thon, K., Rue, H., Skrøvseth, S.O., Godtliebsen, F.: Bayesian multiscale analysis of images modeled as Gaussian Markov random fields. Comput. Stat. Data Anal. 56(1), 49–61 (2012) https://doi.org/10.1016/j.csda.2011.07.009
  * Pasanen et al. (2013) Pasanen, L., Launonen, I., Holmström, L.: A scale space multiresolution method for extraction of time series features. Stat 2(1), 273–291 (2013) https://doi.org/10.1002/sta4.35
  * Repetti et al. (2019) Repetti, A., Pereyra, M., Wiaux, Y.: Scalable Bayesian Uncertainty Quantification in Imaging Inverse Problems via Convex Optimization. SIAM J. Imaging Sci. 12(1), 87–118 (2019) https://doi.org/10.1137/18M1173629
  * Price et al. (2021) Price, M.A., McEwen, J.D., Cai, X., Kitching, T.D., Wallis, C.G.R., the LSST Dark Energy Science Collaboration: Sparse Bayesian mass mapping with uncertainties: Hypothesis testing of structure. Mon. Not. R. Astron. Soc. 506(3), 3678–3690 (2021) https://doi.org/10.1093/mnras/stab1983
  * Price et al. (2020) Price, M.A., McEwen, J.D., Pratley, L., Kitching, T.D.: Spherical Bayesian mass-mapping with uncertainties: Full sky observations on the celestial sphere. Mon. Not. R. Astron. Soc. 500(4), 5436–5452 (2020) https://doi.org/10.1093/mnras/staa3563
  * Pereyra (2017) Pereyra, M.: Maximum-a-Posteriori Estimation with Bayesian Confidence Regions. SIAM J. Imaging Sci. 10(1), 285–302 (2017) https://doi.org/10.1137/16M1071249
  * Kallenberg (2002) Kallenberg, O.: Foundations of Modern Probability, 2nd edn. Springer-Verlag, New York (2002). https://doi.org/10.1007/978-1-4757-4015-8
  * Schervish (1995) Schervish, M.J.: Theory of Statistics. Springer, New York (1995)
  * Alvarez et al. (1993) Alvarez, L., Guichard, F., Lions, P.-L., Morel, J.-M.: Axioms and fundamental equations of image processing. Arch. Rational Mech. Anal. 123(3), 199–257 (1993) https://doi.org/10.1007/BF00375127
  * Florack et al. (1992) Florack, L.M.J., ter Haar Romeny, B.M., Koenderink, J.J., Viergever, M.A.: Scale and the differential structure of images. Image Vis. Comput. 10(6), 376–388 (1992) https://doi.org/10.1016/0262-8856(92)90024-W
  * Duits et al. (2003) Duits, R., Felsberg, M., Florack, L., Platel, B.: $\alpha$ Scale Spaces on a Bounded Domain. In: Griffin, L.D., Lillholm, M. (eds.) Scale Space Methods in Computer Vision, pp. 494–510. Springer, Berlin, Heidelberg (2003). https://doi.org/10.1007/3-540-44935-3_34
  * Brezis (2011) Brezis, H.: Functional Analysis, Sobolev Spaces and Partial Differential Equations. Springer, New York, NY (2011). https://doi.org/10.1007/978-0-387-70914-7
  * Marr and Hildreth (1980) Marr, D., Hildreth, E.: Theory of edge detection. Proc. Roy. Soc. London, Ser. B, Biol. Sci. 207(1167), 187–217 (1980) https://doi.org/10.1098/rspb.1980.0020
  * Voorhees and Poggio (1988) Voorhees, H., Poggio, T.: Computing texture boundaries from images. Nature 333(6171), 364–367 (1988) https://doi.org/10.1038/333364a0
  * Geman and Geman (1986) Geman, D., Geman, S.: Bayesian Image Analysis. In: Bienenstock, E., Soulié, F.F., Weisbuch, G. (eds.) Disordered Systems and Biological Organization, pp. 301–319. Springer, Berlin, Heidelberg (1986). https://doi.org/10.1007/978-3-642-82657-3_30
  * Hanson (1993) Hanson, K.M.: Introduction to Bayesian image analysis. In: Loew, M.H. (ed.) Medical Imaging 1993: Image Processing, vol. 1898. Newport Beach, CA, pp. 716–731 (1993). https://doi.org/10.1117/12.154577
  * Klenke (2014) Klenke, A.: Probability Theory: A Comprehensive Course. Springer, London (2014). https://doi.org/10.1007/978-1-4471-5361-0
  * Kaipio and Somersalo (2005) Kaipio, J., Somersalo, E.: Statistical and Computational Inverse Problems. Springer, New York (2005)
  * Dashti and Stuart (2017) Dashti, M., Stuart, A.M.: The Bayesian Approach to Inverse Problems. In: Ghanem, R., Higdon, D., Owhadi, H. (eds.) Handbook of Uncertainty Quantification, pp. 311–428. Springer International Publishing, Cham (2017). https://doi.org/10.1007/978-3-319-12385-1_7
  * Brooks et al. (2011) Brooks, S.P., Gelman, A., Jones, G.L., Meng, X.-L. (eds.): Handbook of Markov Chain Monte Carlo. Chapman and Hall/CRC, New York (2011). https://doi.org/10.1201/b10905
  * Blei et al. (2017) Blei, D.M., Kucukelbir, A., McAuliffe, J.D.: Variational Inference: A Review for Statisticians. J. Am. Stat. Assoc. 112(518), 859–877 (2017) https://doi.org/10.1080/01621459.2017.1285773
  * Wasserman (2006) Wasserman, L.: All of Nonparametric Statistics. Springer, New York, NY (2006). https://doi.org/10.1007/0-387-30623-4
  * Mammen (1995) Mammen, E.: On Qualitative Smoothness of Kernel Density Estimates. Stat. 26 (1995) https://doi.org/10.1080/02331889508802494
  * Mammen and van de Geer (1997) Mammen, E., van de Geer, S.: Locally adaptive regression splines. Ann. Stat. 25(1), 387–413 (1997) https://doi.org/10.1214/aos/1034276635
  * Davies and Kovac (2004) Davies, P.L., Kovac, A.: Densities, spectral densities and modality. Ann. Stat. 32(3), 1093–1136 (2004) https://doi.org/10.1214/009053604000000364
  * Obereder et al. (2007) Obereder, A., Scherzer, O., Kovac, A.: Bivariate density estimation using BV regularisation. Comput. Stat. Data Anal. 51(12), 5622–5634 (2007) https://doi.org/%****␣ParKirSch23_preprint.bbl␣Line␣625␣****10.1016/j.csda.2007.04.019
  * Hartigan and Hartigan (1985) Hartigan, J.A., Hartigan, P.M.: The Dip Test of Unimodality. Ann. Stat. 13(1), 70–84 (1985) https://doi.org/10.1214/aos/1176346577
  * Davies (1995) Davies, P.L.: Data features. Stat. Neerl. 49(2), 185–245 (1995) https://doi.org/10.1111/j.1467-9574.1995.tb01464.x
  * Davies and Kovac (2001) Davies, P.L., Kovac, A.: Local Extremes, Runs, Strings and Multiresolution. Ann. Stat. 29(1) (2001) https://doi.org/10.1214/aos/996986501
  * Caselles et al. (2015) Caselles, V., Chambolle, A., Novaga, M.: Total Variation in Imaging. In: Scherzer, O. (ed.) Handbook of Mathematical Methods in Imaging, pp. 1455–1499. Springer, New York, NY (2015). https://doi.org/10.1007/978-1-4939-0790-8_23
  * Davies et al. (2009) Davies, P.L., Kovac, A., Meise, M.: Nonparametric regression, confidence regions and regularization. Ann. Stat. 37(5B), 2597–2625 (2009) https://doi.org/10.1214/07-AOS575
  * Rudin et al. (1992) Rudin, L.I., Osher, S., Fatemi, E.: Nonlinear total variation based noise removal algorithms. Phys. D: Nonlinear Phenom. 60(1-4), 259–268 (1992) https://doi.org/10.1016/0167-2789(92)90242-F
  * Bredies and Lorenz (2018) Bredies, K., Lorenz, D.: Mathematical Image Processing, 1st edn. Birkhäuser, New York, NY (2018)
  * Beck and Teboulle (2009) Beck, A., Teboulle, M.: Fast Gradient-Based Algorithms for Constrained Total Variation Image Denoising and Deblurring Problems. IEEE Trans. Image Process. 18(11), 2419–2434 (2009) https://doi.org/10.1109/TIP.2009.2028250
  * Goldfarb and Yin (2005) Goldfarb, D., Yin, W.: Second-order Cone Programming Methods for Total Variation-Based Image Restoration. SIAM J. Sci. Comput. 27(2), 622–645 (2005) https://doi.org/10.1137/040608982
  * Yin et al. (2007) Yin, W., Goldfarb, D., Osher, S.: A comparison of three total variation based texture extraction models. J. Visual Commun. Image Represent. 18(3), 240–252 (2007) https://doi.org/10.1016/j.jvcir.2007.01.004
  * Bot et al. (2009) Bot, R.I., Grad, S.-M., Wanka, G.: Duality in Vector Optimization. Springer, Berlin, Heidelberg (2009). https://doi.org/10.1007/978-3-642-02886-1
  * Beck (2017) Beck, A.: First-Order Methods in Optimization. SIAM, Philadelphia (2017)
  * Nesterov (2005) Nesterov, Y.: Smooth minimization of non-smooth functions. Math. Program. 103(1), 127–152 (2005) https://doi.org/10.1007/s10107-004-0552-5
  * Huber (1981) Huber, P.J.: Robust Statistics. John Wiley & Sons, Inc, New York (1981)
  * Byrd et al. (1995) Byrd, R.H., Lu, P., Nocedal, J., Zhu, C.: A Limited Memory Algorithm for Bound Constrained Optimization. SIAM J. Sci. Comput. 16(5), 1190–1208 (1995)
  * Zhu et al. (1997) Zhu, C., Byrd, R.H., Lu, P., Nocedal, J.: Algorithm 778: L-BFGS-B: FORTRAN subroutines for large-scale bound-constrained optimization. ACM Trans. Math. Softw. 23(4), 550–560 (1997) https://doi.org/10.1145/279232.279236
  * Morales and Nocedal (2011) Morales, J.L., Nocedal, J.: Remark on Algorithm 778: L-BFGS-B, FORTRAN subroutines for large-scale bound constrained optimization. ACM Trans. Math. Softw. 38(1), 7–174 (2011) https://doi.org/10.1145/2049662.2049669
  * Boyd and Vandenberghe (2004) Boyd, S.P., Vandenberghe, L.: Convex Optimization. Cambridge University Press, Cambridge, UK ; New York (2004)
  * Weiss et al. (2009) Weiss, P., Blanc-Féraud, L., Aubert, G.: Efficient Schemes for Total Variation Minimization Under Constraints in Image Processing. SIAM J. Sci. Comput. 31(3), 2047–2080 (2009) https://doi.org/10.1137/070696143
  * Chan et al. (2013) Chan, R.H., Tao, M., Yuan, X.: Constrained Total Variation Deblurring Models and Fast Algorithms Based on Alternating Direction Method of Multipliers. SIAM J. Imaging Sci. 6(1), 680–697 (2013) https://doi.org/10.1137/110860185
  * Hintermüller et al. (2002) Hintermüller, M., Ito, K., Kunisch, K.: The Primal-Dual Active Set Strategy as a Semismooth Newton Method. SIAM J. Optim. 13(3), 865–888 (2002) https://doi.org/10.1137/S1052623401383558
  * Chambolle (2005) Chambolle, A.: Total Variation Minimization and a Class of Binary MRF Models. In: Rangarajan, A., Vemuri, B., Yuille, A.L. (eds.) Energy Minimization Methods in Computer Vision and Pattern Recognition, pp. 136–152. Springer, Berlin, Heidelberg (2005). https://doi.org/10.1007/11585978_10
  * Darbon and Sigelle (2006) Darbon, J., Sigelle, M.: Image Restoration with Discrete Constrained Total Variation Part I: Fast and Exact Optimization. J. Math. Imaging. Vis. 26(3), 261–276 (2006) https://doi.org/10.1007/s10851-006-8803-0
  * Riis et al. (2023) Riis, N.A.B., Alghamdi, A.M.A., Uribe, F., Christensen, S.L., Afkham, B.M., Hansen, P.C., Jørgensen, J.S.: CUQIpy – Part I: Computational Uncertainty Quantification for Inverse Problems in Python. arXiv (2023). https://doi.org/10.48550/arXiv.2305.16949
  * Domahidi et al. (2013) Domahidi, A., Chu, E., Boyd, S.: ECOS: An SOCP solver for embedded systems. In: 2013 European Control Conference (ECC), pp. 3071–3076. IEEE, Zurich (2013). https://doi.org/10.23919/ECC.2013.6669541
  * Meyn and Tweedie (2009) Meyn, S., Tweedie, R.L.: Markov Chains and Stochastic Stability, 2nd edn. Cambridge University Press, Cambridge (2009). https://doi.org/10.1017/CBO9780511626630