The goal of image reconstruction is obtaining an image of the object of interest from indirectly measured, and typically noisy, data. Mathematically, image reconstruction problems are commonly formulated as inverse problems that can be written in the form of operator equations

where is the unknown, is the measurements and is a forward operator that models the data acquisition. In this paper, we consider linear forward operators and assume that equation (1.1) with exact data and operator has a unique solution that we denote by ¯u.

In practice, not only the right-hand side f is noisy, but also the operator A is often not exact as it contains errors that come from imperfect calibration measurements.

Uncertainty in the operator and in the data may be characterised by the inclusions and for some sets . These sets may be referred to as uncertainty sets – a concept widely used in robust optimisation [4]. Given A and F, we would like to find a subset , called the feasible set, that contains the exact solution ¯u. Depending on the particular form of the uncertainty sets A and F and on available additional a priori information about ¯u, the inclusion ¯can be proven for different feasible sets. Two main considerations that affect the choice of a particular feasible set are its size (smaller feasible sets are preferred) and the availability of efficient numerical algorithms for optimisation problems involving the feasible set (therefore, convex feasible sets are preferred).

For ill-posed problems, in general, available feasible sets are too large and contain elements arbitrary far from ¯u. An exception is the case when a compact set that contains the exact solution ¯u is known a priori [24]. In this case, a feasible set of finite diameter can be obtained. In the general case, an appropriate regularisation functional R needs to be minimised on the feasible set to find a stable approximation to ¯u.

This is the idea behind the residual method [14, 13]. Operating in normed spaces, one can define the uncertainty sets as follows:

for an approximate right-hand side , approximate forward operator and approximation errors [14]. Using the information in (1.2), one can define a feasible set as follows [14]:

This set contains all elements of U that are consistent with (1.1) within the tolerances given by (1.2). The inclusion ¯can be easily verified. Unless h = 0 (i.e. the forward operator is exact), the set is non-convex and the residual method results in a non-convex optimisation problem even for convex regularisation functionals.

An alternative approach to modelling uncertainty in A and f using partially ordered spaces was proposed in [15, 16]. Assume that U and F are Banach lattices, i.e. Banach spaces with partial order “”, and that A is a regular operator [21]. Then, uncertainties in A and f can be characterised using intervals in appropriate partial orders, i.e.

where ) is the space of regular operators . Assuming positivity of the exact solution ¯u and using the inequalities in (1.4), we can show that the exact solution ¯u is contained in the following feasible set [15]:

In contrast to (1.3), the set in (1.5) is convex and minimising a convex regularisation functional R on this set results in a convex optimisation problem:

Using the relationship between partial orders and norms in Banach lattices, one can prove the inclusion of the partial-order-based feasible set U (1.5) in the norm-based feasible set for appropriately chosen We briefly review the partial-order-based approach in Section 2.

While convergence of the minimisers of (1.6) to the exact solution can be guaranteed [15], it is not clear, whether the solution to (1.6) actually corresponds to a particular pair (A, f) within the bounds (1.4). It seems more natural to look for approximate solutions in the following set:

i.e. to assume that, while the exact operator and noise-free measurements are unknown, there has to exist at least one pair within the uncertainty bounds (1.4) that exactly explains the solution.

It is not clear a priori, whether the set (1.7) is convex. In this paper we show that the sets (1.7), in fact, coincide, which implies convexity of (see Section 3) and shows that the convex problem (1.6) actually implements the natural formulation (1.7).

It is tempting to add an additional constraint on the operator in (1.7), reflecting additional a priori information about A. For instance, if A is a blurring operator, then (after finite-dimensional approximation) the rows of the matrix A should sum up to one, i.e. we should have that Ae = e for two vectors of ones of appropriate lengths. More generally, one can define an additional linear constraint Av = g for a fixed pair (v, g) to obtain

Unfortunately, even such a simple constraint breaks the convexity of the feasible set. We demonstrate this in Section 4 by explicitly describing the set (1.8) in finite dimensions in the special case when . We also argue that the additional constraint Av = g can still be useful to tighten the bounds if they weren’t carefully chosen initially.

Since the set {A: Av = g} is convex (in A), the analysis of Section 4 shows that convexity of an additional constraint set does not guarantee convexity of the corresponding feasible set in u. Because of this negative result, we confine ourselves to the set (1.5) in our numerical experiments.

In Section 5 we consider an application in image deblurring with uncertainty in the blurring kernel. In many applications, such as astronomy or fluorescence microscopy, the blurring kernel (often referred to as the point spread function) is obtained experimentally by recording light from reference stars [2] or imaging subresolution fluorescent particles [22]. Such blurring kernels inevitably contain errors that can significantly impact the reconstruction. Blind deconvolution [8, 17] aims at reconstructing both the blurring kernel and the image simultaneously, but suffers from severe ill-posedness and non-convexity. The approach we propose takes into account the errors in the available blurring kernel (without attempting to obtain a better estimate of it) while staying within the convex setting.

spaces, endowed with a partial order relation

become Banach lattices, i.e. partially ordered Banach spaces with well-defined suprema and infima of each pair of elements and a monotone norm [21]. If E and F are two Banach lattices, then partial orders in E and F induce a partial order in a subspace of the space of linear operators acting from E to F, namely in the space of regular operators. A linear operator is called regular, if it can be represented as a difference of two positive operators. Positivity of an operator is defined as 0. Partial order in the space of regular operators is introduced as follows: is a positive operator. Every regular operator acting between two Banach lattices is continuous [21].

The framework of partially ordered functional spaces allows quantifying uncertainty in the data f and forward operator A of the inverse problem (1.1) by means of order intervals (1.4). Approximate solutions to (1.1) are the minimisers of (1.6). Convergence of these minimisers to the exact solution of (1.1) is studied as the uncertainty in the data f and forward operator A diminishes. This is formalised using monotone convergence sequences of lower and upper bounds (1.4):

If a sequence of lower (upper) bounds is not monotone, it can always be made monotone by consequently taking the supremum (infimum) of each element in the sequence with the preceding one.

Convergence of the corresponding sequence of minimisers of (1.6) to ¯u is guaranteed by the following theorem [15]:

Theorem 2.1. Let U and F be Banach lattices and F order completeSuppose that the regulariser R satisfies the following assumptions:

is bounded from below on U;

is lower semi-continuous;

• non-empty sub-level sets are strongly sequentially compact.

Then strongly in U.

The assumptions on the regulariser in Theorem 2.1 are rather standard. Conditions of Theorem 2.1 are satisfied, for example, if ). The term can be dropped if boundedness of the -norm can be guaranteed for all (1.5). The term TV (u) can be replaced by any topologically equivalent seminorm, such as ) [5] (see [6] for a proof of topological equivalence) or

Strong compactness of the sub-level sets of R can be replaced by weak compactness if R has the Radon-Riesz property, i.e. that for any sequence weak convergence with convergence of the values ) implies strong convergence . With this modification, Theorem 2.1 admits norms in reflexive Banach spaces as regularisers.

The constraint 0 in (1.5) is important. It can be relaxed to necessarily 0), with some modifications in the formulae [16], provided that the exact solution ¯u satisfies this constraint. If the exact solution may be unbounded from below, the method won’t work.

Let us briefly discuss the inclusion of the partial-order-based feasible set U (1.5) in the norm-based feasible set (1.3). Let us choose

It is easy to verify that we note that, since

and therefore

Since the space of regular operators order complete is a Banach lattice under the so-called [1] and the r-norm is always greater or equal to the operator norm [21], (2.3) implies

The inequality can be shown analogously and (0 follows from (2.1). The proof of the inclusion of (1.5) in (1.3) with as defined in (2.2) can be found in [15, Thm. 2]. Therefore, the partial-order-based feasible set (1.5) is contained in the norm-based feasible set (1.3) if the approximate operator, the approximate right-hand side and the approximation errors are as in (2.2).

Theorem 2.1 guarantees convergence of approximate solutions chosen from the partial-order-based feasible set U (1.5) by minimizing a regulariser over U as in (1.6). It is not clear, however, whether the minimisers solve (1.1) for any particular pair (A, f) within the bounds (1.4). In this section, we give a positive answer to this question for regular integral operators [1] acting between two spaces )-measurable functions, respectively, where S and T are sets, Σ-algebras over these sets and are measures. A linear operator is called an integral operator, if there exists a jointly measurable function ) such that for each

for -almost all . An integral operator is regular if and only if the operator

has range in F [1, Thm. 5.11].

Theorem 3.1. Let U and F in (1.1) be spaces of (-measurable functions, respectively. Let be regular integral operators, be as defined in (1.5). Then for every there exist a regular operator , such that Au = f.

are integral operators, there exist jointly measurable functions and ) such that

for -almost all and by [1, Thm. 5.5] we have -almost all (

is an integral operator and therefore there exists a jointly measurable function ) such that

) (a choice of a jointly measurable ) would do that), but it will suffice for our existence proof. Obviously, such choice of K(s, t) defines an integral operator A that satisfies

for -almost all -almost everywhere, the inequalities above are

or

This system has a solution if the first operand in the supremum in (3.1) is 1 and the first operand in the infimum is . This is indeed the case as a consequence of the conditions (1.4) and the inequalities . Therefore, we can always find a measurable ) satisfying (3.1), for example, by choosing

which is a supremum of two measurable functions and therefore measurable. In the special case we get a unique solution

Hence, we have found a pair (A, f) within the bounds (1.4) for an arbitrary

inclusion holds as well, since for any , with the corresponding pair (have

due to the positivity of u, and hence . Therefore, we have proven the following

Theorem 3.2. Under the assumptions of Theorem 3.1, the sets and (1.7), respectively, coincide.

An immediate consequence of this result is the convexity of the set , since the set U is, obviously, convex. An advantage of the formulation (1.5) is the ease of implementation in an optimisation algorithm. On the other hand, the formulation (1.7) allows to easily include a priori information on the operator A as additional constraints, cf. (1.8).

It is a natural question to ask, whether under some additional constraints on A the feasible set (1.7) remains convex (in u). In this section we answer this question negatively in the case when the additional constraint is linear. We restrict ourselves to the finite-dimensional case when matrix. Note that in the finite-dimensional case partial order in the space of regular operators coincides with the elementwise partial order for matrices, i.e. . Without loss of generality, we also restrict ourselves the special case

Fix a pair (

and consider the set

As noted in the introduction, the additional constraint Av = g can be useful, for example, if the exact forward operator is a convolution operator, i.e. all rows of the matrix A sum up to one. This additional constraint allows us to further restrict the feasible set, while still preserving the inclusion ¯. Intuitively, a tighter feasible set provides more information about the exact solution and can be expected to improve the reconstructions.

While the inclusion of (4.2) in the convex set U (1.5), obviously, still holds, the opposite inclusion does not hold any more. In what follows, we derive an explicit description of the set and argue that this set is not convex. Therefore, the advantages in reconstruction quality offered by using a tighter feasible set come at a price of a significant increase in computational complexity.

Every matrix , can be written as

for each row i = 1, . . . , m. In what follows we will drop the subscript i and consider this system

or, equivalently,

The matrix and the right-hand side in (4.3) have non-negative entries due to (1.5), (4.1) and the inequality . Our goal is to find conditions on u under which this system has a solution 1]. We will use Farkas’ lemma [18] to find out, when the system (4.3) has a positive solution. To find out, when it has a solution 1, we reformulate (4.3) in terms of , which gives us the following system

which also has a positive right-hand side due to (1.5) and (4.1). Combining these two systems,

find conditions under which the system (4.3) has a solution that is simultaneously The system (4.3) has a solution in [0, 1] if and only if the system (4.5) has a solution

solution. Farkas’ lemma gives us the following alternative: either (4.5) has a solution there exists a vector ) such that

and

We can rewrite these conditions equivalently as follows:

The proof will be based on considering various combinations of signs of (separately. First we prove the following

Lemma 4.1. If a solution of the system (4.8)–(4.10) exists, it satisfies the inequality (

Proof. Summing up equations (4.8) and (4.9), we get the following system:

which implies

The coefficients at () in both equations are positive. Therefore, whether , one of the above equations is violated.

Let be the complement ofJ. Inequality (4.10) requires that we choose as small as possible, which is

Substituting this into (4.10), we get

Define

Lemma 4.2. The function as defined in (4.12) has the following properties:

1. is piecewise linear;

5. is monotonically non-decreasing;

6. is continuous;

7. is convex on (0

Proof. 1.–4. Obvious.

5. Every time z crosses a point from left to right, one is replaced by a greater value , and between these points ) is constant, hence the monotonicity of

6. Suspected jumps at are zero, since the summand in (4.12) at

7. Follows from the above.

The subdifferential of ) is given by:

The minimum of ) is obtained at a point such that 0 . Let us permute

the indices so that are sorted in ascending order (the entries in must be re-sorted accordingly). This operation does not affect the products of given by the following condition:

and the minimum of

Note that, although the condition for (4.13) does not contain does depend on u due to the permutation of the indices we made.

Conditions (4.13)-(4.14) define the minimum of the function ). If this minimum is negative, then the system (4.8)-(4.10) has a solution such that and the original system (4.5) has no non-negative solution. In order for the system (4.5) to have a non-negative solution, we must have

Proceeding similarly in the case , we obtain the following conditions:

where is defined by the following condition (the indices are assumed again to be permuted in such a way that are sorted in ascending order):

consists of all (as defined in (1.7)) for which conditions (4.13), (4.15) and (4.16), (4.17) are satisfied.

Remark 4.4. If there was no dependence of , the functions ) as in (4.14) and ) as in (4.16) would be convex and the set would be a difference of the convex set U and two convex sets defined by the inequalities dependence of makes the structure of more complicated. We do not study this structure further in this work.

Figure 1 shows the set in the case when are randomly chosen in the intervals [0, respectively, f is generated using the matrix (a randomly chosen . To visualise the feasible set random points u in the positive quadrant and check conditions (4.13), (4.15) and (4.16), (4.17) (points that satisfy these conditions are shown as blue stars in Fig. 1). In this simple example, the set can be computed analytically as well, the result is shown by two solid lines in Fig. 1. As expected, the result is the same.

The result of Theorem 4.3 gives the impression that the information that A is a convolution matrix (which can be expressed as the condition Ae = e) is of little use for the approach, since including this information in the reconstruction algorithm requires solving a non-convex optimisation problem. However, the additional linear constraint can sometimes be used to tighten the bounds , if they weren’t carefully chosen initially. Indeed, one can attempt finding tighter lower and upper bounds by solving the following optimisation problems:

These optimisation problems are convex and can be efficiently solved in parallel. In the infinite-dimensional case, the analogue of the optimisation problems (4.18) is as follows:

where the inf and sup are taken in the space of regular operators . For these inf and sup to exist, the space of regular operators must be order complete, i.e. any majorised set in it must have a supremum. This is guaranteed when F is order complete [1, Theorem 1.16]. For example, the spaces of measurable functions ) are order complete, whilst the space of continuous functions C(S) is not [21].

Figure 1: The feasible set in a 2D toy problem generated using conditions (4.13), (4.15) and (4.16), (4.17) (blue stars) and using analytic formulas (between red solid lines). As expected, the two sets coincide.

Remark 4.5. An interesting question is, what is the convex hull conv . If it is smaller than U, then the feasible set for u can be tightened while preserving its convexity and better reconstructions can be expected. It is clear that in some situations conv is a strict subset of U, for example, if (cf. (4.19)). However, it is hard to say anything more about conv at the first glance. We leave the study of conv for future work.

Deblurring is widely used to improve the quality of images, for example, in astronomy [23] and fluorescence microscopy [20, 3]. Quite often, we only have an estimate of the blurring kernel, for example, when it is measured experimentally [2, 22] or obtained using simplified models [25]. Therefore, it is necessary to account for the uncertainty in the blurring kernel during the reconstruction. One possibility is to estimate the kernel and the image simultaneously, which is known as blind deblurring [8, 17]. The problem with this approach is that it results in non-convex optimisation problems and is severely ill-posed. Another option is to include the knowledge about the uncertainty in the blurring operator, if such knowledge is available, into the reconstruction process, which is the approach that we pursue.

Deblurring in 1D Let us first consider a simple one-dimensional example to get a feeling for how noise in the operator affects reconstruction and how the proposed approach can alleviate the impact of this noise. Consider the signal u shown in Fig. 2a in blue (dashed line). This signal is convolved with a Gaussian blurring kernel with standard deviation 0.5 and Dirichlet boundary conditions. Then uniform noise with support [to it (i.e. we add 1% noise). The blurred and noisy signal is shown in Fig. 2a in green (solid line). The ground-truth signal is piecewise-constant, suggesting the use of total variation [19] as the regulariser.

First we reconstruct the signal using the exact operator A by solving the optimisation problem (1.6) with . All reconstructions were computed using CVX [12, 11]. The bounds for the right-hand side can be obtained from the noisy signal f as follows: . Note that in this case the problem (1.6) is equivalent

Figure 2: Reconstruction of a piecewise constant signal using total variation. Perfect knowledge of the blurring operator yields nearly perfect reconstruction (b). However, even 10% noise in the blurring operator renders the inversion ill-posed (c). Taking uncertainty in the blurring operator into account yields stable reconstruction, but results in some loss of contrast (d).

Since TV is not strictly convex, we cannot guarantee uniqueness of the solution of (5.1). Non-uniqueness of some TV-based reconstruction models, e.g., the TV model, is a well-known issue [9]. In order to ensure uniqueness we add a small correction term to the regulariser in (5.1). In order to simplify notation we omit this correction term in the statements of optimisation problems involving TV that follow.

As expected, using the exact forward operator and a suitable regulariser we get a nearly perfect reconstruction (Fig. 2b).

Let us assume that only a slightly perturbed version ˜A of the blurring operator A is available:

where are i.i.d. uniform random numbers with support [(i.e. we add 10% noise to the operator). Note that while the true blurring matrix A is typically sparse, this is not any longer true for the perturbed matrix ˜A. It reasonable to assume, however, that all entries smaller than d are pure noise and set them to zero. We also take into account that the rows of the blurring matrix must sum up to one and normalise the rows of ˜A.

Let us solve the optimisation problem (1.6) with as defined above. This is equivalent to solving

The condition that the exact solution ¯u is in the feasible set U (1.5), which is crucial for the convergence proof (Theorem 2.1), does not hold in this case any more and therefore the regularising properties of the approach (1.6) can not be guaranteed. Indeed, we observe that an error in the operator of 10% renders the inversion ill-posed and the reconstruction highly oscillatory (Fig. 2c). This problem can be dealt with by increasing the allowed noise level in the problem (5.3), e.g., by multiplying the right-hand side of (5.3) by a factor C > 1, however, the value of C that will be sufficient to remove oscillations depends on the noise in the operator and is not straightforward to determine.

Let us now acknowledge the fact that the operator ˜A contains errors and derive lower and upper bounds for the unknown true operator from (5.2) as follows:

with . Since we assumed that we know the support of the blurring kernel when we calculated ˜A (while setting all entries in ˜A below a certain threshold to zero), we will also use this information in determining to zero whenever ˜

For the reconstruction we solve the following problem:

The result is shown in Fig. 2d. We observe that the oscillations disappear but there is some loss of contrast compared to the reconstruction using the exact operator (Fig . 2b). These results can be explained as follows. If the feasible set in (1.5) is defined using a noisy operator, it does not contain, in general, the exact solution and therefore no stable approximation of the exact solution can be achieved using elements from this feasible set. Explicitly accounting for the uncertainty in the operator makes the feasible set larger and guarantees the inclusion of the exact solution, enabling stable reconstruction. However, this increased size of the feasible set is also responsible for the loss of contrast observed in Figure 2d as compared to the reconstruction using the exact operator (Figure 2d). The reason for this loss of contrast is that, given the freedom to choose solutions from a larger feasible set, TV seeks to find one with smallest possible contrast.

Deblurring in 2D Let us now turn to two-dimensional images. All images used in this section are greyscale images (with values between 0 and 255) of size 128 128 pixels. Consider first a piecewise constant image shown in Figure 3a. This image is convolved with a Gaussian blur kernel with standard deviation 1. Neumann boundary conditions are used. Then uniform noise with support [= 10 is added to it, which corresponds to 8% noise. The blurred and noisy image is shown in Figure 3b. Reconstruction using the exact operator yields nearly perfect results (Figure 3c). Let us assume, as in the previous section, that only a noisy version ˜A of the forward operator

is available, which we obtain by adding 5% uniform noise to the exact forward operator A:

with . Then we set all entries in ˜A less than d to zero and normalise the rows of ˜A.

Let us reconstruct the image using the noisy operator and isotropic total variation. The result is shown in Figure 3d. We observe numerous artifacts, however, contrast is mainly preserved and corners of the squares are sharp.

Let us now reconstruct the image using interval bounds for the forward operator, which we obtain the same way as in the one-dimensional example. The result is shown in Figure 3e.

Figure 3: Reconstruction of a piecewise-constant image using total variation. Using a noisy operator produces clearly visible artifacts (d). Reconstruction using interval bounds for the operator removes the artifacts at the price of a slight loss of contrast (e, f). Isotropic TV favoures squares with rounded corners (e). Anisotropic TV leaves the corners sharp (f).

The artifacts are removed, but the contrast is slightly reduced, which can be explained the same way as in the one-dimensional example. We also observe that the corners of the squares became rounded, which is the behaviour that we indeed would expect from isotropic TV, which discourages sharp corners.

This illustrates an important feature of the interval-based approach. The feasible set (1.5) obtained using interval bounds for the operator is larger and gives more freedom to the regulariser than a feasible set obtained using a fixed operator (whether exact or noisy). Therefore, features specific to the regulariser are more apparent in the reconstructions in the former case, which explains the rounding of the corners in Figure 3e. To support this idea, let us reconstruct the same image using anisotropic total variation, which does not penalise sharp corners [10]. The result is shown in Figure 3f. Sharp corners are recovered and even the contrast looks a bit better.

Let us now consider another example shown in Figure 4a. This picture is more challenging because it contains much smaller objects, some of which have little contrast from the background, but it is still pieceiwse-constant, making total variation a suitable regulariser.

Already in the reconstruction using the exact operator (Figure 4c) there are some artifacts and the top thin line almost disappeared. The reconstruction quality is even worse when we use the noisy operator (Figure 4d). Artifacts are clearly visible, especially in the bottom line.

Reconstruction using interval bounds for the operator removes the artifacts (Figure 4e), but

Figure 4: Reconstruction of a piecewise-constant image with thin structures using total variation. Even reconstruction using the exact operator produces some artifacts (c). These artifacts are much more clearly present in the reconstruction using a noisy operator (d). Reconstruction using interval bounds for the operator removes the artifacts at the price of some loss of contrast (e, f). The top thin line disappears. Anisotropic TV (f) yields slightly better reconstruction than isotropic TV (e).

causes some lost of contrast, as we already observed in previous examples. The top line, which already had little contrast from the background in the original image, completely disappears in the reconstruction. Using anisotropic total variation further removes the artifacts, but does not recover the top line (Figure 4f).

This loss of contrast and especially the disappearance of some details are, of course, unwanted features of a reconstruction algorithm. However, it is important to note that the solutions shown in Figures 4e and 4f indeed could produce the observed data in Figure 4b (as opposed to the reconstruction in Figure 4d). Therefore, we have no reasons to reject the images in Figures 4e and 4f as plausible reconstructions. By Theorem 3.1 we have that there is a blurring operator A between the interval bounds and a blurry image Au = f, where u is any element of the feasible set (1.5), e.g., the images shown in Figure 4e or 4f. The only way to reject these images would be by narrowing the feasible set by adding more information about the unknown solution or the unknown forward operator. We attempted the latter in Section 4, but came to the conclusion that adding a linear constraint on the forward operator breaks the convexity of the feasible set.

Figure 5: Reconstruction of a natural image using total variation. Reconstruction using a noisy operator produces numerous artifacts (d). In the reconstructions using interval bounds for the operator the artifacts disappear (e, f), but fine details are also removed and the images look cartoon-like.

Deblurring of natural images Let us finally assess the performance of the proposed method on real images. Consider a well-known example, the cameraman (Figure 5a). Let us add blur and noise to it the same way as in the previous examples (Figure 5b).

Reconstruction using the exact operator (Figure 5c) is of reasonable quality, with some loss of fine details, especially in the face of the cameraman and in the camera. Reconstruction using the noisy operator contains numerous artifacts (Figure 5d). In the reconstructions using interval bounds for the operator the artifacts disappear, but also most of the fine details of the image are removed and the image looks very cartoon-like, both with isotropic and anisotropic TV (Figures 5e and 5f).

The reason for the observed behaviour is additional freedom that the larger feasible set gives to the regulariser. Since natural images are rarely piecewise-constant, the piecewise-constant reconstructions that TV promotes don’t look natural. The idea of the approach using interval bounds for the operator is to use less information coming from the forward operator, which is not perfectly known, and instead rely more on the regulariser, which is supposed to reasonably encode prior information about the image. If this is not the case, the approach results in reconstructions that look quite different from the original. Therefore, in order to use the proposed approach on natural images, one needs to try more appropriate regularisers, e.g., TGV [5] or [7]. This step is, however, beyond the scope of the present paper.

In this paper we analysed an approach to image reconstruction problems with uncertainty in the forward operator based on partially ordered spaces. The method is essentially a variant of the residual method with a feasible set based on order intervals. Our main theoretical contribution is the study of this feasible set. It turned out that the feasible set admits two equivalent descriptions, one of which could be modified to include additional a priori constraints on the forward operator. This is especially relevant in deblurring applications, since the rows of a blurring matrix must always sum up to one, which can be expressed as a linear constraint on the matrix. Unfortunately, we came to the conclusion that adding a linear constraint on the forward operator breaks the convexity of the feasible set.

Despite this negative result, we demonstrated in our numerical experiments that the approach is useful in deblurring with an imperfectly known blurring operator, especially in the case when the regulariser well captures qualitative information about the image. Our experiments revealed an important feature of the approach. In the absence of perfect knowledge of the forward operator, the only source of information that can compensate for this lack of knowledge is the regulariser. Therefore, features specific for the regulariser are more apparent in the interval-based method then in standard methods.

[1] Y. A. Abramovich and C. D. Aliprantis. An Invitation to Operator Theory. Vol. 50. Graduate Studies in Mathematics. American Mathematical Society, 2002.

[2] M. Ackermann et al. “Determination of the Point-spread Function for the Fermi Large Area Telescope from On-orbit Data and Limits on Pair Halos of Active Galactic Nuclei”. In: The Astrophysical Journal 765.1 (2013), p. 54.

[3] M Bertero, P Boccacci, G Desider´a, and G Vicidomini. “Image deblurring with Poisson data: from cells to galaxies”. In: Inverse Problems 25.12 (2009), p. 123006. doi: 10.1088/ 0266-5611/25/12/123006.

[4] D. Bertsimas, D. B. Brown, and C. Caramanis. “Theory and Applications of Robust Optimization”. In: SIAM Review 53.3 (2011), pp. 464–501. doi: 10.1137/080734510.

[5] K. Bredies, K. Kunisch, and T. Pock. “Total Generalized Variation”. In: SIAM Journal on Imaging Sciences 3 (2011), pp. 492–526. doi: 10.1137/090769521.

[6] K. Bredies, K. Kunisch, and T. Valkonen. “Properties of : The one-dimensional case”. In: Journal of Mathematical Analysis and Applications 398 (2013), pp. 438–454. doi: 10.1016/j.jmaa.2012.08.053.

[7] M. Burger, K. Papafitsoros, E. Papoutsellis, and C.-B. Sch¨onlieb. “Infimal Convolution Regularisation Functionals of Spaces”. In: Journal of Mathematical Imaging and Vision 55.3 (2016), pp. 343–369. doi: 10.1007/s10851-015-0624-6.

[8] T. F. Chan and C.-K. Wong. “Total variation blind deconvolution”. In: IEEE Transactions on Image Processing 7.3 (1998), pp. 370–375. doi: 10.1109/83.661187.

[9] T. F. Chan and S. Esedoglu. “Aspects of Total Variation Regularized L1 Function Approximation”. In: SIAM Journal on Applied Mathematics 65.5 (2005), pp. 1817–1837. doi: 10.1137/040604297.

[10] S. Esedoglu and S. J. Osher. “Decomposition of images by the anisotropic Rudin-Osher-Fatemi model”. In: Communications on Pure and Applied Mathematics 57.12 (2004), pp. 1609–1626. doi: 10.1002/cpa.20045.

[11] M. Grant and S. Boyd. “Graph implementations for nonsmooth convex programs”. In: Recent Advances in Learning and Control. Ed. by V. Blondel, S. Boyd, and H. Kimura. Lecture Notes in Control and Information Sciences. http://stanford.edu/~boyd/ graph_dcp.html. Springer-Verlag Limited, 2008, pp. 95–110.

[12] M. Grant and S. Boyd. CVX: Matlab Software for Disciplined Convex Programming, version 2.1. http://cvxr.com/cvx. Mar. 2014.

[13] M. Grasmair, M. Haltmeier, and O. Scherzer. “The residual method for regularizing ill-posed problems”. In: Applied Mathematics and Computation 218.6 (2011), pp. 2693 –2710. doi: j.amc.2011.08.009.

[14] V. K. Ivanov, V. V. Vasin, and V. P. Tanana. Theory of Linear Ill-Posed Problems and its Applications. Berlin, Boston: De Gruyter, 2013.

[15] Y. Korolev. “Making use of a partial order in solving inverse problems: II.” In: Inverse Problems 30.8 (2014), p. 085003.

[16] Y. Korolev and A. Yagola. “Making use of a partial order in solving inverse problems”. In: Inverse Problems 29.9 (2013), p. 095012.

[17] D. Perrone and P. Favaro. “A Clearer Picture of Total Variation Blind Deconvolution”. In: IEEE Transactions on Pattern Analysis and Machine Intelligence 38.6 (2016), pp. 1041– 1055. doi: 10.1109/TPAMI.2015.2477819.

[18] R. T. Rockafellar. Convex Analysis. New Jersey: Princeton University Press, 1970.

[19] L. I. Rudin, S. Osher, and E. Fatemi. “Nonlinear total variation based noise removal algorithms”. In: Physica D: Nonlinear Phenomena 60.1 (1992), pp. 259 –268. doi: http: //dx.doi.org/10.1016/0167-2789(92)90242-F.

[20] P. Sarder and A. Nehorai. “Deconvolution methods for 3-D fluorescence microscopy images”. In: IEEE Signal Processing Magazine 23.3 (2006), pp. 32–45. doi: 10.1109/MSP. 2006.1628876.

[21] H. Schaefer. Banach Lattices and Positive Operators. Berlin: Springer, 1974.

[22] P. J. Shaw and D. J. Rawlins. “The point-spread function of a confocal microscope: its measurement and use in deconvolution of 3-D data”. In: Journal of Microscopy 163.2 (1991), pp. 151–165. doi: 10.1111/j.1365-2818.1991.tb03168.x.

[23] J. L. Starck, E. Pantin, and F. Murtagh. “Deconvolution in Astronomy: A Review”. In: Publications of the Astronomical Society of the Pacific 114.800 (2002), p. 1051.

[24] A. N. Tikhonov, A. V. Goncharsky, V. V. Stepanov, and A. G. Yagola. Numerical Methods for the Solution of Ill-Posed Problems. Dordrecht: Kluwer, 1995.

[25] B. Zhang, J. Zerubia, and J.-C. Olivo-Marin. “Gaussian approximations of fluorescence microscope point-spread function models”. In: Appl. Opt. 46.10 (2007), pp. 1819–1829. doi: 10.1364/AO.46.001819.