Space-Variant Image Restoration with Running Sinusoidal Transforms

Various restoration methods (linear, nonlinear, iterative, noniterative, deterministic, stochastic, etc.) optimized with respect to different criteria have been introduced (Bertero & Boccacci, 1998; Biemond et al., 1990; Kundur, & Hatzinakos, 1996; Banham & Katsaggelos, 1997; Jain, 1989; Bovik, 2005; Gonzalez & Woods, 2008). These techniques may be broadly divided in two classes: (i) fundamental algorithms and (ii) specialized algorithms. One of the most popular fundamental techniques is a linear minimum mean square error (LMMSE) method. It finds the linear estimate of the ideal image for which the mean square error between the estimate and the ideal image is minimum. The linear operator acting on the observed image to determine the estimate is obtained on the basis of a priori second-order statistical information about the image and noise processes. In the case of stationary processes and space-invariant blurs, the LMMSE estimator takes the form of the Wiener filter (Jain, 1989). The Kalman filter determines the causal LMMSE estimate recursively. Specialized algorithms can be viewed as extensions of the fundamental algorithms to specific restoration problems. It is based on a state-space representation of the imaging system, and image data are used to define the state vectors. Specialized algorithms can be viewed as extensions of the fundamental algorithms to specific restoration problems. In this paper we deal with restoration of images degraded by space-variant blurs. Basically, all fundamental algorithms apply to the restoration of images degraded by space-variant blurs. However, because Fourier transforms cannot be utilized when the blur is space-variant, space-domain implementations of these algorithms may be computationally formidable due to large matrix operations. Several specialized methods were developed to attack the spacevariant restoration problem. The first class referred to as sectioning is based on assumption that the blur is approximately space-invariant within local regions of the image. Therefore, the entire image can be restored by applying well-known space-invariant techniques to the local image regions. A drawback of sectioning methods is the generation of artifacts at the region boundaries. The second class is based on a coordinate transformation (Sawchuk, 1974), which is applied to the observed image so that the blur in the transformed coordinates becomes space-invariant. Therefore, the transformed image can be restored by a space-invariant filter and then transformed back to obtain the final restored image. However, the statistical properties of the image and noise processes are affected by the


Introduction
Various restoration methods (linear, nonlinear, iterative, noniterative, deterministic, stochastic, etc.) optimized with respect to different criteria have been introduced (Bertero & Boccacci, 1998;Biemond et al., 1990;Kundur, & Hatzinakos, 1996;Banham & Katsaggelos, 1997;Jain, 1989;Bovik, 2005;Gonzalez & Woods, 2008).These techniques may be broadly divided in two classes: (i) fundamental algorithms and (ii) specialized algorithms.One of the most popular fundamental techniques is a linear minimum mean square error (LMMSE) method.It finds the linear estimate of the ideal image for which the mean square error between the estimate and the ideal image is minimum.The linear operator acting on the observed image to determine the estimate is obtained on the basis of a priori second-order statistical information about the image and noise processes.In the case of stationary processes and space-invariant blurs, the LMMSE estimator takes the form of the Wiener filter (Jain, 1989).The Kalman filter determines the causal LMMSE estimate recursively.Specialized algorithms can be viewed as extensions of the fundamental algorithms to specific restoration problems.It is based on a state-space representation of the imaging system, and image data are used to define the state vectors.Specialized algorithms can be viewed as extensions of the fundamental algorithms to specific restoration problems.In this paper we deal with restoration of images degraded by space-variant blurs.Basically, all fundamental algorithms apply to the restoration of images degraded by space-variant blurs.However, because Fourier transforms cannot be utilized when the blur is space-variant, space-domain implementations of these algorithms may be computationally formidable due to large matrix operations.Several specialized methods were developed to attack the spacevariant restoration problem.The first class referred to as sectioning is based on assumption that the blur is approximately space-invariant within local regions of the image.Therefore, the entire image can be restored by applying well-known space-invariant techniques to the local image regions.A drawback of sectioning methods is the generation of artifacts at the region boundaries.The second class is based on a coordinate transformation (Sawchuk, 1974), which is applied to the observed image so that the blur in the transformed coordinates becomes space-invariant.Therefore, the transformed image can be restored by a space-invariant filter and then transformed back to obtain the final restored image.However, the statistical properties of the image and noise processes are affected by the coordinate transformation.In particular, the stationarity in the original spatial coordinates is not preserved in the transform coordinate system.
In this chapter, we carry out the space-variant restoration using running discrete sinusoidal transform coefficients.The running transform is based on the concept of short-time signal processing (Oppenheim & Shafer 1989).A short-time orthogonal transform of a signal x k is defined as where w n is a window sequence, (n,s) represents the basis functions of an orthogonal transform.We use one-dimensional notation for simplicity.Equation ( 1) can be interpreted as the orthogonal transform of x k+n as viewed through the window w n .k s X displays the orthogonal transform characteristics of the signal around time k.Note that while increased window length and resolution are typically beneficial in the spectral analysis of stationary data, for time-varying data it is preferable to keep the window length sufficiently short so that the signal is approximately stationary over the window duration.Assume that the window has finite length around n=0, and it is unity for all n-N 1 , N 2 .Here N 1 and N 2 are integer values.This leads to signal processing in a running window (Vitkus & Yaroslavsky, 1987;Yaroslavsky & Eden, 1996).In other words, local filters in the domain of an orthogonal transform at each position of a moving window modify the orthogonal transform coefficients of a signal to obtain only an estimate of the pixel x k of the window.The choice of orthogonal transform for running signal processing depends on many factors.
We carry out the space-variant restoration using running discrete transform coefficients.The discrete cosine transforms (DCT) and discrete sine transforms (DST) are widely used.This is because the DCT and DST perform close to the optimum Karhunen-Loeve transform (KLT) for the first-order Markov stationary data (Jain, 1989).For signals with the correlation coefficient near to unity, the DCT provides a better approximation of the KLT than the DST.On the other hand, the DST is closer to the KLT, when the correlation coefficient lies in the interval (-0.5, 0.5).Since the KLT is constructed from the eigenvectors of the covariance matrix of data, there are neither single unique transform for all random processes nor fast algorithms.Unlike the KLT, the DCT and DST are not data dependent, and many fast algorithms were proposed.To provide image processing in real time, fast recursive algorithms for computing the running sinusoidal transforms are utilized (Kober, 2004(Kober, , 2007)).We introduce local adaptive restoration of nonuniform degraded images using several running sinusoidal transforms.Computer simulation results using a real image are provided and compared with those of common restoration techniques.

Fast algorithms of running discrete sinusoidal transforms
The discrete cosine and sine transforms are widely used in signal processing applications.Recently, forward and inverse algorithms for fast computing of various DCTs ({DCT-I, DCT-II, DCT-III, DCT-IV) and DSTs (DST-I, DST-II, DST-III, DST-IV) were proposed (Kober, 2004).

www.intechopen.com
Space-Variant Image Restoration with Running Sinusoidal Transforms 5

Discrete sinusoidal transforms
First, we recall the definitions for various discrete sinusoidal transforms.Notation {.} denotes a matrix, the order of which is represented by a subscript.For clarity, the normalization factor 2 N for all forward transforms is neglected until the inverse transforms.The kernel of the orthogonal DCT-I for the order N+1 is defined as 1 cos , where n, s=0,…, N; The kernels of the DCT-II, DCT-III, and DCT-IV for the order N are given as where n, s=0, 1,…, N-1.The kernel of the DST-I for the order N-1 is defined as where n, s=1, 2,…, N-1.The kernels of the DST-II, DST-III, and DST-IV for the order N are given as follows:   where n, s=1, 2,…, N.

www.intechopen.com
Image Restoration -Recent Advances and Applications 6

Fast forward algorithms for computing running discrete sinusoidal transforms
The fast forward algorithms are based on recursive relationships between three subsequent local running spectra.These algorithms for running DCTs (SDCTs) and running DSTs (SDSTs) are based on the second-order recursive equations summarized in Table I.

N Transforms
Recursive equations Table 1.Recursive equations for the computation of forward running sinusoidal transforms.

www.intechopen.com
Space-Variant Image Restoration with Running Sinusoidal Transforms 7 The number of arithmetic operations required for computing the running discrete cosine transforms at a given window position is evaluated as follows.The SDCT-I for the order N+1 with N=N 1 +N 2 requires N-1 multiplication operations and 4(N+2) addition operations.
The SDCT-II for the order N with N= N 1 + N 2 +1 requires 2(N-1) multiplication operations and 2N+5 addition operations.A fast algorithm for the SDCT-III for the order N with N=N 1 + N 2 +1 is based on the recursive equation given in line 3 of Table 1.Next it is useful to represent the equation as where the array   1 ;0 1 . 1 (here [x/y] is the integer quotient) and Eq. ( 10), the number of operations required to compute the DSCT-III can be evaluated as [3/2N] multiplication operations and 4N addition operations.An additional memory buffer of N elements is also required.Finally, the SDCT-IV for the order N with N=N 1 + N 2 +1 requires 3N multiplication operations and 3N+2 addition operations.
The number of arithmetic operations required for computing the running discrete sine transforms at a given window position can be evaluated as follows.The SDST-I for the order N-1 with N=N 1 + N 2 +1 requires 2(N-1) multiplication operations and 2N addition operations.
.N/2-1.Therefore, only N/2-1 multiplication operations are required to compute this term.The total number of multiplications is reduced to 3N/2-2.The SDST-II for the order N with N=N 1 + N 2 +1 requires 2(N-1) multiplication operations and 2N+5 addition operations.Taking into account the property of symmetry of the sine and cosine functions, the SDST-III for the order N with N=N 1 + N 2 +1 requires 2N multiplications and 4N addition operations.However, if N is even, the sum Finally, the SDST-IV for the order N with N=N 1 + N 2 +1 requires 3N multiplication operations and 3N+2 addition operations.The length of a moving window for the proposed algorithms may be an arbitrary integer.

Fast inverse algorithms for running signal processing with sinusoidal transforms
The inverse discrete cosine and sine transforms for signal processing in a running window are performed for computing only the pixel x k of the window.The running signal processing can be performed with the use of the SDCT and SDST algorithms.
The inverse algorithms for the running DCTs can be written as follows. IDCT-I: where N=N 1 +N 2 .If x k is the central pixel of the window; that is, N 1 =N 2 then the inverse transform is simplified to Therefore, in the computation only the spectral coefficients with even indices are involved.
The number of required operations of multiplication and addition becomes one and N 1 +1, respectively.

IDCT-II:
  where N=N 1 +N 2 +1.If x k is the central pixel of the window, that is, N 1 =N 2 then the inverse transform is given by We see that in the computation only the spectral coefficients with even indices are involved.
The computation requires one multiplication operation and N 1 +1 addition operations. IDCT-III: where If N 1 is even, then the computation requires N 1 +1 multiplication operations and 2N 1 addition operations.Otherwise, the complexity is reduced to N 1 multiplication operations and 2N 1 -1 addition operations. IDCT-IV: where N=N 1 +N 2 +1.If x k is the central pixel of the window, that is, N 1 =N 2 then the inverse transform is given by    One multiplication operation and N-1 addition operations are required to perform this computation.
The inverse algorithms for the running DSTs are given as follows. IDST-I: where N=N 1 +N 2 +2.If x k is the central pixel of the window; that is, N 1 =N 2 then the inverse transform is simplified to Therefore, in the computation only the spectral coefficients with odd indices are involved.The complexity is one multiplication operation and N 1 addition operations. IDST-II: where N=N 1 +N 2 +1.If x k is the central pixel of the window; that is, N 1 =N 2 then the inverse transform is given by   In the computation only the spectral coefficients with odd indices are involved.The computational complexity is one multiplication operation and N 1 +1 addition operations.

IDST-III:
   If N 1 is even, then the computation requires N 1 +1 multiplication operations and 2N 1 addition operations.Otherwise, the complexity is reduced to N 1 multiplication operations and 2N 1 -1 addition operations. IDST-IV: where =N 1 +N 2 +1.If x k is the central pixel of the window; that is, N 1 =N 2 , then the inverse transform is given by The complexity is one multiplication operation and N-1 addition operations.

Local image restoration with running transforms
First we define a local criterion of the performance of filters for image processing and then derive optimal local adaptive filters with respect to the criterion.One the most used criterion in signal processing is the minimum mean-square error (MSE).Since the processing is carried out in a moving window, then for each position of a moving window an estimate of the central element of the window is computed.Suppose that the signal to be processed is approximately stationary within the window.The signal may be distorted by sensor's noise.
Let us consider a generalized linear filtering of a fragment of the input one-dimensional signal (for instance for a fixed position of the moving window).Let a=[a k ] be undistorted real signal, x=[x k ] be observed signal, k=1,…, N, N be the size of the fragment, U be the matrix of the discrete sinusoidal transform, E{.} be the expected value, superscript T denotes the transpose.Let aH x  be a linear estimate of the undistorted signal, which minimizes the MSE averaged over the window The optimal filter for this problem is the Wiener filter (Jain, 1989): Let us consider the known model of a linear degradation: The obtained optimal filter is based on an assumption that an input signal within the window is stationary.The result of filtering is the restored window signal.This corresponds to signal processing in nonoverlapping fragments.The computational complexity of the processing is O(N 2 ).However, if the matrix of the optimal filter is diagonal, the complexity is reduced to O(N).Such filter is referred as a scalar filter.Actually, any linear filtering can be performed with a scalar filter using corresponding unitary transforms.Now suppose that the signal is processed in a moving window in the domain of a running discrete sinusoidal transform.For each position of the window an estimate of the central pixel should be computed.Using the equations for inverse sinusoidal transforms presented in the previous section, the point-wise MSE (PMSE) for reconstruction of the central element of the window can be written as follows: where

     
A Al HlXl    is a vector of signal estimate in the domain of a sinusoidal transform,   is a diagonal matrix of the size x NN of the coefficients of an inverse sinusoidal transform (see Eqs. ( 12), ( 14), ( 16), ( 18), ( 20), ( 22), (24), and ( 26)).Minimizing Eq. ( 32), we obtain where is sparse; the number of its non-zero entries is approximately twice less than the size of the window signal.Therefore, the computational complexity of the scalar filters in Eq. ( 33) and signal processing can be significantly reduced comparing to the complexity for the filter in Eq. ( 31).For the model of signal distortion in Eq. ( 30) the filter matrix is given as If the matrices   and UK W U TT aa in Eq. ( 34) are close to diagonal, the matrix of the scalar filter in ( 34) is close to diagonal, and the filter can be written as where , nn Pl Pl P l are diagonal elements of the following matrices For a real symmetric matrix of the covariance function, say K , there exists a unitary matrix U such that UKU T is a diagonal matrix.Actually, it is the KLT.It was shown (Jain, 1989) that some discrete sinusoidal transforms perform close to the optimum KLT for the firstorder Markov stationary data under specific conditions.In our case, the covariance matrices KWU TT aa and WK W T aa are not symmetric.Therefore, under different conditions of degradation different discrete sinusoidal transforms can better diagonalize these matrices.For instance, if a signal has a high correlation coefficient and a smoothed version of the signal is corrupted by additive, weakly-correlated noise, then the matrix For the design of local adaptive filters in the domain of running sinusoidal transforms the covariance matrices and power spectra of fragments of a signal are required.Since they are often unknown, in practice, these matrices can be estimated from observed signals (Yaroslavsky & Eden, 1996).

Computer simulation results
The objective of this section is to develop a technique for local adaptive restoration of images degraded by nonuniform motion blur.Assume that the blur is owing to horizontal relative motion between the camera and the image, and it is approximately space-invariant within local regions of the image.It is known that point spread functions for motion and focus blurs do have zeros in the frequency domain, and they can be uniquely identified by the location of these zero crossings (Biemond et al., 1990).We assume also that the observation noise is a zero-mean, white Gaussian process that is uncorrelated to the image signal.In this case, the noise field is completely characterized by its variance, which is commonly estimated by the sample variance computed over a low-contrast local region of the observed image.To guarantee statistically correct results, 30 statistical trials of each experiment for different realizations of the random noise process were performed.The MSE criterion is used for comparing the quality of restoration.Additionally, a subjective visual criterion is used.In our computer simulation, the MSE is given by where   () , 1 , . .ai i N  is the original image, and   is the restored image.The subjective visual criterion is defined as an enhanced difference between original and restored images.A pixel is displayed as gray if there is no error between the original image and the restored image.For maximum error, the pixel is displayed either black or white (with intensity values of 0 and 1, respectively).First, with the help of computer simulation we answer to the question: how to choose the best running discrete sinusoidal transform for local image restoration?

Choice of discrete sinusoidal transform for local image restoration
The objective of this section is to test the performance of the scalar filter (see Eq. ( 34)) designed with different running sinusoidal transforms for local image restoration while the statistics of the degraded image are varied.In our experiments we used realizations of a wide-sense colored stationary process, which is completely defined by the second-order statistics.The zero-mean process has the bi-exponential covariance function with varying correlation coefficient  .
The generated synthetic image is degraded by running 1D horizontal averaging of 5 pixels, and then a white Gaussian noise with a given standard deviation n  is added.In our tests the window length of 15x15 pixels is used, it is determined by the minimal size of details to be preserved after filtering.Since there exists difference in spectral distributions of the image signal and wide-band noise, the power spectrum of noise can be easily measured from the experimental covariance matrix.We carried out three parallel processing of the degraded image with the use of SDCT-II, SDST-I, and SDST-II transforms.At each position of the moving window the local correlation coefficient is estimated from the restored images.On the base of the correlation value and the standard deviation of noise, the resultant image is formed from the outputs obtained with either SDCT-II or SDST-I, or SDST-II according to Table 2.We also performed local image restoration using only the SDCT.As expected, the result of restoration is slightly worse than that of adaptive local restoration.We see that the proposed

Conclusion
In this chapter we treated the problem of local adaptive technique for space-variant restoring linearly degraded and noisy images.The minimum MSE estimator in the domain of running discrete sinusoidal transforms was derived.To provide image processing at high rate, fast recursive algorithms for computing the running sinusoidal transforms were utilized.Extensive testing using various parameters of degradations (nonuniform motion blurring and corruption by noise) has shown that the original image can be well restored by proper choice of the parameters of the proposed adaptive local restoration algorithm.
r e d i n a m e m o r y b u f f e r o f N elements.From the property of symmetry of the sine function, Fig. 1.(a) Covariance matrix of a noisy signal, (b) DCT-II of the covariance matrix.

154. 2
Local adaptive restoration of real degraded imageA real test aerial image is shown in Fig.2(a).The size of image is 512x512, each pixel has 256 levels of quantization.The signal range is [0, 1].The image quadrants are degraded by running 1D horizontal averaging with the following sizes of the moving window: 5, 6, 4, and 3 pixels (for quadrants from left to right, from top to bottom).The image is also corrupted by a zero-mean additive white Gaussian noise.The degraded image with the noise standard deviation of 0.05 is shown in Fig.2(b).
Fig. 2. (a) Test image, (b) space-variant degraded test image.The results of image restoration by the global parametric Wiener filtering (Jain, 1989) and the proposed method are shown in Figs.3(a) and 3(b), respectively.Figs.3(c) and 3(d) show differences between the original image and images restored by global Wiener algorithm and by the proposed algorithm, respectively.
Fig. 3. (a) Global Wiener restoration, (b) local adaptive restoration in domain of running transforms, (c) difference between the original image and restored image by Wiener filtering, (d) difference between the original image and restored image by proposed algorithm.

Fig. 4 .
Fig. 4. Performance of restoration algorithms in terms of MSE versus standard deviation of additive noise.
+1.If x k is the central pixel of the window; that is, N 1 =N 2 then we can rewrite 

Table 2 .
 and  ) is presented in Table2.So, in similar degradation circumstances and image model, local adaptive restoration yields the best results if the three sinusoidal transforms depending on local statistics of the processed image are used.The decision rule for choosing the best sinusoidal transform at each position of the moving window is given in Table2.Next, we carry out adaptive local restoration with real images.Best local restoration with running discrete sinusoidal transforms versus the model parameters.