IBD Algorithm

The iterative blind deconvolution algorithm proposed by Ayers and Dainty is the most popular method in its class. The basic structure of the algorithm is illustrated in Fig. 2.

Figure 2: IBD Algorithm

First a nonnegative-valued initial estimate $f_0(x,y)$ is input into the iterative scheme [1]. This function is Fourier transformed to give $\hat{F}(u,v)$. From this we form the second function's spectrum, $\tilde{H}(u,v)$. Then we impose image-domain constraint of non-negativity to the inverse transform of this function, by simply putting to zero all points of the function that have a negative value. A positive constrained image $\hat{h}(x,y)$ is consequently formed that is Fourier transformed to give the spectrum $\hat{H}(u,v)$. This is inverted and multiplied by $G(u,v)$ to give the next spectrum estimate $\tilde{F}(u,v)$. Taking inverse Fourier transform of this function and imposing non-negativity constraints complete the loop of a single iteration. The iterative loop is repeated until two positive functions with required convolution $g(x,y)$ is found.

This basic approach has two major problems to deal with:

1) Defining the inverse filter in regions where the function to be inverted has low values, is difficult.

2) We do not have any information at the spatial frequencies where $G(u,v)$ or $F (u,v)$ are zero.

To attack these problems, we change the way we implement the Fourier constraints.At each iteration we have two estimates for $F (u,v)$. $\hat{F}(u,v)$ has a non-negative inverse transform and $\tilde{F}(u,v)$ satisfies the Fourier constraints. So at each iteration those estimates are averaged to form a new estimate Eq.(2). The weight parameter $\beta$ is important for the convergent rate. The regions below the noise level in the convolution, are dealt with by only using the estimate $\hat{F}(u,v)$ Eq.(1).

If the modulus of $\hat{H}(u,v)$ is less than the modulus of $G(u,v)$, then the inverses of the two function estimates are averaged Eq. (3). This takes care of the small values of $G(u,v)$, which are likely to be noisy. All these can be summarized as follows:

if $\vert G(u,v)\vert <$ noise level,

$$\tilde{F}_{i+1}(u,v) = \hat{F}_{i}(u,v)$$ (1)

if $\vert\hat{H}_i(u,v)\vert \geq \vert G(u,v)\vert$ ,

$$\tilde{F}_{i+1}(u,v) = (1-\beta) \hat{F}_{i}(u,v) + \beta \frac{G(u,v)}{\hat{H}_i(u,v)}$$ (2)

if $\vert\hat{H}_i(u,v)\vert < \vert G(u,v)\vert$ ,

$$\frac{1}{\tilde{F}_{i+1}(u,v)} = \frac{(1-\beta)}{ \hat{F}_{i}(u,v)} + \beta \frac{\hat{H}_i(u,v)}{G(u,v)}$$ (3)

Matlab Source Code : ibd.m