Distributed Deblurring of Large Images of Wide Field-Of-View

TL;DR

This paper introduces a distributed image deblurring algorithm that efficiently processes extremely large images by dividing the task among multiple nodes, maintaining quality comparable to centralized methods while reducing resource requirements.

Contribution

It presents a novel distributed deblurring approach for large images, addressing memory limitations and enabling efficient processing across multiple nodes with maintained image quality.

Findings

Achieves deblurring quality similar to centralized methods.

Enables processing of gigapixel images with distributed computing.

Reduces computational resource requirements for large-scale image deblurring.

Abstract

Image deblurring is an economic way to reduce certain degradations (blur and noise) in acquired images. Thus, it has become essential tool in high resolution imaging in many applications, e.g., astronomy, microscopy or computational photography. In applications such as astronomy and satellite imaging, the size of acquired images can be extremely large (up to gigapixels) covering wide field-of-view suffering from shift-variant blur. Most of the existing image deblurring techniques are designed and implemented to work efficiently on centralized computing system having multiple processors and a shared memory. Thus, the largest image that can be handle is limited by the size of the physical memory available on the system. In this paper, we propose a distributed nonblind image deblurring algorithm in which several connected processing nodes (with reasonable computational resources) process simultaneously different portions of a large image while maintaining certain coherency among them to finally obtain a single crisp image. Unlike the existing centralized techniques, image deblurring in distributed fashion raises several issues. To tackle these issues, we consider certain approximations that trade-offs between the quality of deblurred image and the computational resources required to achieve it. The experimental results show that our algorithm produces the similar quality of images as the existing centralized techniques while allowing distribution, and thus being cost effective for extremely large images.

Equations

$${\mathbf{x}}^{\ast }:=\underset{\mathbf{x}}{\mathrm{argmin}}\{{\Psi }_{\text{data}}(\mathbf{y},\mathbf{x})+\lambda {\Psi }_{\text{prior}}(\mathbf{x})\}$$

$$\mathbf{y}=\mathbf{H}\mathbf{x}+\mathbf{\epsilon }$$

$${\mathbf{x}}^{\ast }:=\underset{\mathbf{x}\ge 0}{\mathrm{argmin}}\{\frac{1}{2}\Vert \mathbf{y}-\mathbf{H}\mathbf{x}{\Vert }_{\mathbf{W}}^2+\lambda \varphi (\mathbf{D}\mathbf{x})\}$$

$$\forall (i,j)\in \{1,\dots ,N{\}}^2,\forall \ell \in {\mathcal{P}}_i\cap {\mathcal{P}}_j,{\mathbf{x}}_i(\ell )={\mathbf{x}}_j(\ell ).$$

$$\underset{{\mathbf{x}}_1\cdots {\mathbf{x}}_N}{\mathrm{argmin}}\sum _{i=1}^N{f}_i({\mathbf{x}}_i)\text{s.t. }({\mathbf{x}}_1,\dots ,{\mathbf{x}}_N)\in \mathcal{C}.$$

$${\mathbf{x}}_i\mapsto \frac{1}{2}\Vert {\mathbf{y}}_i-{\mathbf{H}}_i{\mathbf{x}}_i{\Vert }_{{\mathbf{W}}_i}^2,$$

$$f_i({\mathbf{x}}_i)=\frac{1}{2}\Vert {\mathbf{y}}_i-{\mathbf{H}}_i{\mathbf{x}}_i{\Vert }_{{\mathbf{W}}_i}^2+{\lambda }_i{\varphi }_i({\mathbf{D}}_i{\mathbf{x}}_i)+{\iota }_{{\mathbb{R}}_{+}^{{\mathcal{P}}_i}}\left({\mathbf{x}}_i\right)$$

$${\mathrm{prox}}_{{\mathbf{V}}^{-1},h}(\mathbf{v})=\underset{\mathbf{w}\in \mathcal{X}}{\mathrm{argmin}}h(\mathbf{w})+\frac{\Vert \mathbf{w}-\mathbf{v}{\Vert }_{\mathbf{V}}^2}{2}$$

$$\underset{\mathbf{x}\in \mathcal{X}}{\mathrm{argmin}}f(\mathbf{x})+g(\mathbf{x})$$

$${\mathbf{x}}^{(k+1)}$$

$${\mathbf{z}}^{(k+1)}$$

$${\mathbf{u}}^{(k+1)}$$

$${\mathbf{V}}_i{\mathbf{x}}_i:{\mathcal{P}}_i$$

$$\ell$$

$${\mathbf{z}}_i(\ell )=\frac{{\sum }_{i\in {\mathcal{P}}_{\ell }^{-}}{\mathbf{\alpha }}_i(\ell ){\mathbf{x}}_i(\ell )}{{\sum }_{i\in {\mathcal{P}}_{\ell }^{-}}{\mathbf{\alpha }}_i(\ell )}.$$

$$g(\mathbf{x})=\sum _{\ell =1}^n{\iota }_{{\mathcal{C}}_{\ell }}({\mathbf{x}}_{{\mathcal{P}}_{\ell }^{-}}(\ell )).$$

$$\mathbf{z}$$

$$\displaystyle\qquad\qquad\qquad\qquad\Bigl{.}+\frac{1}{2}\sum_{i\in{\mathcal{P}}^{-}_{\ell}}\boldsymbol{\alpha}_{i}(\ell)(\boldsymbol{w}_{i}(\ell)-\boldsymbol{x}_{i}(\ell))^{2}\Bigr{)}$$

$$\bar{\mathbf{z}}(\ell )$$

$$=\frac{{\sum }_{i\in {\mathcal{P}}_{\ell }^{-}}{\mathbf{\alpha }}_i(\ell ){\mathbf{x}}_i(\ell )}{{\sum }_{i\in {\mathcal{P}}_{\ell }^{-}}{\mathbf{\alpha }}_i(\ell )}.$$

$${\mathbf{x}}_i^{(k+1)}$$

$${\bar{\mathbf{z}}}^{(k+1)}(\ell )$$

$${\mathbf{u}}_i^{(k+1)}(\ell )$$

$$\varphi ({\mathbf{D}}_i{\mathbf{x}}_i)=\{\begin{array}{ll}\frac{1}{2}\Vert {\mathbf{D}}_i{\mathbf{x}}_i{\Vert }_2^2& \Vert {\mathbf{D}}_i{\mathbf{x}}_i{\Vert }_2\le \delta \\ \delta (\Vert {\mathbf{D}}_i{\mathbf{x}}_i{\Vert }_2-\frac{\delta }{2})& \text{otherwise}\end{array}$$

Full text

††publicationid: pubid: xxxx–xxxx/xx$00.00 © 2017 IEEE

I Introduction

In many applications, it is essential to have high resolution images for precise analysis and inferences from them. However, a certain amount of degradations (blur and noise) are unavoidable in many of the imaging systems due to several factors, such as the limited aperture size, the involved medium (both atmosphere and optics), and the imaging sensors. With advances in the computational technologies, the digital image restoration techniques, such as image deblurring, have been proven to be an economic way to enhance resolution, signal-to-noise ratio, and contrast of the acquired images. In this paper we can consider the situation where we assume that the information about blur and the statistics of noise is know a priori, and then the restoration is referred to as nonblind deblurring. In general, image deblurring is an ill-posed inverse problem [1, 2]. In a Bayesian setting it can be expressed as maximum-a-posteriori estimation problem [3], which boils down to the following numerical optimization problem:

[TABLE]

where the vectors $\boldsymbol{y}\in\mathbb{R}^{m}$ and $\boldsymbol{x}\in\mathbb{R}^{n}$ represent the observed (acquired blurry and noisy) image, and the unknown crisp image to be estimated, respectively. Two dimensional (2D) images are represented as column vectors by lexicographically ordering their pixels values. The first term $\varPsi_{\text{data}}$ in (1) is called likelihood or data-fidelity term that depends upon the noise and the image formation model. The second term $\varPsi_{\text{prior}}$ is called a priori or regularizer that imposes any prior knowledge on the unknown image $\boldsymbol{x}$. The scalar parameter $\lambda$ keeps trade-off between the likelihood and a priori terms, which generally depends upon the noise level in the observed image.

Blur in an acquired image is characterized by the impulse response of the imaging system, commonly known as point-spread-function (PSF). As far as a narrow field-of-view is concerned, PSF can be considered constant throughout the field-of-view leading to a shift-invariant blur. The blurring operation in this case is modeled as simple a convolution between the sharp image and the PSF, and can be performed efficiently in Fourier domain. In many cases, however, blur varies throughout the field-of-view due to several causes: relative motion between the camera and the scene; moving objects with respect to the background; variable defocussing of non-planar scenes with some objects located in front or behind the in-focus plane; optical aberrations such as space-variant distortions, vignetting or phase aberrations. In here we consider shift-variant blur across the field-of-view of a planar scene, i.e. all the objects in the scene lie in-focus plane; see [4] for discussion on shift-variant blur in the case of non-planar scene. Two subtly different situations can be considered for shift-variant blur in planar scene: i) a large field-of-view is captured into several pieces by multiple imaging systems each capturing a small portion of the whole scene, and ii) the same scene is captured by a single imaging system having wide field-of-view. In the former case, each imaging system can have slightly different blur from each other operating under different settings. Thus, the captured images can have slightly different blurs from each other resulting into piece-wise constant blur in the whole scene. We will refer to this situation as piece-wise constant shift-variant blur. In the latter case, the blur in the whole field-of-view can vary smoothly, e.g. blur due to optical aberrations. We will refer to this situation as smooth shift-variant blur. The shift-invariant and piece-wise constant shift-variant blurs can be considered as special cases of the smooth shift-variant blur. The blurring operation in the smooth shift-variant case cannot be modeled by a simple convolution, and there does not exists any efficient and straight-forward way to perform the operation. However, there are some fast approximations of smooth shift-variant blur operators proposed in [5, 6, 7]; see [8] for detailed comparison between different approximations.

If we approximate the noise in observed image by a non-stationary white Gaussian noise as considered in [9] (see [10, 11, 12] for more refined noise model), then a discrete image formation model can be written as:

[TABLE]

where the matrix $\mathbf{H}\in\mathbb{R}^{m\times n},\ m<n$, denotes the blur operator111Blur operator is a rectangular matrix since size of observed image is restricted by the physical size of the image sensor, i.e., sensor captures little less than what optics can see., and $\boldsymbol{\varepsilon}$ is the zero-mean non-stationary white Gaussian noise, i.e. $\boldsymbol{\varepsilon}(\ell)\sim\mathcal{N}(0,\boldsymbol{\sigma}^{2}(\ell))$ with $\boldsymbol{\sigma}^{2}(\ell)$ denoting the noise variance at $\ell$th component of the vector $\boldsymbol{y}$. We assume that the raw captured image is preprocessed to yield an image that closely follows the above image formation model (2). The preprocessing may includes the correction of the background and of the flat field, the correction of defective pixels and possibly of its correlated noise, and the scaling of the image in photons.

Before we proceed further, let us introduce some more notations. Hereafter, we will use upper-case bold and lower-case bold letters for denoting matrices (and linear operators), and column vectors, respectively. For every set ${\mathcal{T}}$, we denote by $\arrowvert{\mathcal{T}}\arrowvert$ the cardinality of the set, and by $\mathbb{R}^{{\mathcal{T}}}$ the set of functions on ${\mathcal{T}}\to\mathbb{R}$. Let ${\mathcal{X}}$ represents Euclidean space, then we denote by $\langle\,.\,,\,.\,\rangle$ the standard inner product on ${\mathcal{X}}$, and by $\|\,.\,\|$ the Euclidean norm. Letting $\mathbf{V}$ denote a positive definite linear operator of ${\mathcal{X}}$ onto itself, we use the notation $\langle\boldsymbol{x},\boldsymbol{y}\rangle_{\mathbf{V}}=\langle\boldsymbol{x},\mathbf{V}\boldsymbol{y}\rangle$, and we denote by $\|\,.\,\|_{\mathbf{V}}$ the corresponding norm. When $\mathbf{V}$ is a diagonal operator of the form $(\mathbf{V}\boldsymbol{x}):\ell\mapsto{\boldsymbol{\alpha}}(\ell)\ {\boldsymbol{x}}(\ell)$ we equivalently denote $\|\boldsymbol{x}\|_{\mathbf{V}}^{2}$ by $\|\boldsymbol{x}\|_{\boldsymbol{\alpha}}^{2}=\sum_{\ell}\boldsymbol{\alpha}(\ell)\boldsymbol{x}(\ell)^{2}$, where $\boldsymbol{\alpha}=(\alpha_{1},\alpha_{2},\cdots)$. We also denote $\|\boldsymbol{x}\|_{1,\mathbf{V}}=\sum_{\ell}\boldsymbol{\alpha}(\ell)\ |\boldsymbol{x}(\ell)|$ the weighted $L_{1}$-norm, simply noted $\|\boldsymbol{x}\|_{1}$ when $\mathbf{V}$ is the identity. Moreover, if $\mathbf{L}$ is a linear operator, we denote by $\mathbf{L}^{{}^{\mathrm{T}}}$ its adjoint operator.

Considering the above image formation model (2) with a known PSF, the (nonblind) image deblurring problem (1) can be explicitly expressed as:

[TABLE]

where the matrix $\mathbf{W}$ is diagonal with its components given by $\mathbf{W}(\ell,\ell)=1/\boldsymbol{\sigma}^{2}(\ell)$ for observed pixels and $\mathbf{W}(\ell,\ell)=0$ for unmeasured pixels. The function $\phi$ represents some regularizer (e.g., Tikhonov’s $L_{2}$-norm, sparsity promoting $L_{1}$-norm, or Huber’s mixed $L_{1}$-$L_{2}$-norm), and $\mathbf{D}$ is some linear operator, e.g., finite forward difference or some discrete wavelet transform. The formalism (3) of image deblurring problem is referred to as analysis-based approach where the unknown variable is expressed in the image domain itself. An alternative formalism, referred to as synthesis-based approach, is also considered in literature where the unknown variable is in some transformed domain, e.g., coefficients of some discrete wavelet; see [13] for detailed discussion on the two approaches and their comparisons. Without loss of generality, in this paper we consider only the analysis-based approach for expressing the image deblurring problem.

Depending upon the structures of the two terms in the optimization problem (3), the solution can be obtained in a single step i.e., a closed-form solution (e.g., Wiener deconvolution), or one has to rely on iterative solvers, which may require, at each iteration, the same computational expense or even more than the closed-form solution222Applying blur operator $\mathbf{H}$ or its adjoint $\mathbf{H}^{{}^{\mathrm{T}}}$ is an expensive operation, and one may require to apply them several times per iterations of the solver.. The latter one is the frequently occurring case in many applications, e.g., the regularization term is not differentiable (nonsmooth). Though, there exists a vast literature on numerical optimization, but we refer the readers to [14, 15, 16, 17] and the references therein for some of the recent and fast optimization algorithms specially developed for inverse problems in imaging.

I-A Motivation: Deblurring Extremely Large Image

With the advances in imaging technologies, the applications in astronomy, satellite imagery and others are able to capture extremely large image of wide field-of-view. The size of such images can vary from few tens megapixels up to gigapixels. As discussed previously, we consider two imaging situations where such large images are acquired: piece-wise constant and smooth shift-variant blur. In the former situation, the images from multiple narrow field-of-view imaging systems are supposed to be mosaicked together into a single image provided that their relative positions within the whole field-of-view are known. For further simplification, we assume that all the narrow field-of-view imaging systems have the same noise performance. In the both imaging situations, it is highly desirable to obtain a single crisp image from the acquired image(s) for any further applications.

Image deblurring is well studied topic and a vast literature exists [3, 1, 2, 18, 19, 11, 20, 16] for moderate size (few megapixels) images with shift-invariant blur. These methods differ from each other in the way they formulate the data-fidelity term or the regularization term, and the optimization algorithms they propose to solve the resulting problem. The image deblurring problem becomes more complicated when one considers the images suffering from shift-variant blur. Some recent works [5, 6, 7, 8] proposed efficient methods for deblurring images with shift-variant blur. However, all the methods listed above are designed and implemented to work efficiently on a centralized computing system, possibly with multiple processor cores sharing a single memory space. Such a computing system is commonly referred to as multi-threaded shared memory system based on single instruction multiple data (SIMD) architecture. Hereafter, we will refer to all such existing methods by centralized deblurring techniques. In contrast to deblurring a moderate size image, deblurring an extremely large image is a challenging problem since the largest size of the image that can be handled is then limited by the capacity of shared memory available on the centralized system. It is not cost effective to build a centralized system with several processor cores and a huge shared memory when the modern distributed computing systems (consist of several connected processing nodes each having reasonable amount of computational resources) are proving to be far more economical way for solving huge-scale problems than the centralized system. Many domains such as machine learning, data mining [21], and others are already benefiting from distributed computing approach for discovering patterns in huge datasets distributed at different nodes. As per our knowledge, there are very few works, e.g., [22, 23, 24, 25], considering distributed computing approach for image restoration/reconstruction problems. Within this context, we propose a distributed image deblurring algorithm for large images acquired from either of the imaging situations mentioned above.

I-B Our contributions

In this paper, we propose a distributed algorithm for deblurring large images, which needs less frequent inter-node communication. Our algorithm can handle the acquired image(s) from both the imaging situations: smooth shift-variant and piece-wise constant shift-variant blur. For the first imaging situation, we consider splitting the large image into sufficiently small overlapping blocks and then deblur them simultaneously on the several processing nodes while maintaining certain coherencies among them so as to obtain a single crisp image without any transition artifacts among the deblurred blocks. For the second imaging situation, we assume that the relative positions of narrow images in the whole field-of-view are known, and there are certain overlaps among them. If a narrow image is not sufficiently small to fit in the memory of a node, then we can split it further into smaller overlapping blocks. Similar to the first situation, we deblur the small images simultaneously to obtain a single crisp image. To do so, we reformulate the image deblurring problem (1) into a distributed optimization problem with consensus, and then present an efficient optimization method to solve it. Our distributed deblurring algorithm is rather generic in the sense that it can handle different situations such as shift-invariant or shift-variant blur with any possible combination of the data-fidelity and regularization terms. Depending on the structures of the data-fidelity and regularizer, we can select any fast optimization algorithm for solving the local deblurring problems at processing nodes. By several numerical experiments, we show that our algorithm is cost effective in term of computational resources for large images, and is able obtain the similar quality of deblurred images as the existing centralized deblurring techniques that are applicable only to images of moderate size.

The remaining parts of the paper is organized as follows. In Section II, we discuss the difficulties associated with distributed formulations of image deblurring problem, and the possible approaches to overcome it. In Section III, we present the distributed formulation and an efficient distributed algorithm to solve it. We discuss the criteria for selecting the different parameters associated with our approach. In Section IV, we discuss the implementation details and present several numerical experiments where we compare the results obtained from different approaches. Finally, in Section V, we conclude our work with possible future enhancements.

II Distributed Computing Approach for Image Deblurring

A generic approach for dealing with large-scale problem is the “divide and conquer” strategy, in which a large-scale problem is decomposed into smaller subproblems in such way that solving them and assembling theirs results would produce the expected final result or, at least, reasonably close to it. Distributed computing approach has emerged as a framework of such a strategy. Distributed computing systems are consist of several processing nodes, each having a reasonable amount of computational resources (in term of memory and processor cores), connected together via some communication network so that nodes can exchange messages to achieve a certain common goal. They are built based on the multiple instructions multiple data (MIMD) architecture. Nodes in a distributed system may not be necessarily located at the same physical location, so a high speed communication links may not be always feasible among them. Thus, an efficient distributed algorithm is the one which is computation intensive rather than communication demanding. In many applications such as machine learning, data mining, distributed approaches have become de-facto standard for efficiently estimating extremely large number of parameters using huge dataset distributed on the different nodes. Taking such an inspiration, one can devise a distributed strategy for deblurring extremely large images. A possible approach is to use the distributed array abstraction available on modern distributed systems, and reimplement the standard centralized deblurring techniques using distributed arrays instead of a shared memory array. However, the bottleneck of such an approach would be extensive data communication among the nodes at each iteration of the optimization algorithm. To overcome this, a straightforward approach would be to split the given observed image $\boldsymbol{y}$ into $N$ smaller contiguous blocks $\boldsymbol{y}_{i}\in\mathbb{R}^{m_{i}},i=1,2,\cdots,N$, and deblur them independently to obtain the deblurred blocks $\hat{\boldsymbol{x}}_{i}\in\mathbb{R}^{n_{i}},i=1,\cdots,N$ that can be merged to obtain the single deblurred image $\hat{\boldsymbol{x}}\in\mathbb{R}^{n}$. However, there may arise several issues from both the theoretical and practical point-of-views as discussed below.

II-A Issues with Distributed Approach for Image Deblurring

For sake of simplicity, let us first consider the shift-invariant image deblurring problem, and point out some of the major issues:

Practically, it is infeasible to explicitly create the blur operator matrix $\mathbf{H}$. For this reason blur operation is efficiently performed in Fourier domain using Fast Fourier Transforms (FFT) algorithm. Thus, there is no straight forward way to split the matrix $\mathbf{H}$ into block matrices $\boldsymbol{H}_{i}\in\mathbb{R}^{m_{i}\times n_{i}}$ so that they could be used to perform independent deblurring of the observed blocks $\boldsymbol{y}_{i}$ such that the final result would be equivalent to using the original matrix $\mathbf{H}$ on whole image $\boldsymbol{y}$. Any approximation of block matrices $\boldsymbol{H}_{i}$ chosen not appropriately may create artifacts such as nonsmooth transition among the deblurred blocks. 2. 2.

Blur operation performed in Fourier domain assumes circular boundary condition. The circulant blur assumption in image deblurring process can lead to severe artifacts due to discontinuity at boundaries caused by the periodic extension of image [16]. Thus, deblurring the individual blocks $\boldsymbol{y}_{i}$ can introduce boundary artifacts into them that will eventually worsen the quality of final deblurred image. 3. 3.

Splitting the image into blocks can also raise issues with the regularization term. Many frequently used regularizers are not separable, i.e., its value on the whole image is not equivalent to the sum its value on the blocks. Such an approximation can result into structural incoherencies among the deblurred blocks.

The issues mentioned above complicate further the shift-variant image deblurring problem, since there is no any efficient and straightforward way to formulate shift-variant blur operator. Considering the above issues associated with the naïve “split, deblur independently, and merge” approach, it is obvious that one needs to approximate the problem (3) in a way so that it can be solved efficiently in a distributed manner, yet obtain the solution that remains reasonably close to the solution of the original problem (3). In the next section, we describe how to tackle the above issues, and then we will present our distributed image deblurring algorithm.

II-B Tackling the Issues in Distributed Image Deblurring

As discussed above, it is nontrivial problem to explicitly create a blur operator and then decompose it into block matrices to operate them independently for distributed image deblurring. However, we can formulate an approximation of the block matrices $\boldsymbol{H}_{i}$ such that when they are used in deblurring process, certain homogeneities among the deblurred blocks are imposed. Let us consider the smooth shift-variant blur case. Provided that we are able to sample $N$ local PSFs $\boldsymbol{h}_{i},i=1,\cdots,N$ at regular grid points within the field-of-view, the approximation of shift-variant blur operator proposed in [6, 7] suggests an interesting idea to formulate the block matrices $\boldsymbol{H}_{i}$. Their approximations of the shift-variant blur operator is based upon the idea that a PSF at any point within field-of-view can be well approximated by linear combination (interpolation) of the PSFs sampled at the neighboring grid points. Said so, the shift-variant blurring operation can be written as: $\boldsymbol{y}=\mathbf{R}\sum_{i=1}^{N}\boldsymbol{Z}_{i}\boldsymbol{H}_{i}\boldsymbol{W}_{i}\boldsymbol{C}_{i}\boldsymbol{x}$, where $\boldsymbol{C}_{i}$ are chopping operators that select overlapping blocks from the image $\boldsymbol{x}$, $\boldsymbol{W}_{i}=\operatorname{diag}(\boldsymbol{\omega}_{i})$ are interpolation weights corresponding to each blocks, $\boldsymbol{H}_{i}$ are blur operators corresponding to the sampled local PSFs $\boldsymbol{h}_{i}$, $\boldsymbol{Z}_{i}$ are operators that zero-pads the corresponding blurred blocks keeping it in the same relative position with respect to the whole image, and $\mathbf{R}$ is a chopping operator that restrict the final blurred image to sensor size. Interpolation weights $\boldsymbol{\omega}_{i}$ are non-zero only within certain locality depending upon the order of interpolation. This approximation allows to perform shift-variant blurring combining several local shift-invariant blurring, which can be performed efficiently in Fourier domain. This approximation renders smooth variation of blur throughout the whole image. Higher accuracy in blurring operation can be achieved by denser sampling of PSFs with the field-of-view. The image deblurring results presented in [6, 8] suggest that the first-order interpolation (e.g., $\boldsymbol{\omega}_{i}$ is a ramp within the range $[0,1]$ between the two adjacent grid points for 1D signal as illustrated in Fig.1) is sufficient to achieve reasonably good quality of deblurred images. Using first-order interpolation, the shift-variant blur operator is only four times expensive than the shift-invariant blur operation on the same size of image. However, this fast approximation of shift-variant blur operator is efficient only for a centralized multi-threaded shared memory implementation. It is not directly applicable for distributed setting as it will lead to intensive data communication among the nodes each time the blur operator or its adjoint is used by an iterative solver. Nevertheless, this approximation suggests the following idea for distributed image deblurring: i) split the observed image $\boldsymbol{y}$ into $N$ overlapping blocks $\boldsymbol{y}_{i}$, ii) generate 2D first-order interpolation weights $\boldsymbol{\omega}_{i}$ corresponding to each block, iii) provided with $N$ PSFs $\boldsymbol{h}_{i}$ sampled locally within each block, distribute the observed blocks, the PSFs, and the interpolation weights among the $N$ nodes, and iv) then on each node perform local image deblurring while maintaining certain consensus among the overlapping pixels of the adjacent blocks. The consensus operation can be a weighted averaging among the overlapping pixels. Using the interpolation weights $\boldsymbol{\omega}_{i}$ for averaging operation will indirectly impose the smooth variation of blur across whole observed image. Without loss of generality, the aforementioned strategy is also applicable to the cases when an image suffers from shift-invariant or piece-wise constant shift-variant blur since these two operations can be well approximated by the above fast shift-variant blur operator for the smooth blur variation.

As pointed out above, performing local deblurring in Fourier domain may lead to boundary artifacts, thus to avoid such artifacts in the deblurred blocks, we borrow the idea from [16]. We express the local blur operator as: $\boldsymbol{H}_{i}=\boldsymbol{\mathcal{C}}_{i}\boldsymbol{\mathcal{H}}_{i},i=1,\cdots,N$, where $\boldsymbol{{\mathcal{H}}}_{i}\in\mathbb{R}^{n_{i}\times n_{i}}$ are circular convolution matrices formed from PSFs $\boldsymbol{h}_{i}$, and $\boldsymbol{\mathcal{C}}_{i}\in\mathbb{R}^{m_{i}\times n_{i}}$ are chopping operators, which restrict the convolution results to the valid regions. Expressing the local blur operators in this form allows an efficient non-circulant blur operation in Fourier domain while suppressing boundary artifacts in the deblurred patches.

Concerning the last issue about approximating regularizer on the whole image by sum of it on blocks, it depends upon the structure of regularizer. If the regularizer is separable in image domain, e.g., $\phi(\mathbf{D}\boldsymbol{x})=\|\boldsymbol{x}\|_{1}\ \text{or}\ \|\boldsymbol{x}\|^{2}_{2}$, then the aforesaid approximation holds. However, for many of the frequently used regularizers, e.g., $\|\mathbf{D}\boldsymbol{x}\|^{2}_{2}$ or $\|\mathbf{D}\boldsymbol{x}\|_{1}$ or $\|\mathbf{D}\boldsymbol{x}\|_{2}$ where $\mathbf{D}$ represents finite-forward difference operator or some discrete wavelet transform, the aforementioned approximation does not hold. This can lead to incoherencies among the deblurred blocks, which can render nonsmooth transitions among them. However, as shown hereafter a sufficient amount of overlaps among the adjacent blocks with aforementioned consensus imposed on the overlapping pixels during deblurring process will limit the incoherencies, and suppress any nonsmooth transition in the final deblurred image.

III The Proposed Approach

Considering all the ideas developed in Section II, we propose a generic framework for distributed image deblurring applicable to images suffering from both the shift-invariant and shift-variant blur. We make some reasonable approximations, and reformulate the original problem (3) into the distributed optimization problem presented below. The resulting optimization problem is solved in a distributed manner by Douglas-Rachford (D-R) splitting algorithm that first appeared in [26].

III-A General Setting

Consider a distributed computing system with a set of $N$ nodes having peer-to-peer bidirectional connections among them, or at least connection topology shown in Fig. 2. Given an observed image $\boldsymbol{y}\in\mathbb{R}^{m}$, we split it into $N$ blocks $\boldsymbol{y}_{i}\in\mathbb{R}^{m_{i}},\forall i=1,2,\cdots,N$, with certain overlaps among them. We generate 2D first-order interpolation weights $\boldsymbol{\omega}_{i}$ corresponding to the observed blocks as shown in Fig. 3(e). Provided that we are able to sample $N$ PFSs $\boldsymbol{h}_{i}$ within the blocks, and $\boldsymbol{H}_{i}\in\mathbb{R}^{m_{i}\times n_{i}}$ be the corresponding blur operators, we distribute the observed blocks, the corresponding PSFs, and the interpolation weights among the nodes. We, then, seek to distributively estimate the whole unknown crisp image $\boldsymbol{x}\in\mathbb{R}^{n}$. Let us denote by ${\mathcal{P}}_{1},\dots,{\mathcal{P}}_{N}$ a collection of $N$ subsets of $\{1,\dots,n\}$. For every $i=1,\dots,N$, we assume that $i$th compute node is in charge of estimating the components of $\boldsymbol{x}$ corresponding to the indices ${\mathcal{P}}_{i}$. The subsets ${\mathcal{P}}_{1},\dots,{\mathcal{P}}_{N}$ are overlapping. Hence, different nodes handling a common component of $\boldsymbol{x}$ must eventually agree on the value of the latter. Formally, we introduce the product space ${\mathcal{X}}:=\mathbb{R}^{{\mathcal{P}}_{1}}\times\cdots\times\mathbb{R}^{{\mathcal{P}}_{N}}$, and we denote by ${\mathcal{C}}$ the set of vectors $(\boldsymbol{x}_{1},\dots,\boldsymbol{x}_{N})\in{\mathcal{X}}$ satisfying the restricted consensus condition

[TABLE]

Moreover, we assume that every $i$th node is provided with a local convex, proper and lower semicontinuous function $f_{i}:\mathbb{R}^{{\mathcal{P}}_{i}}\to(-\infty,+\infty]$. We consider the following constrained minimization problem on $\mathbb{R}^{{\mathcal{P}}_{1}}\times\cdots\times\mathbb{R}^{{\mathcal{P}}_{N}}$:

[TABLE]

For our image deblurring problem, the local cost function $f_{i}$ is composed of the local data-fidelity term

[TABLE]

for some positive definite $\mathbf{W}_{i}(\ell,\ell)=1/{\boldsymbol{\sigma}_{i}^{2}(\ell)}$ as in (3), and a regularizer $\phi_{i}(\mathbf{D}_{i}\ \boldsymbol{x}_{i})$ with positivity constraint on $\boldsymbol{x}_{i}$, i.e., $\boldsymbol{x}_{i}\in\mathbb{R}^{{\mathcal{P}}_{i}}_{+}$. If ${\mathcal{A}}$ is a set, the notation $\iota_{{\mathcal{A}}}$ stands for the indicator function of the set ${\mathcal{A}}$, equal to zero on that set and to $+\infty$ elsewhere. Thus, the local cost function needed to be minimized at each of the nodes have the form:

[TABLE]

where $\phi_{i}$ are convex, proper and lower semicontinuous functions and $\mathbf{D}_{i}$ are linear operators on $\mathbb{R}^{{\mathcal{P}}_{i}}$.

III-B Optimization Algorithm

Before we present our distributed optimization algorithm for solving problem (4), let us introduce one more notation. For any convex, proper and lower semicontinuous function $h:{\mathcal{X}}\to(-\infty,+\infty]$, we introduce the proximity operator

[TABLE]

for every $\boldsymbol{v}\in{\mathcal{X}}$.

For solving (4), we consider D-R Splitting algorithm, thus we reformulate (4) as

[TABLE]

where $g=\iota_{{\mathcal{C}}}$ is the indicator function of ${\mathcal{C}}$ and $f(\boldsymbol{x})=\sum_{i}f_{i}(\boldsymbol{x}_{i})$ for every $\boldsymbol{x}=(\boldsymbol{x}_{1},\dots,\boldsymbol{x}_{N})$ in ${\mathcal{X}}$. Let us equip the Euclidean space ${\mathcal{X}}$ with the inner product $\langle\,.\,,\,.\,\rangle_{\mathbf{V}}$ for some positive definite linear operator $\mathbf{V}:{\mathcal{X}}\to{\mathcal{X}}$. Let $\rho^{(k)}$ be a sequence in $]0,2[$, and $\boldsymbol{\epsilon}^{(k)}_{1}$ and $\boldsymbol{\epsilon}^{(k)}_{2}$ be sequences in ${\mathcal{X}}$, then the D-R splitting algorithm writes as:

[TABLE]

If the following holds:

the set of minimizers of (4) is non-empty, 2. 2.

$\boldsymbol{0}\in\mathrm{ri}(\mathrm{dom}(f)-{\mathcal{C}})$, where $\mathrm{ri}$ represents relative interior, 3. 3.

$\sum_{k}\rho^{(k)}(2-\rho^{k})=+\infty$, 4. 4.

and $\sum_{k}\|\boldsymbol{\epsilon}^{(k)}_{1}\|_{2}+\|\boldsymbol{\epsilon}^{(k)}_{2}\|_{2}<+\infty$,

then the iterates $\boldsymbol{u}^{(k)}$ converges weakly to some point in ${\mathcal{X}}$ and $\boldsymbol{x}^{(k)}$ converge to a minimizer of (4) as $k\to\infty$; see [27, Corollary 5.2] for the proof. The parameter $\rho^{(k)}$ is referred to as a relaxation factor, which can be tuned to improve the convergence. The sequences $\boldsymbol{\epsilon}^{(k)}_{1}$ and $\boldsymbol{\epsilon}^{(k)}_{2}$ allow some perturbations in the two $\operatorname{prox}$ operations, which is very useful in the cases when the $\operatorname{prox}$ operations do not have closed-form solutions and have to rely on some iterative solvers.

From now onward, we assume that $\mathbf{V}$ is a diagonal operator of the form $\mathbf{V}\boldsymbol{x}=(\mathbf{V}_{1}\boldsymbol{x}_{1},\dots,\mathbf{V}_{N}\boldsymbol{x}_{N})$, where for every $i$,

[TABLE]

where $\boldsymbol{\alpha}_{i}(\ell)$ is a positive coefficient to be specified later. For every $\ell\in\{1,\dots,n\}$, we introduce the set ${\mathcal{P}}_{\ell}^{-}=\{i\,:\,\ell\in{\mathcal{P}}_{i}\}$.

Lemma III.1**.**

For every $\boldsymbol{x}\in{\mathcal{X}}$, the quantity $\boldsymbol{z}=\operatorname{prox}_{\mathbf{V}^{-1},\iota_{{\mathcal{C}}}}(\boldsymbol{x})$ is such that for every $i\in\{1,\dots,N\}$ and every $\ell\in{\mathcal{P}}_{i}$,

[TABLE]

Proof.

For every $\boldsymbol{x}\in{\mathcal{X}}$ of the form $\boldsymbol{x}=(\boldsymbol{x}_{1},\dots,\boldsymbol{x}_{N})$, and every $\ell\in\{1,\dots,n\}$, we use the notation $\boldsymbol{x}_{{\mathcal{P}}^{-}_{\ell}}(\ell)=(\boldsymbol{x}_{i}(\ell)\,:\,i\in{\mathcal{P}}^{-}_{\ell})$. We denote by ${\mathcal{C}}_{\ell}$ the linear span of the vector $(1,\dots,1)$ in ${\mathbb{R}}^{|{\mathcal{P}}^{-}_{\ell}|}$, where $|{\mathcal{P}}^{-}_{\ell}|$ is the cardinality of ${\mathcal{P}}^{-}_{\ell}$. The function $g=\iota_{\mathcal{C}}$ writes

[TABLE]

For an arbitrary $\boldsymbol{x}\in{\mathcal{X}}$, the quantity $\boldsymbol{z}=\operatorname{prox}_{\mathbf{V}^{-},g}(\boldsymbol{x})$ is given by

[TABLE]

Clearly, for every $\ell\in\{1,\dots,n\}$, the components $\boldsymbol{z}_{i}(\ell)$ are equal for all $i\in{\mathcal{P}}^{-}_{\ell}$, to some constant w.r.t. $i$, say $\bar{\boldsymbol{z}}(\ell)$. The later quantity is given by

[TABLE]

∎

By the above lemma, the D-R splitting algorithm can be explicitly written as follows. For every $i\in\{1,\dots,N\}$ and every $\ell\in\{1,\dots,n\}$,

[TABLE]

The resulting algorithm is a synchronous distributed algorithm without any explicit master node. It is depicted in Algorithm 1. Hereafter, we will refer to it by proposed deblurring method. Given some initial guess of the local solutions at each node, the first step of Algorithm 1 execute the Local-Solver in parallel on all the nodes to obtain local deblurred blocks. Then all the nodes synchronize at the start of the second step, and then they exchange the overlapping pixels with its adjacent nodes to distributively perform the weighted averaging of those pixels. Then, the last step of the algorithm is executed in parallel on all the nodes.

The convergence speed of proposed deblurring method is dependent upon the parameters $\boldsymbol{\alpha}_{i}$, and like in other optimization algorithms e.g., ADMM [28], selecting the optimal values of $\boldsymbol{\alpha}_{i}$ for fast convergence is a tedious task. We select $\boldsymbol{\alpha}_{i}=\gamma\;\boldsymbol{\omega}_{i}$ for $\gamma>0$, so that we can tune $\gamma$ for fast convergence, and as well impose the smooth variation of the blur among adjacent blocks.

III-C Size of Image Blocks and the Extent of the Overlaps

The size of observed image blocks and the extent of overlaps among them have impact upon the computational cost and the quality of final deblurred image. Thus, they must be chosen appropriately depending upon the situations: shift-invariant, piece-wise constant shift-variant, and smooth shift-variant blur. Without loss of generality, we assume that all local PSFs have the same size.

In the case of smooth shift-variant blur, the size of image blocks and the extent of overlaps depend upon how densely the grid points, thus the local PSFs, can be sampled within the field-of-view. Denser the grid points, better is the approximation of shift-variant blur in the image. If we consider $N$ equidistant grid points within the field-of-view, then the support of image blocks extend (in both horizontal and vertical directions) over three consecutive grid points, except for the image blocks at the boundaries of image. Thus, each block overlaps half of its size from top, bottom, left and right with its adjacent blocks, i.e., a pixel in the observed image appears into four adjacent blocks. As discussed previously, this extent of overlaps are necessary to impose smooth variation of blur among the image blocks. See the illustration in Fig. 1, which demonstrate for 1D signal the scheme of splitting the observed image into blocks, and the extent of overlaps among them. Similarly, Fig. 5(a) shows $5\times 5$ grid points overlaid upon the observed image where the local PSFs are sampled, and Fig. 5(c) shows resulting overlapping image blocks.

In the case of shift-invariant blur, again we can consider $N$ equidistant grid points within the field-of-view, and split the observed image into $N$ overlapping blocks. However, in this case, we do not need the same extent of overlaps among the blocks as in the case of smooth shift-variant blur, but a smaller extent of overlaps. In fact, it is sufficient to have extent of overlaps slightly larger than half the size of the local PSF. In this case, overlaps among the blocks are required not to impose smooth variation of blur (as the PSF same over whole field-of-view), but to minimize any discrepancies in the deblurred image due to approximation introduced in the regularization term. Similarly, in the case of piece-wise constant shift-variant blur, we can consider the same splitting scheme as in the case of shift-invariant blur for the same reason.

In the case of proposed deblurring method, the maximum number blocks we can split the image is dependent upon the number of available nodes to process in parallel. In such a situation, the larger the number of available nodes, the larger the number of blocks we can chose, thus, smaller the size of blocks. From a computational point of view, the smaller the size of the blocks, the smaller the memory and computation time required by nodes to execute the Local-Solver333The size of vectors affects the efficiency of Fast Fourier Transform algorithm and other operations. Smaller vector fit well into the different levels of data-caches in processor so that it consumes lesser clock cycles to execute the operations compared to the situation when the vectors are larger to fit into the data-caches leading to many cache misses and eventually requiring several cycles to fetch them from main memory.. Moreover, the smaller the size of blocks, the smaller is the extent of overlaps among them, thus, the smaller the amount of data to be exchanged among the nodes. With all these advantages from a computational point-of-view, we may intend to split the observed image into as many blocks as the number of available nodes, given that we are provided with as many local PSFs within the field of view. However, the size of the local PSF suggests a lower limit on the size of blocks. We should not select the size of blocks smaller than the size of the PSF, otherwise the local deblurring problem at each node will be highly underdetermined (more number of unknown pixels at boundaries of the blocks need to be estimated), which would eventually deteriorate the quality of the final deblurred image. Thus, to achieve better quality deblurred image, we should select the size of blocks at least twice the size of PSF.

IV Numerical Experiments and Results

Our algorithm depicted in Algorithm 1 can be seen as general framework for distributed image deblurring. Depending upon the situation, we can derive a particular instance from it, e.g., shift-invariant or shift-variant deblurring, and depending upon the structures of data-fidelity and regularizer, we can select any efficient optimization method as the Local-Solver. In order to validate the proposed deblurring method, we performed two numerical experiments, first, considering the simpler case of shift-invariant deblurring, and then more difficult case of the shift-variant deblurring.

Moreover, to evaluate the performance of the proposed deblurring method in terms of quality of deblurred image, we compared it with two other possible approaches. First obvious choice was centralized deblurring as a reference method that solves the original problem (3). Since centralized deblurring method would require a huge memory for a large image, thus we selected images of reasonable sizes for our experiments. We chose the naïve “split, deblur independently, and merge” approach as the second method. In oder to minimize incoherencies among the deblurred blocks obtained by the second method, we split the observed image into overlapping blocks, and blended the locally deblurred blocks into single image by weighted averaging of the overlapping regions by using the same 2D first-order interpolation weights as in proposed method. Hereafter, we will refer to this approach as independent deblurring method. The independent deblurring method is not intended to solve correctly either of the optimization problems (3) or (4), but it is a straight forward way to deblur large images with minimal computational resources. Similarly, in the case of shift-variant image deblurring, we consider comparing our algorithm with the centralized and independent deblurring methods. For the centralized deblurring method, we used the fast shift-variant blur operator based on first-order PSF interpolation as described in [8].

Let us first present the general settings of our experiments, and later more specific details. For the case of shift-invariant deblurring, we considered “Lena” image, which we resized to $1024\times 1024$ pixels, and extended its dynamic range linearly to have maximum intensity up to 6000 photons/pixels. We refer to this image as the reference image. The reference image was blurred with an Airy disk PSF (of size $201\times 201$ pixels) formed due to a circular aperture of radius 6 pixels, and was corrupted with white Gaussian noise of variance $\sigma^{2}=400$ photons/pixels to obtain the observed image. To study the impact of factors discussed in Section III-C, we considered different cases by varying number of blocks and the extent of overlaps among them. Figure 3 shows an instance when observed image is split into $3\times 3$ blocks with overlaps of $100\times 100$ pixels among the adjacent blocks.

For the shift-variant deblurring case, we considered “Barbara” image, which was resized it to $1151\times 1407$ pixels, and its dynamic range was extended linearly to have maximum intensity up to 6000 photons/pixels. As above, we call this image as reference image. We generated $9\times 9$ normalized Gaussian PSFs each of size $201\times 201$ pixels with the central PSF having full-width-half-maximum (FWHM) $=3.5\times 3.5$ pixels, and linearly increasing FWHM in the radial direction up to $16.5\times 10.5$ pixels for the PSF at extreme corner of the reference image as depicted in Fig. 4. These PSFs, somehow, mimic the shift-variant blur due to optical aberrations (coma) [29]. We blurred the reference image with this shift-variant PSFs using shift-variant blur operator based on PSF interpolation described in [8]. We obtained the final observed image by adding white Gaussian noise of variance $\sigma^{2}=400$ photons/pixels to the blurred image. As in shift-invariant deblurring case, we considered different scenarios by varying the number of blocks. Figure 5 shows an experimental setup when using only $5\times 5$ gird of PSFs, i.e., when observed image split is into $5\times 5$ overlapping blocks.

In both the cases, we deliberately selected low level noise in the observed image so that any incoherency or artifact arising in the deblurred image due to the approximations we made was not superseded by the strong regularization level required for low signal-to-noise ratio in observed image. Also, we selected sufficiently large size of PSFs so that they are not band-limited, which is the case in many real imaging systems.

IV-A Choice of Regularization Term

For our experiments, we selected the regularization function $\phi$ to be Huber loss and $\mathbf{D}$ to be circular forward finite difference operator, so the regularizer is written as

[TABLE]

We chose this regularizer for two reasons: i) it is smooth (differentiable), which makes the functions $f_{i}$ in (4) smooth so that we can choose any fast optimization algorithm such as quasi-Newton methods (e.g., BFGS class methods [30]) as the Local-Solver, and ii) it behaves in between Tikhonov and the total-variation regularization, depending upon the value of $\delta$. Thus, it is able to preserve the sharp structures in the images while avoiding stair-case artifacts in smoothly varying regions usually rendered by total-variation regularizer [9]. However, our proposed deblurring method is not restricted to any particular choice of regularizer; depending upon the application one can chose any regularizer and any fast optimization algorithm for local deblurring problem.

IV-B Implementation Details

Since the optimization problems arising in all deblurring methods were smooth, thus we chose a quasi-Newton method called limited-memory variable metric with bound-constraint (VMLM-B)444An open source implementation of VMLM-B in C programming is available at https://github.com/emmt/OptimPack. [31]. VMLM-B is a variant of the standard limited-memory BFGS with bound-constraint (LBFGS-B) [32], but with some subtle differences. VMLM-B does not require any manual tunning of parameters to achieve fast convergence (the only parameter step-length at each iteration of VMLM-B is estimated by line-search methods satisfying Wolfe’s conditions). Thus, this left us with only a single parameter $\gamma$ to be tunned to achieve faster convergence of our algorithm. Moreover, using the same algorithm (VMLM-B) for solving the optimization problems in all three deblurring methods ensured a fair comparison among them in terms of quality of deblurred images and the computational expenses.

In all our experiments, we set $\rho^{(k)}=1$. After few trials, we found that $\gamma=0.001$ results in fast convergence of our algorithm. Since D-R splitting algorithm converges even when the $\operatorname{prox}$ operations are carried out inexactly, thus to speed-up our algorithm, we allowed the Local-Solver to be less accurate during the initial iterations, and then more accurate as the iterations progress [28, Section 3.4.4]. To do so, we allowed VMLM-B to perform 10 inner iterations at the beginning, and increased it by 10 at every next iteration of main loop. To speed-up further, we did warm-start of the VMLM-B at every iterations of our algorithm by supplying $\boldsymbol{x}^{(k)}$ as an initial guess for next $\boldsymbol{x}^{(k+1)}$ estimation. All the results presented below were obtained after 25 iterations of the proposed algorithm. It was observed that 25 iterations were generally sufficient, and more than 25 iterations did not bring any noticeable difference in the results. For the centralized and independent deblurring methods, we allowed VMLM-B to perform a maximum of 1000 iterations or until it satisfied its own stopping criterion based on gradient convergence or progress in the line-search method. It was noticed that, generally, 500 to 600 iterations were sufficient and further iterations did not bring any noticeable difference.

All the three deblurring methods including the centralized, independent and the proposed were implemented in high-level dynamic programming language “Julia”555http://julialang.org/ [33]. The proposed distributed deblurring is implemented using Message Passing Interface (MPI) based on Open MPI library666https://www.open-mpi.org/. The source codes for all the demonstrations shown in this paper will be freely available at https://github.com/mouryarahul/ImageReconstruction.

IV-C Results and Discussion

As mentioned above, we conducted two different experiments, one for shift-invariant, and another for smooth shift-variant blurred image. To compare the quality of the deblurred images obtained from all the three methods we considered two image quality metrics: signal-to-noise ratio (SNR), and Structural Similarity Index (SSIM) [34]. In all the cases, we heuristically fixed the regularization parameter $\delta=100$, and performed deblurring for different values of $\lambda$ in a sufficiently large range to see if one of the methods produces better quality image for certain range of $\lambda$ than others. The plots (SNR vs $\lambda$ and SSIM vs $\lambda$) in the Fig. 3(i–j) show the results obtained from shift-invariant deblurring when observed image was split into $3\times 3$ blocks. To study the influence of size of blocks and extent of overlaps on the quality of deblurred image, we conducted several trails with different settings. The plots in the Fig. 3(k–l) show the influence of extent of overlaps on the image quality. The plots shows that our algorithm performed slightly better than the centralized deblurring in terms of both the SNR and SSIM. This could be due to the approximation introduced in the regularization term in the proposed deblurring method, e.g., sum of the regularizer on the blocks may be favorable than the regularizer on whole for some images depending upon the contents in the image. Moreover, the proposed deblurring may also be benefited by the explicit overlaps among the adjacent blocks; each node can result into slightly different values of the overlapping pixels, which eventually leads to better estimate of those pixels after averaging the estimates from different nodes. This is indicated by the fact that there is slight increase in both the SNR and SSIM with increase in the extent of overlaps as seen in the plots Fig. 3(k–l). It also indicates that an overlap equal to half the size of PSF is sufficient, and a larger overlap did not produce much improvement, but, of course, it did increase the computational and communication cost. We also noticed that the SNR and SSIM of deblurred image obtained from proposed deblurring is slightly less dependent upon the extent of overlaps compared to that obtained by independent deblurring.

The plots in Fig. 6(a–b) compare the results obtained by the three deblurring methods for the case of smooth shift-variant blurred image. First, we noticed that the image quality obtained by all the three methods improves drastically with the increase in density of grid of PSFs sampled within field-of-view. The SNR and SSIM of deblurred image obtained using only $3\times 3$ PSFs is significantly lower than the cases when $5\times 5$, $6\times 6$ and $8\times 8$ grids of PSFs are used. This is due to the fact that the observed image was simulated using a finer $9\times 9$ grids of PSFs, and the coarser grids of PSFs are less accurate in capturing smooth variation of blur than a finer grid of PSFs. We also noticed that unlike the shift-invariant deblurring case, the values of SNR and SSIM obtained by the proposed deblurring is slightly lower than that obtained by the centralize deblurring. This could be due to the fact that both the centralized and proposed deblurring benefits from the explicit overlaps among the block, but some information is always lost at the boundaries of the deblurred blocks in the case of proposed deblurring, which is not the case for centralized deblurring. Also, there must be some influence of the approximation introduced into the regularization term for the case of proposed deblurring. We observed from the plots in Fig. 3(i–j) and Fig. 6(a–b) that three methods do not attained the highest SNR or SSIM at the same value of the regularization parameter $\lambda$.

As expected, we observed that the naïve independent deblurring method performed significantly lower than the other two methods. As pointed out above, this is due to the fact that the method is not intended to solve correctly either the original problem (3) or the distributed formulation (4), but it is the crudest and computationally cheapest way to perform image deblurring by splitting image into pieces.

In our experiments, we observed that the first step, i.e., the Local-Solver, of the proposed deblurring algorithm is the one which consumed the significant computation time among the three steps; Local-Solver took 600 to 800 times more computation time than the consensus step (including communication time among the nodes). This should be true, in general, for many other local deblurring algorithms devised using different combination of data-fidelity and regularization terms depending upon the applications. Thus, the proposed deblurring algorithm is efficient in the sense that it is computation intensive rather than being communication intensive.

For the small or moderate size of images as we considered in our experiments, the proposed deblurring algorithm is computationally more expensive than the centralized and independent deblurring methods; it takes at least 10 to 15 times more computation time than the latter ones. However, in the case of extremely large images for which centralized deblurring is practically not feasible, computational expenses of proposed deblurring is justified by the better quality of deblurred image, which cannot be achieved by the computationally cheaper independent deblurring. Some more simulations performed with different set of images and PSFs are presented in Appendix References, and the results suggest similar conclusions as above.

V Conclusion

In this paper we have proposed a distributed image deblurring algorithm for large images for which the existing centralized deblurring methods are practically inapplicable due to requirement of huge physical memory on a single system. The proposed algorithm is rather a generic framework for distributed image deblurring for it can handle different imaging situations such as deblurring a single large image suffering from shift-invariant blur, a wide field-of-view captured into multiple narrow images by different imaging systems with slightly different PSFs, and a large image of a wide field-of-view suffering from smoothly varying blur captured by a single imaging system. Depending upon the application, one can easily adapt it to include different data-fidelity and regularization terms, and then select any fast optimization algorithm to solve the local deblurring problem at the different nodes of a distributed computing system. Our algorithm is efficient in the sense that it is computation intensive rather than being communication intensive. We showed by experimental results on simulated observed images that the proposed deblurring algorithm produce almost similar or little lower quality (measured in term of SNR and SSIM) of deblurred images than that obtained by centralized deblurring. But, this small compromise in the quality of deblurred image is trade-off by the cost effectiveness of our distributed approach for large images, which is practically not feasible for the centralized deblurring methods. Moreover, we compared the proposed deblurring to a naïve and computationally cheaper independent deblurring, and showed the latter always performed significantly lower than the former. Thus, when high accuracy is desirable, e.g., in astronomical application, the proposed deblurring should be preferred over independent deblurring method, of course at expense of extra computational cost.

In this paper, we considered nonblind image deblurring, i.e., the PSFs were known a priori, however, in real imaging scenario calibrating PSFs accurately is a tedious and challenging task. A more practical way would be to follow a blind image deblurring approach which is able estimate the PSFs from the observed image(s), and recover the single crisp image. Thus, the next perspective step would be to extend the proposed algorithm toward distributed blind image deblurring method which should be able to estimate distributively the PSFs imposing certain regularity among them, and eventually estimate the unknown crisp image.