Sage Journals: Discover world-class research

Abstract

The segmentation of blurred images is of great importance. There have been several recent pieces of work to tackle this problem and to link the areas of image segmentation and image deconvolution in the case where the blur function κ is known or of known type, such as Gaussian, but not in the case where the blur function is not known due to a lack of robust blind deconvolution methods. Here we propose two variational models for simultaneous reconstruction and segmentation of blurred images with spatially invariant blur, without assuming a known blur or a known blur type. Based on our recent work in blind deconvolution, we present two solution methods for the segmentation of blurred images based on implicitly constrained image reconstruction and convex segmentation. The first method is aimed at obtaining a good quality segmentation while the other is aimed at improving the speed while retaining the quality. Our results demonstrate that, while existing models are capable of segmenting images corrupted by small amounts of blur, they begin to struggle when faced with heavy blur degradation or noise, due to the limitation of edge detectors or a lack of strict constraints. We demonstrate that our new algorithms are effective for segmenting blurred images without prior knowledge of the blur function, in the presence of noise and offer improved results for images corrupted by strong blur.

Keywords

Constrained image reconstruction deconvolution denoising convex segmentation

Introduction

Image segmentation is an important technique in image processing which aims to capture either all of the objects of an image^1–5 or only some of them (selective segmentation^6–9). Variational models that partition images based on intensity often employ edge detection techniques to aid the segmentation and some can handle fuzzy boundaries.^10,11 Many of these models can deal with the presence of noise, but blur proves to be more problematic and most variational models struggle to capture all of the required objects, particularly in cases where there is a reliance on the edge detector.

Work in the segmentation of blurred images is at an early stage but there exist models such as those presented in the literature^12–15 which aim to segment blurred images based on Mumford–Shah^4,16–18 or Chan–Vese^4,19,20 segmentation and total variation image restoration.^21,22 In Bar et al.,¹² the authors attempt to segment blurred images by forming a joint functional incorporating segmentation and image reconstruction, assuming the blur type is known (in the so-called semi-blind method). A framework of alternate minimisation is adopted such that the image is simultaneously restored and segmented. Similar techniques can also be found in other work such as Jung et al.¹⁴ As an alternative approach, some models such as those given in the literature^13,15 treat the problems of image reconstruction and segmentation separately in a two-stage model by first restoring the image and then segmenting the restored data.

The main contribution of this work is the proposal of two new models which incorporate blind deconvolution (with implicitly constrained image reconstruction) and convex segmentation. Here, the former task is particularly important in astronomical imaging and medical imaging and offers advantages over hard constraints such as scaling or truncation, while the latter provides a global minimiser to the segmentation problem with no reliance on the initial guess of objects. In a similar manner to Bar et al.,¹² we form a joint functional and proceed to segment the image using alternate minimisation, although it is feasible to use a two-stage approach derived from the joint model (as tested and compared later). We also present an accelerated solution method for acceleration and convergence while sacrificing only a small amount of quality. Our tests will show that related models which do not impose constraints for the restoration (and the restored blur kernel for the blind case) do not perform as well when the underlying blur is heavy.

The rest of this paper is organised as follows. In the ‘Existing methods’ section, we review efforts to tackle the problem of segmenting images which have been corrupted by blur as well as introducing relevant segmentation and image reconstruction models. In the ‘Two-stage restoration and segmentation of images with unknown blur’ section, we introduce two new two-stage models for the cases of segmenting blurred images in the presence of Gaussian noise and Poisson noise incorporating implicitly constrained deblurring and convex segmentation. In the ‘Segmentation of images corrupted by unknown blur’ section, we introduce our main model which is a new joint model for the segmentation of blurred images. In the ‘A relaxed model for the segmentation of images with unknown blur’ section, we introduce our relaxed model for the segmentation of blurred images using alternate direction methods for joint convergence. In the ‘Experimental results’ section, we present experimental results. In the final section, we present the conclusions of this work.

Before proceeding, we make the assumption that the image blur is spatially invariant. As such, we may model the blurred image z as the convolution of the true image u with a kernel function k (also referred to as the blur function or point spread function (PSF)) with the possibility of some additive noise η

z (x, y) = [k * u] (x, y) + η (x, y)

(1)

where the operation of convolution is denoted by * and given by

[k * u] (x, y) = \int - \infty \infty \int - \infty \infty k (x - x', y - y') u (x', y') d x' d y'

Existing methods

In recent years, several approaches have been developed to tackle the problem of accurately segmenting images which have been corrupted by blur. Such approaches most commonly involve image reconstruction to restore the ‘true image’ and segmentation. Models can be classed as two stage in which reconstruction of the corrupted image is carried out, followed by segmentation of the restored image.^13,15 In contrast, there also exist joint models which attempt to deal with the tasks of reconstruction and segmentation simultaneously by minimising a joint functional.^12,14 In this section, we review some existing methods for two stage and joint segmentation of blurred images.

Segmentation of blurred images

Bar et al. showed in their 2004 paper¹² that the two problems of segmentation and image restoration could be coupled together and hence solved at the same time. Both the case of non-blind deconvolution where the kernel is known and the case of semi-blind deconvolution assuming that the blur function is of Gaussian type, leaving the argument σ of the Gaussian function describing the blur to be found, are considered. The problem was solved by minimising the joint functional

f BSK (u, k σ, v) = 12 \int Ω (k σ * u - z) 2 d Ω + β \int Ω v 2 | \nabla u | 2 d Ω + α \int Ω ɛ | \nabla v | 2 + (v - 1) 24 ɛ d Ω + γ \int Ω | \nabla k σ | 2 d Ω

(2)

which is dependent on the image u, the edge integration map v and the kernel function parameter σ. Minimising with respect to the arguments, the authors simultaneously recover and segment the image. A special case exists in the case of known blur where minimisation with respect to the kernel width is not necessary, nor is the final term of the functional.

Jung et al.¹⁴ presented models for joint multi-phase segmentation, deblurring and denoising of images by considering the region-based active contours without edges model and gave the joint formulation as the minimisation of the energy functional

f JCSV (c 1, c 2, φ) = \int Ω | z - k * (c 1 H (φ) + c 2 (1 - H (φ))) | 2 d Ω + μ \int Ω | \nabla H (φ) | d Ω

for two-phase segmentation. A similar model is also presented to allow for the segmentation into m distinct levels. This is then solved by alternate minimisation of the arguments, solving the H¹ gradient descent

φ t (x, t) = - \nabla H 1 f JCSV (c 1, c 2, φ) .

Paul et al.²³ further generalised this idea to image restoration.

Reddy et al.¹⁵ approach this problem in the blind case, where the blur function is an unknown Gaussian using the Chan–Vese model

f RCR = μ Length (Γ) + ν Area (inside (Γ)) + λ 1 \int inside (Γ) | u (x, y) - c 1 | 2 d Ω + λ 2 \int outside (Γ) | u (x, y) - c 2 | 2 d Ω

where Γ is the contour of the segmentation,

u (x, y) = | k (x, y) * z (x, y) | 2

is the square of the convolution of the received data z and a Gaussian kernel, c₁ and c₂ are constants which approximate the image intensities inside and outside of the contour Γ and

μ, ν, λ 1, λ 2

are non-negative parameters. The authors adopt a two-stage approach, first applying a deblurring algorithm, followed by segmentation of the reconstructed image.

Chan et al.¹³ recently presented a two-stage convex method for segmenting blurred images which have been corrupted by either Poisson or multiplicative Gamma noise. Their technique is to extract a smooth image u from the received image z and then to threshold u to reveal segmentation features. The functional to be minimised, given the blurring operator $A$ , is

f CYZ (u) = \int Ω | \nabla u | d Ω + μ 2 \int Ω | \nabla u | 2 d Ω + λ \int Ω A u - f \log A u d Ω

which has a unique solution and can be solved by split-Bregman²⁴ or Chambolle–Pock.²⁵

In this paper, we shall consider both joint and two-stage approaches using convex segmentation and implicitly constrained deblurring in an attempt to improve the quality of the results. Below, we first review the relevant segmentation and image reconstruction techniques.

Two-stage restoration and segmentation of images with known blur

We build a two-stage model for segmenting blurred images in the non-blind case, assuming that we can model precisely the blur degradation, by first deblurring the image and subsequently finding the segmentation of the result. Assuming that we know the function causing the degradation of the image and assuming the presence of additive white Gaussian noise, we may attempt to recover the hidden sharp image by solving an ROF-type minimisation problem²¹ of a functional consisting of a deconvolution fitting term and a regularisation term such as that given by the total variation semi-norm. That is, we solve the minimisation problem

min u (x) {\int Ω ([κ * u] (x) - z (x)) 2 d x + α \int Ω | \nabla u (x) | d x}

(3)

where

κ (x)

is the PSF which describes the blur degradation of the sharp image and is assumed to be known,

z (x)

is the known blurred and noisy image and

u (x)

is the unknown sharp image which is to be restored; * represents the operation of convolution and α is a regularisation parameter which measures the trade-off between data fitting and regularisation. The function

u (x)

which minimises the functional of (3) is the restored image in which we require features to be distinguished using techniques such as Chan–Vese segmentation.⁴ This is formulated as a variational minimisation problem given by

min φ (x), c 1, c 2 {λ 1 \int Ω (c 1 - u (x)) 2 H (φ) d x + λ 2 \int Ω (c 2 - u (x)) 2 (1 - H (φ)) d x + \int Ω | \nabla H (φ (x)) | d x}

(4)

where

c 1, c 2 \in ℝ

are, respectively, the average intensity values inside and outside of the contour defined by the Heaviside H of the indicator function

φ, u (x)

is the image to be segmented which is assumed to be known and

λ 1, λ 2

adjust the weights of the data fitting and contour constraint. The segmentation of blurred images is thus achieved by first solving the problem (3) to obtain the restored image

u (x)

and then segmenting this image by solving the problem (4). This technique may achieve a good result; however, the blur function must be already known and the segmentation result of the non-convex problem (4) is heavily dependent on the restored image which solves the deconvolution problem.

In the following section, we extend this two-stage idea to the blind case, where we do not know the blur function. We attempt to improve the result by building in transformations which allow for the image intensities to be constrained implicitly and by adopting a segmentation technique which is capable of finding the global minimum of the problem.

Segmentation of images corrupted by unknown blur

In this section, we consider the case of segmenting images with unknown blur which allows for more generality in cases where the blur cannot be estimated or found. We build refined blind two-stage models and an improved joint model for solving such problems.

Two-stage restoration and segmentation of images with unknown blur

Building on the previous section, we construct a two-stage model for segmenting blurred images in the blind case. That is, we first aim to reconstruct the sharp image from the corrupted received data without knowledge of the blur function and proceed to segment the result. This is a popular formulation and has been used in work such as Chan et al.¹³ and Reddy et al.¹⁵

Proceeding in the case of blind deconvolution, we may attempt to restore the image and blur function simultaneously. Following work such as Chan and Wong²⁶ and You and Kaveh²⁷ we deblur the image by solving the regularised joint minimisation problem

\arg min {f CW (u, κ) = \int Ω ([κ * u] (x) - z (x)) 2 d Ω + R 1 (u (x)) + R 2 (κ (x))}

(5)

where R₁ and R₂ are regularisation functions which enforce smoothness constraints on the image and blur function, respectively. A solution is sought subject to the constraints

κ (x) \geq 0, u (x) \geq 0, \int Ω κ (x) d Ω = 1, κ (x) = κ (- x)

(6)

which are employed in an attempt to find a unique solution to this non-jointly-convex problem and are imposed explicitly at each outer iteration of an alternate minimisation scheme. Chen et al.²⁸ proposed an improved way of imposing such constraints by enforcing two of them implicitly in the functional, resulting in the minimisation of the functional

f CHWZ (u, κ, ψ, ω) = \int Ω (κ * u - z) 2 d Ω + R 1 (τ a (ψ)) + R 2 (τ b (ω)) + A a (u, ψ; ϕ 1, γ 1) + A b (κ, ω; ϕ 2, γ 2)

(7)

where the functions R₁ and R₂ denote total variation regularisation for the image and blur function, respectively

R 1 (τ a (ψ)) = \int Ω | \nabla τ a (ψ) | d Ω, R 2 (τ b (ω)) = \int Ω | \nabla τ b (ω) | d Ω

and

A

denotes the ADMM term which is used to penalise the distance between the image (respectively blur) function and transformed ψ (respectively ω) functions and which is given by

A a (u, ψ; ϕ 1, γ 1) = γ 1 \int Ω (u - τ a (ψ)) 2 d Ω + < ϕ 1, u - τ a (ψ) >

(8)

Alternate minimisation of the functional is achieved by solving the resulting Euler Lagrange equations

E CHWZ 1 (u, κ, ψ, ϕ 1) = 0, E CHWZ 2 (u, ψ) = 0, E CHWZ 1 (κ, u, ω, ϕ 2) = 0, E CHWZ 2 (κ, ω) = 0

where the functions

E CHWZ 1

and

E CHWZ 2

are given by

E CHWZ 2 (u, ψ) = α 1 \partial τ a \partial ψ \nabla \cdot (\nabla τ a (ψ) | \nabla τ a (ψ) |) + γ 1 (τ a (ψ) - u) \partial τ a \partial ψ - ϕ 1 \partial τ a \partial ψ

(10)

where

k † (x, y) = k (- x, - y)

. We can solve equation (9) efficiently using Fourier transforms and equation (10) by gradient descent techniques. We solve for the blur function k and ω similarly. We recover the image by alternately minimising equation (7), solving equations (9), (10) and the associated Euler Lagrange equations for the blur function, followed by updating the dual function

ϕ 1

until an acceptable tolerance is reached.

Once such u is computed, we then proceed to segment this restored image u by solving equation (11), which allows for the global minimum of the segmentation problem to be found,² replacing the data function z with the restored image u, we have

min 0 \leq ν \leq 1 {μ \int Ω | \nabla ν | d Ω + λ 1 \int Ω | u - c 1 | 2 ν d Ω + λ 2 \int Ω | u - c 2 | 2 (1 - ν) d Ω}

(11)

To enforce the constraint $ν \in [0, 1]$ , we introduce a regularised version of the penalty function, $ς (ν)$ , given by

ς ɛ (ν) = b ɛ (ν) H ɛ (b ɛ (ν)), b ɛ (x) = \sqrt{(2 x - 1) 2 + ɛ} - 1

(12)

where

H ɛ

is the Heaviside approximation defined above. This function and the effect of the ɛ parameter can be seen in Figure 1. When computing the global minimiser,

ν *

, thresholding the function at any value

p \in (0, 1)

gives the contour of the object, Γ _p .

Incorporating this idea into our two-stage approach, we obtain the segmentation of the restored image by solving the problem (11) and enforce the constraint $ν \in [0, 1]$ using the penalty function $ς ɛ (ν)$ given above in equation (12). We then obtain the contour of the object by thresholding the global minimiser $ν *$ at a value $p \in (0, 1)$ . We present the joint segmentation model as

min ν, c 1, c 2 {μ \int Ω | \nabla ν | d Ω + λ 1 \int Ω | u - c 1 | 2 ν d Ω + λ 2 \int Ω | u - c 2 | 2 (1 - ν) d Ω + σ \int Ω b ɛ H ɛ (b ɛ) d Ω}

(13)

where

b ɛ = b ɛ (ν)

. We derive the Euler Lagrange equation by minimising equation (13) with respect to ν. This is given by

E cs (ν) = 0

where

E cs (ν) = μ \nabla \cdot (\nabla ν | \nabla ν |) + λ 1 (u - c 1) 2 - λ 2 (u - c 2) 2 + 2 σ (2 ν - 1) b ɛ + 1 (H ɛ (b ɛ) δ ɛ (b ɛ) b ɛ)

(14)

where

δ ε (ν)

is a smooth approximation of the delta function depending on the parameter ɛ. We solve this equation using additive operator splitting (AOS)^29,30 as follows. Letting

f = - λ ((u - c 1) 2 - (u - c 2) 2) - - 2 σ (2 ν - 1) b (ν) + 1 (H ɛ (b ɛ (ν)) δ ɛ (b ɛ (ν)) b ɛ (ν))

and denoting

W = | \nabla ν | - 1

, we can write

ν t = - E 1 (ν)

in the form

\partial ν \partial t = μ (\partial x (W \partial x ν) + \partial y (W \partial y ν)) + f

After discretisation, we rewrite in the matrix-vector form ( $\overset{\land}{ν} n = ν n + τ f n$ )

ν n + 1 = 12 \sum ℓ = 12 (I - 2 τ μ A ℓ (ν n)) - 1 \overset{\land}{ν} n

(15)

Here, $A ℓ$ is the diffusion quantity in the ℓ direction ( $ℓ = 1, 2$ for x and y directions, respectively). Keeping the other arguments fixed and minimising with respect to c₁ and c₂, we have, respectively, the equations given below

c 1 (u, ν) = \int Ω u ν d Ω \int Ω ν d Ω, c 2 (u, ν) = \int Ω u (1 - ν) d Ω \int Ω (1 - ν) d Ω

(16)

In order to solve this model, we make an initial estimate of the image, which is typically the received data since it is the closest approximation we have. Using this information, we obtain the approximation u of the true image and using this we proceed with alternate minimisation of equation (13) until we reach an acceptable tolerance. Our overall algorithm is presented in Algorithm 1.

Algorithm 1

Segmentation of blurred images: $ν (ℓ) \to A 1 T (ν (0), z, maxit)$

1: $u (0) \to z, ψ (0) \to ξ a (u (0))$

2: $ϕ 1 (0) \to 1$

3: for $ℓ \to 1 : maxit$ do

4: Update $u (ℓ)$ by solving (9)

5: Update $ψ (ℓ)$ by solving (10)

6: Update $ϕ 1 (ℓ) \to ϕ 1 (ℓ - 1) + γ 1 (u (ℓ) - τ a (ψ (ℓ)))$

7: Update $κ (ℓ)$

8: Update $ω (ℓ)$

9: Update $ϕ 2 (ℓ) \to ϕ 2 (ℓ - 1) + γ 2 (κ (ℓ) - τ b (ω (ℓ)))$

10: end for

11: $u \to τ a (ψ ℓ)$

12: for $ℓ \to 1 : maxit$ do

13: Calculate $c 1 (ℓ) \to c 1 (u, ν (ℓ - 1)), c 2 (ℓ) \to c 2$ $(u, ν (ℓ - 1))$ using (16)

14: Update $ν (ℓ)$ by solving $E cs (ν (ℓ - 1)) = 0$ using (14)

15: end for

Considering now the case of Poisson noise being present in the image, we make an alteration to our deblurring algorithm to take this into account. We thus attempt to restore the true image from the corrupted image by solving the Robust Richardson Lucy problem introduced above for the image, employing the function $Φ (s) = 2 \sqrt{s + β}$ .³¹ Solving this problem, we obtain an approximation of the true image. In this two-stage setting, once we have obtained the approximation of the image, we proceed with the segmentation as described in the Gaussian case above. This is outlined in Algorithm 2 below. It can be noted that while this restoration method provides a restriction on the lower bound of the restored image u, it does not provide an upper limit. We may obtain this by a projection $ℙ$ of the restored data onto the ideal range at each iteration.

Algorithm 2

Segmentation of blurred images: $ν (ℓ) \to A 1 T (ν (0), z, maxit)$

1: $u (0) \to z$

2: for $ℓ \to 1 : maxit$ do

3: Update the image $u (ℓ)$

4: Update the blur function $κ (ℓ)$

5: Update transformed blur function $ω (ℓ)$

6: Update $ϕ 2 (ℓ) \to ϕ 2 (ℓ - 1) + γ 2 (κ (ℓ) - τ b (ω (ℓ)))$

7: end for

8: $u \to u ℓ$

9: for $ℓ \to 1 : maxit$ do

10: Calculate $c 1 (ℓ) \to c 1 (u, ν (ℓ - 1)), c 2 (ℓ) \to c 2$ $(u, ν (ℓ - 1))$ using (16)

11: Update $ν (ℓ)$ by solving $E cs (ν (ℓ - 1)) = 0$ using (14)

12: end for

We will demonstrate in Result set 1 in the ‘Experimental results’ section that deblurring considerations are important for obtaining a good segmentation of a blurred image. We also demonstrate in Result sets 1 and 2 that the advanced techniques which we describe in this section are able to give improvements over similar techniques.

In the following section, we introduce our approach to the problem of segmenting blurred images by forming joint models which aim to reduce computation time by simultaneously deblurring and segmenting the image.

A new joint model for the segmentation of images with unknown blur

We now construct a joint variational model for the segmentation of blurred images which includes terms designed for the segmentation of an image and which takes into account the possible presence of blur and noise corruption. By the minimisation of this single energy functional, an image may be simultaneously restored and segmented, thus providing an accurate segmentation of the blurred image. There are two key formats for forming such a joint functional in the literature. First, we may replace the image in the deblurring problem with the binary segmentation and attempt to restore this while recovering the average intensities. While this may provide good results, Paul et al.²³ have shown that this may not be robust. We therefore opt for the second approach as follows. We form the joint model by replacing the received data term $z (x)$ in the segmentation functional (13) by the restored image function $u (x)$ and add the constraint that this function should satisfy the deconvolution minimisation problem (5) and associated constraints (6). Imposing this constraint by incorporating the terms into the existing functional, we form the new joint minimisation model. We present this joint approach as the minimisation of the functional

f (u, c 1, c 2, φ, κ) = λ 1 \int Ω | u - c 1 | 2 H ɛ (φ) d Ω + λ 2 \int Ω | u - c 2 | 2 (1 - H ɛ (φ)) d Ω + μ \int Ω | \nabla H ɛ (φ) | d Ω + 12 | | k * u - z | | L 2 (Ω) 2 + α 1 \int Ω | \nabla u | d Ω + α 2 \int Ω | \nabla k | d Ω

(17)

subject to the constraints (6). Here, the restored image u provides the intensity and spatial information for the segmentation terms following the coefficients

λ 1, λ 2

. Following the ideas for implicitly constrained deconvolution discussed in the previous section, we introduce the functions

ψ, ω

in order to implicitly apply constraints on the intensity values of the image and blur functions. We incorporate the transform functions

τ a (ψ), τ b (ω)

into the functional with terms penalising the distance of the image and blur functions from these terms, resulting in the problem of minimising the functional

f (u, c 1, c 2, φ, κ, ψ, ω, ϕ 1, ϕ 2) = λ 1 \int Ω | τ a (ψ) - c 1 | 2 H ɛ (φ) d Ω + λ 2 \int Ω | τ a (ψ) - c 2 | 2 (1 - H ɛ (φ)) d Ω + μ \int Ω | \nabla H ɛ (φ) | d Ω + 12 | | k * u - z | | L 2 (Ω) 2 + α 1 \int Ω | \nabla τ a (ψ) | d Ω + A a (u, ψ; ϕ 1, γ 1) + α 2 \int Ω | \nabla τ b (ω) | d Ω + A b (κ, ω; ϕ 2, γ 2)

(18)

subject to unit integral and symmetry constraints on the blur function and where

A a

and

A b

are given by equation (8). Finally, we use the idea of finding the global minimum of the segmentation problem. In order to achieve this, we replace the regularised Heaviside of the indicator function with the general function ν and add the condition that this function must have values in the range [0, 1]. This is enforced by incorporating the function

ς ɛ

(12) into the functional (18). Given these considerations, we form our joint functional of global-minimum segmentation and implicitly constrained deconvolution as

f J 1 (ν, c 1, c 2, u, ψ, k, ω; ϕ 1, ϕ 2) = μ \int Ω | \nabla ν | d Ω : segmentationregularisation + λ 1 \int Ω | τ a (ψ) - c 1 | 2 ν d Ω + λ 2 \int Ω | τ a (ψ) - c 2 | 2 (1 - ν) d Ω : segmentationfitting + α 1 \int Ω | \nabla τ a (ψ) | d Ω + A a (u, ψ; ϕ 1, γ 1) : image regularisation + constraint + α 2 \int Ω | \nabla τ b (ω) | d Ω + A b (κ, ω; ϕ 2, γ 2) : blur regularisation + constraint + σ \int Ω b ɛ (ν) H ɛ (b ɛ (ν)) d Ω ︸ segmentationconstraint + 12 | | k * u - z | | L 2 (Ω) 2 ︸ deconvolutionfitting

(19)

In order to solve this model, we derive the partial differential equations defined by the first-order optimality conditions. Solving the resulting equations, we aim to minimise the joint functional. We will take each of the arguments in turn.

Segmentation indicator function ν

Minimising the functional (19) with respect to ν, fixing the other arguments, we derive the Euler Lagrange equation from $\partial f J 1 / \partial ν$ which is given by

μ \nabla \cdot (\nabla ν | \nabla ν |) + λ 1 (τ a (ψ) - c 1) 2 - λ 2 (τ a (ψ) - c 2) 2 + 2 σ (2 ν - 1) b ɛ + 1 (H ɛ (b ɛ) δ ɛ (b ɛ) b ɛ) = 0

(20)

where

b ɛ = b ɛ (ν)

and

δ ε (ν)

is a smooth approximation to the delta function. We solve this equation for the function ν, fixing the remaining functions and parameters, using AOS.^29,30 This method is known to be faster than gradient descent methods such as time marching. Letting

f = - λ ((u - c 1) 2 - (u - c 2) 2) - - 2 σ (2 ν - 1) b (ν) + 1 (H ɛ (b ɛ (ν)) δ ɛ (b ɛ (ν)) b ɛ (ν))

and denoting

W = | \nabla ν | - 1

, we can write

ν t = - E 1 (ν)

in the form

\partial ν \partial t = μ (\partial x (W \partial x ν) + \partial y (W \partial y ν)) + f

After discretisation, we rewrite in the matrix-vector form ( $\overset{\land}{ν} n = ν n + τ f n$ )

ν n + 1 = 12 \sum ℓ = 12 (I - 2 τ μ A ℓ (ν n)) - 1 \overset{\land}{ν} n

(21)

Region average intensity values c₁, c₂

Keeping the other arguments fixed and minimising with respect to c₁ and c₂, we have, respectively, the equations given by $\partial f J 1 / \partial c 1$ and $\partial f J 1 / \partial c 2$ which can be evaluated directly to give the average intensities inside and outside of the segmentation contour

c 1 (ψ, ν) = \int Ω τ a (ψ) ν d Ω \int Ω ν d Ω, c 2 (ψ, ν) = \int Ω τ a (ψ) (1 - ν) d Ω \int Ω (1 - ν) d Ω

(22)

Image function u

Minimising now with respect to u, we have the equation

k † * (k * u - z) + γ 1 (u - τ a (ψ)) + ϕ 1 = 0

(23)

which contains the key deconvolution component. This can be re-written with the left-hand side as a convolution of u as

[k † * k + δ γ 1] * u = k † * z + γ 1 τ a (ψ) - ϕ 1

It is important to note that after the discretisation of this equation, the term $k † * k$ along with the operation of convolution defines a symmetric positive definite (spd) matrix. Put briefly, if k and u are the discretised k and u, respectively, then we have

\bar{k} \circ k \circ u = Hu

for an spd matrix H where ^ denotes the operation of discrete convolution. We can solve this problem using the preconditioned conjugate gradient (pcg) method with a preconditioner, following the idea of Vogel and Oman,³² given by

P = (\tilde{k} * \tilde{k} + χ I) 12 (χ I + (γ 1 + ϕ 1)) (\tilde{k} * \tilde{k} + χ I) 12

(24)

where

\tilde{k}

is a circulant approximation to k.

Transformed image function ψ

Minimising equation (19) with respect to the function ψ, we obtain

E ψ (ψ) = 2 λ 1 (τ a (ψ) - c 1) ν \partial τ a (ψ) ψ + 2 λ 2 (τ a (ψ) - c 2) (1 - ν) \partial τ a (ψ) ψ + α 1 \partial τ a (ψ) ψ \nabla \cdot (\nabla τ a (ψ) | \nabla τ a (ψ) |) - γ 1 (u - τ a (ψ)) \partial τ a (ψ) ψ - ϕ 1 \partial τ a (ψ) \partial ψ

(25)

which may be solved by gradient descent. Letting

ψ = ψ (x; t)

, we solve the problem

ψ t (x, t) = - E ψ (ψ (x, t)) s . t . ψ t (x, t) | t = 0 = ψ (0) (x)

Discretising this equation by forward differences in terms of time t and rearranging, we have

ψ (x, t + 1) = ψ (x, t) - Δ t E ψ (ψ (x, t))

Beginning with the initial estimate of ψ at t = 0 which is determined by the inverse transform of the received data z in the first instance and the latest approximation in subsequent alternate minimisation iterations, we evolve in time until the stopping criteria are met. That is, until the L²-norm of the residual is below a certain tolerance

‖ ψ (x, t + 1) - ψ (x, t) ‖ L 2 = Δ t ‖ E ψ (ψ (x, t)) ‖ L 2 < ψ tol

Point spread (blur) function k

Minimising now with respect to the blur function k, we have the equation for the blur function

u † * (u * k - z) + γ 2 (k - τ b (ω)) + ϕ 2 = 0

(26)

which may be solved for k in a similar manner to equation (23).

Transformed point spread (blur) function ω

Minimising with respect to ω, we obtain

α 2 \partial τ b (ω) ω \nabla \cdot (\nabla τ b (ω) | \nabla τ b (ω) |) - γ 2 (k - τ b (ω)) \partial τ b (ω) \partial ω - ϕ 2 \partial τ b (ω) \partial ω = 0

(27)

which may be solved using a gradient descent scheme.

Overall algorithm

In order to solve this model, we make an initial estimate of the image, which is typically the received data since it is the closest approximation we have. We also make an estimate of the PSF based on visual observation of the received image. Using this information, we obtain the initial estimates of ψ and ω and then calculate the first estimates of c₁ and c₂. We next update the image, ψ, the PSF, ω, the function ν, $ϕ 1$ and $ϕ 2$ . We proceed to iterate until we reach an acceptable tolerance. Our algorithm is presented in Algorithm 3.

Algorithm 3

Segmentation of blurred images: $ν (ℓ) \to A 1 J (ν (0), k (0), z, maxit)$

1: $u (0) \to z, ψ (0) \to ξ a (u (0)), ω (0) \to ξ b (k (0))$

2: $ϕ 1 (0) \to 1, ϕ 2 (0) \to 1$

3: for $ℓ \to 1 : maxit$ do

4: Calculate $c 1 (ℓ) \to c 1 (ψ (ℓ - 1), ν (ℓ - 1)), c 2 (ℓ) \to c 2$ $(ψ (ℓ - 1), ν (ℓ - 1))$ using (22)

5: Update $u (ℓ)$ by solving (23)

6: Update $ψ (ℓ)$ by solving (25)

7: Update $k (ℓ)$ by solving (26)

8: Update $ω (ℓ)$ by solving (27)

9: Update $ν (ℓ)$ by solving (20)

10: Update $ϕ 1 (ℓ) \to ϕ 1 (ℓ - 1) + γ 1 (u (ℓ) - τ a (ψ (ℓ)))$

11: Update $ϕ 2 (ℓ) \to ϕ 2 (ℓ - 1) + γ 2 (k (ℓ) - τ b (ω (ℓ)))$

12: end for

We demonstrate in Result set 2 of the ‘Experimental results’ section that segmenting a blurred image by the joint approach can offer improved results over the corresponding two-stage method. Furthermore, in Result set 3, we show that our method offers improved results over comparable methods. In the following section, we consider an alternative joint method which aims to improve the speed of obtaining a solution.

A relaxed model for the segmentation of images with unknown blur

In order to speed up the above model, we consider a way to simplify the equation for ψ through relaxation of the functional. We introduce into the segmentation fitting terms a new variable $ϖ (x, y)$ which should be equal to the sharp image at convergence. We achieve this by adding distance measures which penalise the difference between the restored image and the function ϖ, thus driving them close together. We add penalisation parameters so that the influence of this relationship can be tweaked. We present our joint problem functional as

max ϕ 1, ϕ 2, ζ min ν, c 1, c 2, ϖ, u, ψ, k, ω {f J 2 (ν, c 1, c 2, ϖ, u, ψ, k, ω; ϕ 1, ϕ 2, ζ) = μ \int Ω | \nabla ν | d Ω + λ 1 \int Ω | ϖ - c 1 | 2 ν d Ω + λ 2 \int Ω | ϖ - c 2 | 2 (1 - ν) d Ω + σ \int Ω b ɛ (ν) H ɛ (b ɛ (ν)) d Ω + 12 | | k * u - z | | L 2 (Ω) 2 + α 1 \int Ω | \nabla τ a (ψ) | d Ω + γ 12 | | u - τ a (ψ) | | L 2 (Ω) 2 + < ϕ 1, u - τ a (ψ) > + υ 2 | | ϖ - τ a (ψ) | | L 2 (Ω) 2 + < ζ, ϖ - τ a (ψ) > + α 2 \int Ω | \nabla τ b (ω) | d Ω + γ 22 | | k - τ b (ω) | | L 2 (Ω) 2 + < ϕ 2, k - τ b (ω) >}

(28)

Figure 1.

Illustration of the continuous approximation $ς ɛ$ to the piecewise linear function ς. For lower ɛ, the approximation is very close to ς. (a) $ς (x)$ , (b) $ς ɛ (x), ɛ = 10 - 2$ and (c) $ς ɛ (x), ɛ = 10 - 3$ .

In order to solve this model, we derive the partial differential equations defined by the first-order optimality conditions. We will take the Euler–Lagrange equations for each argument in turn. Minimising with respect to ν, we obtain the equation derived from $\partial f J 2 / \partial ν$

E 1 (ν) = μ \nabla \cdot (\nabla ν | \nabla ν |) + λ 1 \int Ω | ϖ - c 1 | 2 d Ω - λ 2 \int Ω | ϖ - c 2 | 2 d Ω + 2 σ (2 ν - 1) b ɛ + 1 (H ɛ (b ɛ) + δ ɛ (b ɛ) b ɛ)

(29)

where

b ɛ = b ɛ (ν)

. We find a solution to

E 1 (ν) = 0

(29) using AOS as described earlier.

Minimising $f J 2$ with respect to u, we obtain the equation

k † * k * u + (γ 1 + ϕ 1) u = k † * z + (γ 1 + ϕ 1) τ a (ψ)

(30)

where

k † (x, y) = k (- x, - y)

. This can be solved quickly using pcg.³³

Minimising $f J 2$ now with respect to ψ, we obtain the equation

E 2 (ψ) = α R (ψ) - τ a' (ψ) (γ 1 + ϕ 1) (u - τ a (ψ)) + τ a' (ψ) (υ (ϖ - τ a (ψ))) - ζ † τ a (ψ))

(31)

where

R (ψ)

is the derivative of the regularisation term

\int Ω | \nabla τ a (ψ) | d Ω

. We can solve

E 2 (ψ) = 0

(31) using semi-implicit time marching,

ψ t = - E 2 (ψ)

by discretising the time step.

Now, minimising $f J 2$ with respect to ϖ, we obtain

E 3 (ϖ) = 2 λ 1 (ϖ - c 1) ν + 2 λ 2 (ϖ - c 2) (1 - ν) + u (ϖ - τ a (ψ)) + ζ † (ϖ - τ a (ψ))

Note that we can solve the sub-problem $E 3 (ϖ) = 0$ directly by giving the solution to

(2 λ 1 ν + 2 λ 2 (1 - ν) + u + ζ †) ϖ = 2 λ 1 c 1 ν + 2 λ 2 c 2 (1 - ν) + (u + ζ †) τ a (ψ)

(32)

To obtain the average intensities, we minimise $f J 2$ with respect to c₁ and c₂, we obtain equations which can be evaluated directly

c 1 (ϖ, ν) = τ a (\int Ω ϖ ν d Ω \int Ω ν d Ω), c 2 (ϖ, ν) = τ a (\int Ω ϖ (1 - ν) d Ω \int Ω 1 - ν d Ω)

(33)

Finally, minimising with respect to k and ω, we obtain equations (26) and (27), respectively.

In order to solve the model (28), we make an initial estimate of the image, which we allow to be given by the received data z since this is the closest approximation to the true data that we have. Alternatively, if we know or can make an estimate of the PSF, we may solve a Tikhonov model²² and attempt to use this as the initial estimate based on visual judgement. We then calculate the initial estimate of ψ as the inverse transform of the initial estimate of the image. Similarly, in the blind case, we make an initial estimate of the PSF based on visual observation and compute its inverse transform function. We next make an initial estimate of the contour, obtaining the initial estimate of ν. Using these and equation (33), we make the initial estimates of c₁ and c₂. We then proceed to solve the model (28), alternately minimising with respect to the arguments. The final segmentation is then given by the contour Γ _p derived from the final function ν. We present this algorithm in Algorithm 4 below.

Algorithm 4

Segmentation of blurred images: $ν (ℓ) \to A 2 J (ν (0), k (0), z, maxit)$

1: $u (0) \to z$ , $ϖ (0) \to u (0)$ , $ψ (0) \to ξ a (u (0))$ , $ω (0) \to ξ b (k (0))$

2: $ϕ 1 (0) \to 1, ϕ 2 (0) \to 1, ζ (0) \to 1$

3: for $ℓ \to 1 : maxit$ do

4: Calculate $c 1 (ℓ) \to c 1 (ϖ (ℓ - 1), ν (ℓ - 1)), c 2 (ℓ) \to c 2$ $(ϖ (ℓ - 1), ν (ℓ - 1))$ using (33)

5: Update $u (ℓ)$ using (30)

6: Update $ψ (ℓ)$ using (31)

7: Update $k (ℓ)$

8: Update $ω (ℓ)$

9: Update $ν (ℓ)$ using (29)

10: Update $ϖ (ℓ)$ using (32) 11: Update $ϕ 1 (ℓ) \to ϕ 1 (ℓ - 1) + γ 1 (u (ℓ) - τ a (ψ (ℓ)))$ 12: Update $ϕ 2 (ℓ) \to ϕ 2 (ℓ - 1) + γ 2 (u (ℓ) - τ b (ω (ℓ)))$ 13: Update $ζ (ℓ) \to ζ (ℓ - 1) + υ (τ a (ψ (ℓ)) - w (ℓ))$ 14: end for

Experimental results

Attempting to segment a blurred image with a segmentation technique (such Chan–Vese⁴) is sufficient to obtain a close result if the degradation is not strong but as the amount of corruption increases, it is not possible to obtain a good result because such models are not designed with blur degradation in mind. Meanwhile, the work of Bar et al.¹² is capable of segmenting blurred images where the corruption is small but begins to struggle to obtain good quality results in the presence of significant blur degradation or noise.

In this section, through experiments, we demonstrate that Algorithms 1 and 2 offer improvements over competing models for blurred images. We also show that Algorithm 3 is capable of obtaining a good quality result with the possibility of slow convergence while Algorithm 4 converges faster to a similar, if slightly lower, quality.

We present results of segmenting the following images (see Figure 2) with the addition of varying levels of Gaussian blur and noise: Im₁: Text (Figure 2(a)), Im₂: Cells medical (Figure 2(b)), Im₃: Box-Triangle (Figure 2(c)), Im₄: QR Code (Figure 2(d)), Im₅: Fingerprint (Figure 2(e)), Im₆: Tree (Figure 2(f)).

Figure 2.

Images used for test examples: (a) Im₁, (b) Im₂, (c) Im₃, (d) Im₄, (e) Im₅ and (f) Im₆.

We denote by $Γ Mod 1$ and $ν Mod 1$ , respectively, the contour obtained and the segmented area obtained by solving Mod₁. The notation is similar for the other models.

Experiments were carried out using Matlab R2013a on a HPE-595uk with an Intel(R) Core(TM) i7-2600 processor and 16GB RAM.

Models

In order to compare our results with competing and other relevant models, we define the following models to be tested in this section:

Mod₁: The Chan–Vese segmentation model (CV).⁴

Mod₂: The two-stage model by standard TV deblurring followed by CV segmentation.

Mod₃: The two-stage model by standard deblurring for Poisson noise followed by CV segmentation.

Mod₄: The Bar et al. model¹² using equation (2) – without constraints on $k σ, u$ .

$New G T$ : Algorithm 1 – a two-stage model given by implicitly constrained blind deblurring (7) for Gaussian noise followed by convex segmentation (13).

$New P T$ : Algorithm 2 – adapted Algorithm 1 for Poisson noise.

New ^J : Algorithm 3 – the joint minimisation model (19) for blind deblurring and convex segmentation, with built-in constraints on k, u.

$New R J$ : Algorithm 4 – the relaxed joint minimisation model (28) from (19).

Measuring error

In order to make a numerical evaluation of our model, we require a ground truth. For the artificial images Im₁ and Im₃, we already know the contour and consider this to be the ground truth segmentation. For the remaining images, we estimate the true contour by assuming that the segmentation of the true (uncorrupted) image is correct (see Figure 3) and we consider methods of measuring the accuracy of the final contour Γ, that is how close it is to the segmentation by Chan–Vese of the true image given by $Γ true$ . To this end, we compare the final contours and indicator functions. We compare results using the following error measures, of which Er₁ and Er₂ tend to zero and Er₃ and Er₄ tend to one as the segmentation of the blurred image tends towards the segmentation of the true image. Throughout, we denote the true, received and restored image u as $u true, u rec$ and $u rest$ , respectively, with similar notation for the remaining functions. To distinguish models, we also denote by $Γ New J$ the contour obtained by solving model New ^J and adopt similar notation for the remaining models and functions.

L2 area-based difference gives the L₂ norm of the difference in segmented images. It measures the closeness of the final indicator functions

Er 1 = | | ν rest - ν true | | 22

where

ν rest

corresponds to the function achieved from solving the model.

Contour difference gives the L₂ norm of the difference between final contours

Er 2 = | | Γ rest - Γ true | | 22

Letting the set of points which are considered to be inside in the contour be given by

S rest = {(x, y) \in Ω | ν rest (x, y) > 10 - 1} S true = {(x, y) \in Ω | ν true (x, y) > 10 - 1}

we define Er₃ (Tanimoto Coefficient)^34–36 and Er₄ (Dice Similarity Coefficient)^34,37 as

Er 3 = N (S rest \cap S true) N (S rest \cup S true) and Er 4 = 2 N (S rest \cap S true) N (S rest) + N (S true)

respectively where

N (S)

denotes the number of elements in S. It is clear that as the restored segmentation tends towards the true segmentations, both error values tend towards one.

Figure 3.

Segmentation of images Im₁–Im₆ using model Mod₁: (a) Im₁ segemented by Mod₁, (b) Im₂ segemented by Mod₁, (c) Im₃ segemented by Mod₁, (d) Im₄ segemented by Mod₁, (e) Im₅ segemented by Mod₁ and (f) Im₆ segemented by Mod₁.

Result sets

We group our experimental results by the following result sets:

RS₁: Result set 1 consists of images corrupted by blur with the assumption that Gaussian noise is present. We illustrate the performance of Mod₁ to segment the image and consider it against the performance of Mod₂ and $New G T$ . We see in Table 1 and Figure 4 that while Mod₁ can give a reasonable result, it is not reliable for segmenting blurred images. We also demonstrate in Table 1 that the result can be improved by using the two-stage model Mod₂ and further enhanced by the constrained model $New G T$ . We also demonstrate in Table 2 that this idea carries over to the case of Poisson noise corruption in the image. It can be seen that model $New P T$ offers an improvement over the two-stage Poisson model Mod₃.

Table 1.

Result set 1. Error values for Im₁–Im₆ corrupted by Gaussian blur and segmented by Mod₁. Values in bold indicate the lowest error achieved for each image.

		Small Gaussian Blur, σ = 9			Large Gaussian Blur, σ = 19
Image	Model	Er ₁	Er ₂	Er ₃	Er ₁	Er ₂	Er ₃
Im₁	Initial	101.16	64.34	0.13	101.16	64.34	0.13
	Mod₁	101.86	86.19	0.37	141.48	82.47	0.22
	Mod₂	15.32	22.54	0.96	44.82	65.81	0.85
	$New G T$	13.09	19.85	0.97	17.24	25.46	0.95
Im₂	Initial	104.31	32.03	0.23	104.31	32.03	0.23
	Mod₁	32.27	36.30	0.89	63.59	37.38	0.67
	Mod₂	18.99	23.73	0.97	14.73	28.28	0.89
	$New G T$	7.27	12.74	0.99	10.50	16.45	0.99
Im₃	Initial	109.65	34.84	0.42	109.65	34.84	0.42
	Mod₁	16.79	23.81	0.99	32.10	36.70	0.95
	Mod₂	1.62	0.00	1.00	1.71	0.00	1.00
	$New G T$	1.57	0.00	1.00	1.56	0.00	1.00
Im₄	Initial	138.19	97.49	0.19	138.19	97.49	0.19
	Mod₁	81.87	100.09	0.74	142.16	91.71	0.45
	Mod₂	21.69	30.98	0.89	48.71	68.60	0.88
	$New G T$	17.94	25.92	0.98	20.29	29.17	0.98
Im₅	Initial	109.07	104.47	0.29	109.07	104.47	0.29
	Mod₁	120.87	109.63	0.49	134.30	109.44	0.42
	Mod₂	44.96	66.06	0.87	46.33	65.98	0.86
	$New G T$	16.11	24.54	0.98	27.18	38.99	0.95
Im₆	Initial	145.97	70.51	0.27	145.97	70.51	0.27
	Mod₁	69.01	75.80	0.85	84.28	75.87	0.79
	Mod₂	46.22	67.20	0.88	35.91	52.68	0.95
	$New G T$	18.01	26.50	0.99	27.74	40.31	0.97

Note: In many cases, the competition is close but $New G T$ obtains the same or improved error values over competing models in all cases.

Table 2.

Result set 1. Error values for Im₁–Im₆ corrupted by Gaussian blur and segmented by Mod₃ and $New P T$ . Values in bold indicate the lowest error achieved for each image.

		Small Gaussian Blur, σ = 9			Large Gaussian Blur, σ = 19
Image	Model	Er ₁	Er ₂	Er ₃	Er ₁	Er ₂	Er ₃
Im₁	Mod₁	101.86	86.19	0.37	141.48	82.47	0.22
	Mod₃	72.97	77.32	0.54	96.44	78.63	0.39
	$New P T$	67.52	75.39	0.57	90.31	74.70	0.43
Im₂	Mod₁	32.27	36.30	0.89	63.59	37.38	0.67
	Mod₃	11.74	18.63	0.98	25.72	33.06	0.92
	$New P T$	11.74	18.47	0.98	25.50	32.94	0.92
Im₃	Mod₁	16.79	23.81	0.99	32.10	36.70	0.95
	Mod₃	12.29	16.94	0.99	12.49	17.49	0.99
	$New P T$	11.97	16.94	0.99	12.28	17.29	0.99
Im₄	Mod₁	81.87	100.10	0.74	142.16	91.72	0.45
	Mod₃	57.22	81.51	0.86	110.57	89.48	0.61
	$New P T$	55.51	79.34	0.86	110.51	89.40	0.61
Im₅	Mod₁	120.87	109.63	0.49	134.30	109.44	0.42
	Mod₃	102.98	103.05	0.57	104.25	103.54	0.56
	$New P T$	102.95	103.05	0.57	104.21	103.51	0.56
Im₆	Mod₁	69.01	75.80	0.85	84.29	75.87	0.79
	Mod₃	50.96	69.64	0.91	57.27	70.41	0.89
	$New P T$	50.96	69.64	0.91	57.27	70.38	0.89

Note: The competition is close for most examples, but overall $New P T$ outperforms Mod₃.

RS₂: In Result set 2, we consider the comparison between the results obtained by solving the two-stage model

New G T

with those obtained as the solutions of the joint models New ^J and

New R J

. We can see in Table 3 that there is little advantage in considering the problems of deblurring and segmentation separately. In most cases, the joint models achieve better results than the two-stage model while the two-stage model obtains an improved result in only one case.

Table 3.

Result set 2. Error values given by Er₁ for Im₁–Im₄ corrupted by Gaussian blur and segmented by $New G T, New J$ and $New R J$ . Values in bold indicate the lowest error achieved for each image.

Image	Initial	$New G T$	New ^J	$New R J$
Im₁	101.1621	13.0888	15.3386	16.5861
Im₂	104.3123	7.2681	5.8182	6.9165
Im₃	109.6528	1.5731	1.0923	1.5249
Im₄	138.1870	17.9387	14.3647	17.6317

Note: For Im₁, $New G T$ outperforms the other models but in the remaining cases New ^J and $New R J$ obtain improved results.

RS₃: Result set 3 consists of images corrupted by small and large amounts of blur as well as noise. We demonstrate in Table 4 and Figures 4 to 6 that model Mod₁ is sometimes sufficient to obtain a fairly close result but misses a considerable amount of detail while Mod₄ can give an improvement over this. We also demonstrate that the new joint models New ^J and

New R J

are capable of segmenting these examples and offer further enhancement over Mod₄. Further to this, we demonstrate in Table 5 and Figure 7 that, when the level of blur is larger, New ^J offers further improvements over Mod₄. Finally, we demonstrate in Table 6 and Figure 8 that, as the level of noise is increased, New ^J continues to perform well.

Table 4.

Result sets 3, 4. Error values and cpu times for images Im₁–Im₆, ‘Circles’ and ‘Knee’ images corrupted by small Gaussian blur. Values in bold indicate the lowest error and least CPU time achieved for each image.

Model	Er ₁	Er ₂	Er ₃	Er ₄	cpu
Im₁ corrupted by small Gaussian blur
Initial	101.16	64.34	0.13	0.23
Mod₁	101.86	86.19	0.37	0.54	317
Mod₄	94.94	85.99	0.40	0.58	1362
New ^J	15.34	23.19	0.96	0.98	559
$New R J$	16.59	24.54	0.95	0.98	443
Im₂ corrupted by small Gaussian blur
Initial	104.31	32.03	0.23	0.37
Mod₁	32.27	36.30	0.89	0.94	527
Mod₄	23.02	31.94	0.94	0.97	861
New ^J	5.82	11.14	0.99	0.99	546
$New R J$	6.92	12.33	0.99	0.99	427
Im₃ corrupted by small Gaussian blur
Initial	109.65	34.84	0.42	0.59
Mod₁	16.79	23.81	0.99	0.99	319
Mod₄	13.64	19.44	0.99	0.99	634
New ^J	1.09	0	1	1	550
$New R J New R J$

Note: For Im₁, $New G T$ outperforms the other models but in the remaining cases New ^J and $New R J$ obtain improved results.

Figure 4.

Result sets 1, 3. Illustration of the performance of the Mod₁ for Im₁ corrupted by Gaussian blur: (a) initial contour, (b) segmentation given by Mod₁, (c, d) segmentation given by New ^J . Mod₁ gives a rough segmentation while the spaces between the letters which are hidden by the blur are successfully segmented using New ^J .

Figure 5.

Result sets 3, 4. Illustration of the performance of the New ^J for Im₂ corrupted by Gaussian blur: (a) received data, (b, c) segmentation using New ^J , (d) the difference between the segmentation using New ^J and using Mod₁. The segmentation is closer to the true edge using New ^J while Mod₁ also captures the blurred edge.

Figure 6.

Result sets 3, 4. Illustration of the performance of New ^J for (top-bottom) Im₄, Im₃, Im₅ and Im₆ corrupted by Gaussian blur. The edges hidden by blur are successfully segmented by New ^J which cannot be segmented by Mod₁.

Figure 7.

Result set 3. Illustration of the performance of the New ^J for (top-bottom) Im₁, Im₄, Im₂ and Im₆ corrupted by strong Gaussian blur. New ^J is capable of segmenting edges in these challenging cases which cannot be segmented by Mod₁.

Figure 8.

Result set 3. Illustration of the performance of the New ^J for (top-bottom) Im₁, Im₃, Im₄ and Im₅ corrupted by Gaussian blur and noise. The edges hidden by blur are successfully segmented by New ^J which cannot be segmented by Mod₁.

New R J

Table 5.

Result set 3. Error values and cpu times for images Im₁–Im₆, ‘Circles’ and ‘Knee’ images corrupted by strong Gaussian blur. Values in bold indicate the lowest error or least CPU time achieved for each image.

Model	Er ₁	Er ₂	Er ₃	Er ₄	cpu
Im₁ corrupted by large Gaussian blur
Initial	101.16	64.34	0.13	0.23
Mod₁	141.48	82.47	0.22	0.36	520
Mod₄	151.08	80.40	0.21	0.35	1479
New ^J	35.68	54.83	0.80	0.89	557
$New R J$	44.93	62.98	0.75	0.86	549
Im₂ corrupted by large Gaussian blur
Initial	104.31	32.03	0.23	0.37
Mod₁	63.59	37.38	0.67	0.80	522
Mod₄	47.17	36.88	0.79	0.88	950
New ^J	6.79	13.86	0.99	0.99	566
$New R J$	10.17	16.31	0.99	0.99	434
Im₃ corrupted by large Gaussian blur
Initial	109.65	34.84	0.42	0.59
Mod₁	32.10	36.70	0.95	0.97	418
Mod₄	28.31	36.76	0.96	0.98	730
New ^J	11.96	16.88	0.99	0.99	552
$New R J$	12.08	16.97	0.99	0.99	232
Im₄ corrupted by large Gaussian blur
Initial	138.19	97.49	0.19	0.32
Mod₁	142.16	91.72	0.45	0.62	525
Mod₄	140.11	91.25	0.49	0.66	1477
New ^J	25.03	38.54	0.96	0.98	589
$New R J$	30.42	43.59	0.95	0.98	444
Model	Er ₁	Er ₂	Er ₃	Er ₄	cpu
Im₅ corrupted by large Gaussian blur
Initial	109.07	104.47	0.29	0.45
Mod₁	134.30	109.44	0.42	0.59	523
Mod₄	126.19	110.54	0.46	0.63	937
New ^J	74.10	91.39	0.71	0.83	564
$New R J$	74.77	92.98	0.69	0.81	334
Im₆ corrupted by large Gaussian blur
Initial	145.97	70.51	0.27	0.43
Mod₁	84.29	75.87	0.79	0.88	520
Mod₄	76.36	76.09	0.82	0.90	940
New ^J	40.19	58.20	0.94	0.97	574
$New R J$	42.54	60.44	0.94	0.97	539
‘Circles’ corrupted by large Gaussian blur
Initial	65.21	28.28	0.69	0.82
Mod₁	114.90	32.76	0.47	0.64	525
Mod₄	124.12	33.08	0.44	0.61	1364
New ^J	26.29	29.82	0.93	0.96	543
$New R J$	24.81	29.75	0.95	0.97	430
‘Knee’ corrupted by large Gaussian blur
Initial	120.54	60.70	0.32	0.49
Mod₁	89.11	63.17	0.70	0.82	522
Mod₄	87.84	63.11	0.72	0.84	1148
New ^J	47.29	56.39	0.89	0.94	548
$New R J$	43.29	55.53	0.90	0.95	539

Note: In all cases, New ^J and $New R J$ achieve improved results and competition is close between New ^J and $New R J$ . For most cases, the cpu time is lower for $New R J$ with the exception of three examples which have slightly lower cpu time for Mod₁ with deteriorated results.

Table 6.

Result set 3. Error values and cpu times for Im₁, Im₃–Im₅ corrupted by Gaussian blur and zero-mean Gaussian noise of variance 0.005. Values in bold indicate the lowest error and least CPU time achieved for each image.

Model	Er ₁	Er ₂	Er ₃	Er ₄	cpu
Im₁ corrupted by Gaussian blur and noise
Initial	101.16	64.34	0.13	0.23
Mod₁	101.84	86.26	0.37	0.54	319
Mod₄	94.93	86.12	0.40	0.58	1899
New ^J	13.86	19.85	0.97	0.98	607
$New R J$	14.90	20.95	0.96	0.98	452
Im₃ corrupted by Gaussian blur and noise
Initial	109.65	34.84	0.42	0.59
Mod₁	16.86	23.94	0.99	0.99	317
Mod₄	14.32	21.54	0.99	0.99	633
New ^J	1.75	0	1	1	512
$New R J$	2.18	0	1	1	242
Model	Er ₁	Er ₂	Er ₃	Er ₄	cpu
Im₄ corrupted by Gaussian blur and noise
Initial	138.19	97.49	0.19	0.32
Mod₁	81.91	99.99	0.74	0.85	530
Mod₄	65.71	84.35	0.82	0.90	958
New ^J	26.84	36.25	0.96	0.98	640
$New R J$	23.90	32.14	0.97	0.99	342
Im₅ corrupted by Gaussian blur and noise
Initial	109.07	104.47	0.29	0.45
Mod₁	120.86	109.64	0.49	0.65	547
Mod₄	113.25	109.98	0.52	0.68	877
New ^J	41.61	62.13	0.87	0.93	600
$New R J$	45.32	63.95	0.86	0.93	575

Note: In all cases, New ^J and $New R J$ achieve improved results. Cpu time is lower for $New R J$ in two cases. In the remaining cases, it is lower for Mod₁ and closely followed by $New R J$ which achieved significantly improved results.

RS₄: Result set 4 demonstrates the ability of $New R J$ and compares the performance of this model with New ^J . We can see in in Tables 4 –6 and Figure 9 that it is generally the case that $New R J$ is faster than New ^J while New ^J obtains better results.

Figure 9.

Result set 4. Images corrupted by Gaussian blur segmented using $New R J$ .

Conclusions

We have proposed a new model for the effective segmentation of blurred images in the blind case where the blur function is unknown (New ^J ) and presented results demonstrating its ability to capture edges which are blurred and difficult to segment closely as well as edges that are hidden by blur. We have also presented an accelerated model ( $New R J$ ) which is also capable of achieving good results with similar examples. We have also demonstrated that these joint models can obtain improved performance over comparable two-stage models. This model can be further extended to the semi-blind case where some information about the blur function may be assumed to be known^38–43 which allows for increases in speed and to work with multi-channel images. This model will be considered for selective segmentation⁹ and vessel segmentation techniques among others.

Footnotes

Acknowledgements

The first and the third authors thank the UK EPSRC for an industrial CASE studentship (voucher 10002441) jointly with St Paul’s Eye Unit, Royal Liverpool University Hospital. The second author thanks the University of Liverpool for a GTA award.

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: The third author acknowledges the partial support from the UK EPSRC grant EP/K036939/1.

References

Bresson

Esedoglu

Vandergheynst

. Fast global minimization of the active contour/snake model. J Math Imaging Vis 2007; 28: 151–167.

Chan

Esedoglu

Nikilova

. Algorithms for finding global minimizers of image segmentation. SIAM J Appl Math 2006; 66: 1632–1648.

Chan

Sandberg

Vese

. Active contours without edges for vector-valued images. J Vis Commun Image R 2000; 11: 130–141.

Chan

Vese

. Active contours without edges. IEEE T Image Process 2001; 10: 266–277.

Kao

Gore

. Minimization of region-scalable fitting energy for image segmentation. IEEE T Image Process 2008; 17: 1940–1949.

Badshah

Chen

. Image selective segmentation under geometric constraints using an active contour approach. Commun Comput Phys 2009; 7: 759–778.

Gout

Le Guyader

Vese

. Segmentation under geometrical conditions using geodesic active contours and interpolation using level set methods. Numer Algorithms 2005; 39: 155–173.

Nguyen

TNA

Cai

Zhang

. Robust interactive image segmentation using convex active contours. IEEE T Image Process 2012; 21: 3734–3743.

Rada

Chen

. A new variational model with dual level set functions for selective segmentation. Commun Comput Phys 2012; 12: 261–283.

10.

Barchiesi

Kang

. A variational model for infinite perimeter segmentations based on lipschitz level set functions: denoising while keeping finely oscillatory boundaries. Multiscale Model Simul 2010; 8: 1715–1741.

11.

Rada

Chen

. On a variational model for selective image segmentation of features with infinite perimeter. J Math Res Appl 2013; 33: 253–272.

12.

Bar L, Sochen N and Kiryati N. Variational pairing of image segmentation and blind restoration. In: Computer vision – ECCV 2004. Berlin: Springer, 2004, pp.166–177.

13.

Chan

Yang

Zeng

. A two-stage image segmentation method for blurry images with Poisson or multiplicative gamma noise. SIAM J Imaging Sci 2014; 7: 98–127.

14.

Jung M, Chung G, Sundaramoorthi G, et al. Sobolev gradients and joint variational image segmentation, denoising, and deblurring. In: IS&T/SPIE electronic imaging. San Jose, CA: International Society for Optics and Photonics, 2009, pp.72460I–72460I.

15.

Reddy R, Chandra M and Rao R. Segmentation of blurred images using improved Chan-Vese snake model. In: Int Conf Signal Processing Communication Computing and Networking Technologies (ICSCCN), IEEE, New York, NY, USA, 2011, pp.502–505.

16.

Ambrosio

Tortorelli

. On the approximation of free discontinuity problems. B Unione Mat Ital 1992; 7: 105–123.

17.

Kim J, Tsai A, Cetin M, et al. A curve evolution-based variational approach to simultaneous image restoration and segmentation. In: IEEE Int Conf Image Processing (ICIP), IEEE, New York, NY, USA, 2002, volume 1, pp.I–109–112.

18.

Mumford

Shah

. Optimal approximation by piecewise smooth functions and associated variational problems. Commun Pure Appl Math 1989; 42: 577–685.

19.

Chung G. A unifying framework for piecewise-constant image segmentation and deblurring. PhD thesis, University of California, Los Angeles, 2007.

20.

Chung G and Vese LA. Energy minimization based segmentation and denoising using a multilayer level set approach. In: Energy minimization methods in computer vision and pattern recognition (EMMCVPR). Berlin: Springer, 2005, pp.439–455.

21.

Rudin

Osher

Fatemi

. Nonlinear total variation based noise removal algorithms. Phys D 1992; 60: 259–268.

22.

Vogel

. Computational methods for inverse problems, Philadelphia, PA: SIAM, 2002.

23.

Paul

Cardinale

Sbalzarini

. Coupling image restoration and segmentation: a generalized linear model/Bregman perspective. Int J Comput Vis 2013; 104: 69–93.

24.

Goldstein

Osher

. The split Bregman method for l1-regularized problems. SIAM J Imaging Sci 2009; 2: 323–343.

25.

Chambolle

Pock

. A first-order primal-dual algorithm for convex problems with applications to imaging. J Math Imaging Vis 2011; 40: 120–145.

26.

Chan

Wong

. Total variation blind deconvolution. IEEE T Image Process 1998; 7: 370–375.

27.

You

Kaveh

. A regularization approach to joint blur identification and image restoration. IEEE T Image Process 1996; 5: 416–428.

28.

Chen

Harding

Williams

. A new study of blind deconvolution with implicit incorporation of non-negativity constraints. Int J Comput Math 2015, pp. 860263.

29.

Neittaanmäki

Tai

. A parallel splitting-up method for partial differential equations and its applications to Navier–Stokes equations. Math Modell Numer Anal 1992; 26: 673–708.

30.

Weickert

Romeny

BTH

Viergever

. Efficient and reliable schemes for nonlinear diffusion filtering. IEEE T Image Process 1998; 7: 398–410.

31.

Persch

Elhayek

Welk

. Enhancing 3-d cell structures in confocal and sted microscopy: a joint model for interpolation, deblurring and anisotropic smoothing. Meas Sci Technol 2013; 24: 125703.

32.

Vogel CR and Oman ME. Fast total variation-based image reconstruction. In: Proc. 1995 ASME design engineering conferences, 1995, volume 3, part C, pp.1009–1015.

33.

Chen

. Matrix preconditioning techniques and applications, Cambridge, UK: Cambridge University Press, 2005.

34.

Crum

Camara

Hill

. Generalized overlap measures for evaluation and validation in medical image analysis. IEEE T Med Imaging 2006; 25: 1451–1461.

35.

Duda

Richard

. Hart, pattern classification and scene analysis, New York: John Wiley and Sons, 1973.

36.

Jaccard

. Nouvelles recherches sur la distribution florale. Bull Soc Vaudoise Sci Nat 1908; 44: 223–270.

37.

Dice

. Measures of the amount of ecologic association between species. Ecology 1945; 26: 297–302.

38.

Almeida

MSC

Almeida

. Blind and semi-blind deblurring of natural images. IEEE T Image Process 2010; 19: 36–52.

39.

Bar

Sochen

Kiryati

. Semi-blind image restoration via Mumford–Shah regularization. IEEE T Image Process 2006; 15: 483–493.

40.

Makni S, Ciuciu P, Idier J, et al. Semi-blind deconvolution of neural impulse response in FMRI using a Gibbs sampling method. In: IEEE Int Conf Acoustics Speech and Signal Processing (ICASSP), IEEE, New York, NY, USA, 2004, volume 5, pp. V–601–604.

41.

Park

Dobigeon

Hero

. Semi-blind sparse image reconstruction with application to MRFM. IEEE T Image Process 2012; 21: 3838–3849.

42.

Pillonetto

Cobelli

. Identifiability of the stochastic semi-blind deconvolution problem for a class of time-invariant linear systems. Automatica 2007; 43: 647–654.

43.

Sarri P, Thomas G, Sekko E, et al. Myopic deconvolution combining Kalman filter and tracking control. In: IEEE int conf acoustics speech and signal processing, IEEE, New York, NY, USA, 1998, volume 3, pp.1833–1836.

An effective variational model for simultaneous reconstruction and segmentation of blurred images

Abstract

Keywords

Introduction

Existing methods

Segmentation of blurred images

Two-stage restoration and segmentation of images with known blur

Segmentation of images corrupted by unknown blur

Two-stage restoration and segmentation of images with unknown blur

Algorithm 1

Algorithm 2

A new joint model for the segmentation of images with unknown blur

Segmentation indicator function ν

Region average intensity values c1, c2

Image function u

Transformed image function ψ

Point spread (blur) function k

Transformed point spread (blur) function ω

Overall algorithm

Algorithm 3

A relaxed model for the segmentation of images with unknown blur

Algorithm 4

Experimental results

Models

Measuring error

Result sets

Conclusions

Footnotes

Acknowledgements

Declaration of conflicting interests

Funding

References

Region average intensity values c₁, c₂