1 Introduction
Sharpening of blurred images by deconvolution has been a longstanding object of research efforts, and a variety of algorithms for this illposed inverse problem is available which differ widely in their generality, restoration quality and computational expense, with articulate tradeoffs between restoration quality and computational cost. In the case of spaceinvariant blur, where each point of the unobservable (latent) sharp image is smeared out in the image plane in the same way, a typical blur model reads as
(1) 
where the function denotes the observed unsharp image, a pointspread function (PSF) that acts here via convolution , and some additive noise. Dependent on the imaging process, the latter may be replaced by other types of noise such as Poisson or impulse noise.
A socalled nonblind deconvolution problem in which the blurred image and the pointspread function are available as inputs to an algorithm consists then in determining an image such that
(2) 
Methods to solve this task range from the fast and simple Wiener filter [17] which is essentially a regularised inversion of the convolution via the Fourier domain, up to expensive iterative methods [1, 4, 7, 8, 11, 13, 14, 15]. Some of these algorithms require computation via the Fourier domain, while others can be implemented in the spatial domain. Only the latter ones can be extended straightforwardly to the more general setting of spatially variant blur.
Except for the Wiener filter, all commonly used algorithms are iterative. The computational cost of these algorithms is dominated by convolution operations or Fourier transforms that need to be performed in each iteration step, and the total number of iterations needed to achieve satisfactory restoration quality. For many, though not all, of the iterative methods, the expensive step in each iteration is a convolution of the current approximation
with the PSF . This applies, e.g., to total variation deconvolution [3, 13] and similar variational methods [1, 18], RichardsonLucy deconvolution [8, 11] and methods derived thereof [4, 15], but not, e.g., to the methods from [7, 14].Goal and contributions.
We aim at raising the efficiency of iterative deconvolution algorithms in which convolution with is carried out in each iteration step. To this end, we consider algorithmic improvements that speed up the computation of convolutions in the spatial domain, without the use of Fourier transforms. This is interesting for two reasons: First, dramatic speedups can be achieved by convolution implementations tailored to specific PSFs. This has been demonstrated in a recent paper [16] for the particularly favourable case of linear motion blur aligned with an axis direction of the sampling grid, in which a fast boxfilter algorithm clearly outperforms Fourier convolutions. Second, unlike Fourierbased algorithms, spatial domain convolution can be generalised to spatially variant blur settings.
As our first approach, we investigate efficient algorithms for spatial convolution with specific types of spaceinvariant PSFs that often occur in practical situations. Thereby, we generalize the approach of [16]. As our second approach, we address the possibility to restrict convolution to subsets of image pixels in order to save computation cost on those pixels which do not change anymore.
By describing our work as algorithmic optimisations, we follow the widespread use of this term in computer science for algorithmic modifications that allow to perform a computation in a more resourceefficient way (although no optimum in the mathematical sense is usually achieved), see e.g. [5] or for the slightly more general term program optimisation (which also includes code optimisation) [12, p. 84].
Related work.
Efficient box filtering goes back to McDonnell [9], but seems to have received little attention in the image deconvolution context. Efficient convolution techniques, especially for kernels with uniform intensity, have also been investigated in [10]. For efforts to accelerate RL deconvolution see e.g. [2, 6].
Structure of the Paper.
We recall RichardsonLucy deconvolution as the underlying iterative method in Section 2. The first approach to algorithmic improvements, addressing specifically structured convolution kernels, is the topic of Section 3. The second approach, restricting convolution computation to subsets, is considered in Section 4. The achievable speedups are demonstrated by experiments in Section 5, followed by conclusions in Section 6.
2 RichardsonLucy Deconvolution
RichardsonLucy deconvolution (RL) [8, 11] is widely used because of its simplicity and favourable results for moderate noise. Starting from , it uses the iteration
(3) 
to create a sequence of images of increasing sharpness. Here, denotes the adjoint PSF, geometrically obtained by reflecting about the origin. The number of iterations until stopping is the single parameter of the method.
Although we present only tests with RL here, our algorithmic improvements can be applied equally to modified RL methods with additional regularisers or robustified data terms [4, 15], with similar results (compare [16] for linear motion blur).
Boundary treatment.
Convolution with (spatially invariant) PSF always transfers information across image boundaries. In turn, missing information in deconvolution typically leads to artifacts propagating from the boundary into the interior of the image. There are two basic ways to handle these boundary problems in computing convolutions. Either, convolution is computed on a reduced domain for which all input data are available. This, however, is not an option in iterative procedures. Alternatively, the image has to be continued (padded). Convolution via discrete Fourier transforms implicitly introduces periodic continuation. In spatial domain convolution, an often reasonable compromise between suppression of artifacts and computational effort is to replace every missing input pixel (outside the image bounds) with the closest image pixel, which means a constant continuation along lines perpendicular to the image boundary. We will implement this boundary treatment in all of our convolution procedures acting in the spatial domain.
3 Fast Convolution with Special PointSpread Functions
(a)  (b)  (c)  (d)  (e) 
Linear motion blur.
The simplest case we mention is linear motion blur in scanline () direction, see Figure 1(a). This particular PSF has already been considered in [16]. Here, convolution can efficiently be implemented using a fast box filter [9] that acts in each scanline via a sliding window. For an image and an pixel PSF, the complexity of this filter is , as contrasted to for naive spatial convolution or for an FFTbased method.
The same result, with basically the same complexity, can be obtained as follows: one computes first in each scanline of a discrete image an array of cumulated sums ; then each pixel of the convolution result is given by subtracting two array entries such as , compare Fig. 2(a). Each of the two steps has linear complexity in the number of pixels; however, to implement the desired boundary conditions, the image must be extended by the PSF size in direction, implying again overall complexity.
Rectangular blur.
As our next test case we consider a 2dimensional PSF consisting of a rectangle with constant density aligned with the scanlines, see Fig. 1(b). While this PSF type is of less practical importance, it can serve as an intermediate step towards the more realistic PSFs considered later on. A straightforward transfer of the sliding window idea uses the separability of this PSF: an intermediate image is computed by the boxfilter method in direction, followed by another boxfilter step in direction. The overall complexity of this approach for an image and an pixel PSF is .
An alternative generalisation starts from the abovementioned cumulated sum approach. In a first step, a cumulated sum array is computed, in which contains the sum of all greyvalues left and above position . This step has linear complexity in the number of pixels. In a second step, one computes each pixel of via , compare Fig. 2(b). Again, for the desired boundary conditions, the cumulated image must be computed in size , yielding the same overall complexity as above. One checks easily, however, that less operations per pixel are needed.
Defocus blur.
A frequent source of blur in image acquisition is defocussing which, for a circularshaped camera aperture, leads to a discshaped point spread with constant density. To apply the sliding window approach, note that shifting by one pixel to the right removes from the mask just one (no longer straight) line of pixels at the left boundary while adding one line at the right boundary. If the PSF is pixels high, this update requires operations. For a PSF enclosed in an box, the complete sliding window summation thus takes time. This algorithm can be applied to any convex shape with uniform intensity. We will refer to it as generic box filter. The cumulated sum approach does not reduce complexity in this situation.
Sparse blur.
If images are degraded by motion blur not aligned to the sensor grid, or irregular motion e.g. due to camera shake, representing the pointspread function in an axisparallel rectangle means that most pixels will be zero. Direct summation over such an enclosing rectangle is therefore inefficient.
An alternative representation of such a PSF is a list of tuples where each tuple contains coordinates and intensity of one support pixel of the PSF. By summation over this list, convolution is computed in time where is the number of support pixels. This procedure will be denoted as list filter. Unlike in the previous settings, we do not specifically consider PSFs with uniform intensity on all support pixels, since sparse blurs rarely satisfy such a condition.
4 Selective Convolution
In RL deconvolution each single iteration acts locally: It measures the error of the reblurred image compared to the observed image (by the quotient ), and redistributes this error by another blurring step (with the adjoint PSF ). Other iterative deconvolution methods that use convolution with the given PSF work in a similar way. Such a process can take many iterations to finally transport the error corrections to their correct locations, and achieve a good overall reconstruction.
Observation shows, however, that often the sharpness of the iterates improves dramatically during the first few iteration steps, and large parts of the image (typically, those with simpler structures) are well reconstructed already at this point. It is only a minority of the pixels in more complex structured image regions that cause the demand for many iterations. Based on this observation, it appears attractive to focus the computation to those pixels which really need the work, and to spare those pixels which are not or almost not changed during later iterations.
We implemented this idea in a very simple form: In each iteration, the absolute changes of all pixel values are checked, and those pixels whose change is below a given threshold are marked inactive in subsequent iterations, and thereby excluded from the expensive convolutions. For RL, this comes down to
if is active,  (4)  
if is inactive.  (5) 
Every ten iterations, all pixels are set back to active, thus one full iteration step is carried out, allowing to bring pixels back into the update process which still need small updates under the influence of more active pixels in their neighbourhood.
While more complicated adaptive update rules can be conceived, and are worth further research, an advantage of the present very simple rule is that the selection process itself requires almost no computational effort, which will also be demonstrated experimentally in the next section.
Selective convolution cannot be combined with those fast convolution techniques from Section 3 that rely on sliding windows or cumulated sums. However, it can be combined with the list filter.
5 Experimental Evaluation
All algorithms were implemented in C, and compiled with gcc 4.6.3 at optimisation level O2. Computations were run singlethreaded on a PC with an AMD Phenom X4 quadcore 64bit CPU clocked at 3.00 GHz under Ubuntu Linux 12.04. As this is not a realtime environment, interference of other processes in the system causes stochastic variations in runtime. Therefore each runtime was averaged from a series of 100 program runs, which reduced the standard deviation to
seconds or less in each series.Specific convolution operators.
In our first experiment, Table 1, we measure the computation time of 100 RL iterations with the different convolution implementations from Section 3. For reference, naive spatial convolution (direct summation over an axisparallel rectangle enclosing the PSF support) and an implementation via Fast Fourier Transform are included. The moderate PSF dimensions of and pixels correspond to common practical use cases.
PSF type  

PSF size  
Naive spatial  0.48  0.72  3.74  10.90  3.74  10.90  3.74  10.90 
Fourier  1.74  1.75  1.75  1.75  1.75  1.74  1.74  1.75 
List  0.50  0.78  3.14  10.73  0.51  0.80  1.97  7.37 
Generic box  0.17  0.17  0.93  1.56  0.91  1.48  0.93  1.54 
Box 2D sliding  0.13  0.13  0.70  1.19  —  —  —  — 
Box 2D cumul.  0.21  0.22  0.21  0.23  —  —  —  — 
Box 1D  0.09  0.09  —  —  —  —  —  — 
For each experiment, a test image was generated by convolving the cameraman image, see Fig. 3(a), with one of the following eight point spread functions: 1D box representing a linear motion blur in direction, 9 or 17 pixels long; 2D squareshaped box kernel of or pixels; diagonal () line kernel with or pixels, as a simple representative of a sparse point spread function; defocus kernel of or pixels diameter. Defocus blurring with diameter is shown in Fig. 3(b). As the convolution was carried out via the Fourier domain, periodic boundary conditions led to slight wraparound artifacts. Each test image was RLdeconvolved with its matching PSF, testing all those convolution routines which could handle the respective PSF.
As expected, in all cases the most specific applicable convolution algorithms lead to the fastest deconvolution. The generic box algorithm performs favourable for all “massive” convex 2D PSFs of constant density, while sparse PSFs are well treated by the list filter. Both implementations outperform FFTbased methods in their respective application areas. However, for the 2D box and defocus PSF of edge length/diameter the FFT method is close to break even with the generic box implementation.
Threshold  Omitted  SNR (dB)  time  speedup  

pixels (%)  vs. orig.  vs. ref.  (s)  
Reference  0.00  13.31  —  3.74  1.00 
0.00  13.31  —  3.77  0.99  
12.54  13.30  45.33  3.34  1.12  
22.32  13.31  42.32  2.99  1.25  
34.58  13.30  39.54  2.55  1.47  
53.69  13.25  34.74  1.84  2.03  
66.34  13.11  30.82  1.36  2.75  
75.57  12.84  27.18  1.01  3.70  
83.48  12.23  22.67  0.71  5.27 
a  b  c 

d  e 
Selective convolution.
In our second experiment, Table 2 and Fig. 3, we performed RL deconvolution with the naive spatial convolution but the convolution was carried out only for part of the pixels, as described in Section 4. The cameraman image convolved with the defocus PSF of diameter served as test case. For this test case 100 iterations of RL lead to a visible sharpening.
In Table 2, we report for different values of the thinning threshold the ratios of omitted pixels in convolutions, two signaltonoise ratios, run times (again averaged from 100 runs each) and speedup factors. The SNR value w.r.t. the original cameraman image allows to assess the reconstruction quality for each threshold, while SNR w.r.t. the reference image (100 iterations of unthinned RL) measures deviations introduced by the thinning.
The first line of the table contains the reference values for the RL method without thinning. Using the implementation with thinning but with the threshold set to zero (next line) reveals that the cost of the additional logic for thinning is not more than about . Raising the threshold up to speeds up the computation up to about 2.75 times while the reconstruction quality is almost not affected, as can be seen from the SNR figures, and is visually confirmed in Fig. 3(d). With larger thresholds and thus more aggressive thinning, deconvolution results are visibly affected, see Fig. 3(e).
6 Conclusions
We have demonstrated that the computational expense of iterative deconvolution methods can be substantially reduced if the point spread function is of one of several specific types that occur often in practice. The way to achieve this is to use specialised convolution operators. For small blurs, computation can be several times as fast as with Fourierbased convolution. Moreover, we have seen that some speedup (in our example about 2.5) can also be achieved by selectively performing convolution operations in later iteration steps only for image pixels where still significant change takes place.
From the deconvolution point of view, our experimental setting is prone to overassess reconstruction quality. In the context of the present paper, which is not at all about reconstruction quality but only concerned with algorithmic optimisations, this is not an issue, however.
A deeper investigation of the selective convolution approach including refined thinning rules, a broader experimental evaluation, and combination with the list filter is the subject of ongoing work. GPUbased parallelisation of the specialised convolution filters is another interesting topic for further research, as well as the integration of the techniques in blind or semiblind deconvolution frameworks.
References

[1]
L. Bar, N. Sochen, and N. Kiryati.
Image deblurring in the presence of saltandpepper noise.
In R. Kimmel, N. Sochen, and J. Weickert, editors,
Scale Space and PDE Methods in Computer Vision
, volume 3459 of Lecture Notes in Computer Science, pages 107–118. Springer, Berlin, 2005.  [2] D. S. C. Biggs and M. Andrews. Acceleration of iterative image restoration algorithms. Applied Optics, 36(8):1766–1775, Mar 1997.
 [3] T. F. Chan and C. K. Wong. Total variation blind deconvolution. IEEE Transactions on Image Processing, 7:370–375, 1998.
 [4] N. Dey, L. BlancFéraud, C. Zimmer, Z. Kam, J.C. OlivoMarin, and J. Zerubia. A deconvolution method for confocal microscopy with total variation regularization. In Proc. IEEE International Symposium on Biomedical Imaging (ISBI), April 2004.
 [5] V. P. Gerdt. On an algorithmic optimization in computation of involutive bases. Programming and Computer Software, 28(2):62–65, 2002.
 [6] T. J. Holmes and Y.H. Liu. Acceleration of maximumlikelihood image restoration for fluorescence microscopy and other noncoherent imagery. Journal of the Optical Society of America A, 8(6):893–907, 1991.
 [7] D. Krishnan and R. Fergus. Fast image deconvolution using hyperlaplacian priors. In Advances in Neural Information Processing Systems, pages 1033–1041, 2009.
 [8] L. B. Lucy. An iterative technique for the rectification of observed distributions. The Astronomical Journal, 79(6):745–754, June 1974.
 [9] M. J. McDonnell. Boxfiltering techniques. Computer Graphics and Image Processing, 17(1):65–70, 1981.
 [10] D. M. Mount, T. Kanungo, N. S. Netanyahu, C. Piatko, R. Silverman, and A. Y. Wu. Approximating large convolutions in digital images. IEEE Transactions on Image Processing, 10(12):1826–1835, 2001.
 [11] W. H. Richardson. Bayesianbased iterative method of image restoration. Journal of the Optical Society of America, 62(1):55–59, 1972.
 [12] R. Sedgewick. Algorithms. AddisonWesley, Reading, Massachusetts, 1984.
 [13] C. R. Vogel and M. E. Oman. Fast, robust total variationbased reconstruction of noisy, blurred images. IEEE Transactions on Image Processing, 7:813–824, 1998.
 [14] Y. Wang, J. Yang, W. Yin, and Y. Zhang. A new alternating minimization algorithm for total variation image reconstruction. SIAM Journal on Imaging Sciences, 1(3):248–272, 2008.
 [15] M. Welk. Robust variational approaches to positivityconstrained image deconvolution. Technical Report 261, Department of Mathematics, Saarland University, Saarbrücken, Germany, March 2010.
 [16] M. Welk, P. Raudaschl, T. Schwarzbauer, M. Erler, and M. Läuter. Fast and robust linear motion deblurring. Technical Report arXiv:1212.2245 [cs.CV], 2012.

[17]
N. Wiener.
Extrapolation, Interpolation and Smoothing of Stationary Time Series with Engineering Applications
. MIT Press, Cambridge, MA, 1949.  [18] Y.L. You and M. Kaveh. Anisotropic blind image restoration. In Proc. 1996 IEEE International Conference on Image Processing, volume 2, pages 461–464, Lausanne, Switzerland, September 1996.
Comments
There are no comments yet.