1.

Introduction

Nowadays, a great number of image filters intended for removal of noise of a given type exists (e.g., for additive,¹ Poisson,² speckle,³ and impulse⁴ noise). Other filters are capable to remove different types of mixed noise, such as additive and impulsive,^{5, 6} multiplicative and impulsive,⁷ multiplicative and additive,⁸ etc., under assumption that a mixed noise type and its basic parameters are known in advance or preestimated. In order to perform filtering in an efficient way, one has to have a preliminary knowledge on noise type and statistics.^{5, 9, 10} On the other hand, such a priori information may not always be available in many typical practical situations that often arise in remote sensing, hyperspectral imaging, ultrasound medical diagnostics, cDNA imaging, etc. ^{11, 12, 13, 14, 15}

For many imaging systems, noise is not a stationary process and its statistics and properties could be quite dissimilar in different fragments of a particular image. In practice, dissimilarities in noise statistics can be due to several basic reasons. One reason is the change of imaging conditions or considerably different distances to particular parts of an observed scene as in side-look radar imaging.¹⁶ Another reason could be uncontrollable influence of several different sources of noise. For instance, in maritime radars, interference arises due to sea clutter that depends on wind speed, direction of sea waves, incidence angle,¹⁶ etc. A third reason could be nonlinear amplification regulations in input circuits before image digitalization as in ultrasound medical devices. Finally, for some imaging systems, adequate models of noise present in the acquired or transferred images have not been established or commonly accepted yet. ^{17, 18, 19, 20}

Therefore, there are typical and atypical situations concerning properties of noise present in images. Typical situations considered thoroughly in image-filtering literature are characterized by the following: noise type is known or its characteristics are known a priori or can be estimated for a given type of noise. On the contrary, less-studied situations addressed in this paper are the following:

These shortcomings led to design locally adaptive filters (see Refs. 7, 18, 24, 25, 26 and references therein), including several approaches that employ partial differential equations, anisotropic diffusion, and total variation minimization.^{27, 28, 29} They have demonstrated considerable improvement of performance in comparison to nonadaptive nonlinear filters due to exploiting different mechanisms of adaptation to image local content. However, some of them (e.g., Refs. 18, 24, 25), require a priori information on noise type and statistics, other ones are rather complicated.

There are also many techniques that perform adaptive filtering in the domain of an orthogonal transform. To name a few, let us mention discrete wavelet transform and discrete cosine transform (DCT)–based filters for the additive noise case,^{1, 30, 31, 32} different filters designed for speckle removal,^{3, 33, 34} and other types of noise.^{2, 35, 36, 37} However, for most of these filters, it is supposed that noise type is known and its parameters are either known in advance or preestimated with appropriate accuracy by preliminary “global” analysis of an image at hand. ^{36, 38}

If noise statistics is unknown and noise is nonstationary, the following set of questions arises:

The paper is organized as follows. In Section 2, basic principles of DCT-based filtering are described. Main advantages of this approach are mentioned and demonstrated. The ways how to adapt this technique to signal-dependent noise are shown. Section 3 deals with description of two proposed modifications of locally adaptive DCT-based filter. Preliminary testing for simulated images and noise is carried out in Section 4, which also contains data for comparison of the designed methods performance to one of the state-of-the-art methods. Limitations of the method for spatially correlated noise are demonstrated. Real-life image-processing examples are given in Section 5. Finally, we draw conclusions and present possible directions of future work.

2.

Basic Principles of Adaptive DCT-Based Filtering

All methods of transform-based signal-image-filtering (denoising) rely on the same basic assumption that a signal in transform domain has more sparse representation than noise.^{30, 33, 35, 36, 37, 38, 39, 40} Then, if a component in orthogonal transform domain has a rather large absolute value, it most probably corresponds to information content and should be preserved. If a component has amplitude close to zero (less than a threshold), then it most probably relates to noise and has to be removed (suppressed). Then, filtering can be performed using different orthogonal transforms (wavelets, discrete cosine, Haar, etc.). A choice of basis functions as well as filtering efficiency depend on how high are energy compaction properties of a transform for a given signal (image).

The Karhunen–Loeve transform possesses the best signal decorrelation property in least-squares terms.^{41, 42} Unfortunately, its coefficients are data dependent and, thus, no fast algorithm of the Karhunen–Loeve transform exists. DCT is almost as efficient as the Karhunen-Loeve transform in terms of signal decorrelation and energy compaction, especially in the case of signals modeled by Markov processes.^{41, 42} There are also fast algorithms for calculating the DCT.^{42, 43} Because of these properties, DCT is widely used in image and video compression.^{44, 45}

Let us recall an operation principle of DCT-based denoising.^{31, 33, 37, 42} In general form, the two-dimensional DCT coefficients of an $N_1\times N_2$ matrix $\mathbf{A}$ may be defined as

After the thresholding, the set of changed coefficients { $D_T(k,l,n,m)$ : $k=1,2,\dots ,N$ , $l=1,2,\dots ,N$ } is obtained for each image block defined by the indices $n$ and $m$ . Next, the inverse DCT is applied to each block of thresholded coefficients and the preliminary filtered image values $I_{\mathrm{pf}}(q,s,n,m)$ are obtained for $q=n,n+1,\dots ,n+N-1$ , $s=m,m+1,\dots ,m+N-1$ . Because now almost all image pixels belong to $N\times N$ different image blocks, we have for each image pixel $(q_o,s_o)$ usually $N\times N$ preliminary filtered pixel values $I_{\mathrm{pf}}(q_o,s_o,n,m)$ , where $n=q_o,q_o-1,\dots ,q_o-N+1$ and $m=s_o,s_o-1,\dots ,s_o-N+1$ . Only the pixels near the image borders are exceptions because, for them, the number of such output estimates is smaller (minimally, one for all four image corners).

Then, at the last step, all output estimates for each pixel must be combined in order to obtain the final filtered image $I_{\mathrm{f}}$ . The simplest way to do this is to average these estimates. The resulting value for the pixel $(q_o,s_o)$ is then

A threshold $T(k,l,n,m)$ may vary, depending on spatial coordinates defined by indices $n$ and $m$ , or may be frequency dependent⁴⁶ in the case when noise is spatially correlated with known or preestimated spatial correlation characteristics. We assume that there is no a priori information about noise spatial correlation. Thus, we consider only threshold dependence on $n$ and $m$ .

Concerning the potential of DCT-based filtering with fully overlapping blocks, we would like to recall that it provides performance (in terms of output MSE or PSNR) better or comparable to the best wavelet-based denoising techniques for pure additive and pure multiplicative noise cases.^{31, 33} Furthermore, the main advantage of DCT-based filtering for the considered case of nonstationary noise is that, being carried out blockwise, it can be easily adapted to local statistical properties of noise.^{33, 47} To show how this can be done, let us consider the following two image models. The first one is a the case of film grain noise for which a general observation model is

Let us confirm this fact with two examples. Figure 1 presents the plots of output MSE on $\beta$ for the DCT filter described above ( $8\times 8\text{pixel}$ blocks with full overlapping). These plots have been obtained for two standard test images (Lenna and Baboon) corrupted by additive white Gaussian noise [ $\gamma =0$ in Eq. 3] with variances ${\sigma }_n^2=50$ and ${\sigma }_n^2=100$ , and by Poisson noise. In the former case, $T(n,m)$ is set as $\beta {\sigma }_n$ . In the latter case, $T(n,m)=\beta {\left[\widehat{I}(n,m)\right]}^{1∕2}$ .

Fig. 1

Output MSE versus $\beta$ for two standard test images: (a) Lenna and (b) Baboon, corrupted by additive white Gaussian noise with variances 50 (dotted line) and 100 (dashed line) and by Poisson noise (solid line).

As it is seen, all dependencies have minimums that are observed for $\beta$ of $\sim 2.6$ for the image Lenna and $\sim 2.3$ for the image Baboon. For other standard test images (Goldhill, Peppers, Barbara) minimums are observed for $\beta$ approximately equal to 2.6. Similar to other filtering algorithms, noise suppression is less efficient if a processed image has a more complex structure and/or the noise variance is smaller. Note that for additive noise with ${\sigma }_n^2=100$ , the values of output peak signal-to-noise ratio (PSNR) provided by the DCT filter with $\beta =2.6$ are equal to 35.40, 34.42, and $34.61\mathrm{dB}$ for the images Lenna $(512\times 512)$ , Barbara $(512\times 512)$ , and Peppers $(256\times 256)$ , respectively. This is better than for one of the best state-of-the-art Kervrann’s filter (see Ref. 26, Table IV, $\mathrm{PSNR}=35.18$ , 33.79, and $34.07\mathrm{dB}$ , respectively).

Thus, DCT-based filtering has a high potential. However, such performance of the DCT filter is observed under condition of known noise type and statistics. And we are interested in the case of non-stationary noise with unknown characteristics.

As is seen from Eq. 4, the product ${\sigma }_n\left[\widehat{I}(n,m)\right]$ is, in fact, an estimate of noise standard deviation (SD) $\widehat{\sigma }(n,m)$ in a given block. For $\gamma =0$ (i.e., for pure additive noise), one gets exactly ${\sigma }_n$ . The threshold is fixed and equal to $\beta {\sigma }_n$ . Similarly, for pure multiplicative case $\widehat{\sigma }(n,m)={\sigma }_{\mu }(n,m)\widehat{\overline{I}}(n,m)$ where ${\sigma }_{\mu }$ is the multiplicative noise normalized SD. As an estimate of local mean $\widehat{\overline{I}}(n,m)$ , one can use the DCT spectral coefficient $D(1,1,n,m)$ .

Thus, generalizing the presented results, the conclusion is as follows. To set a local threshold, one needs a local estimate of noise standard deviation $\widehat{\sigma }(n,m)$ . Consider the following image additive observation model

3.

Locally Adaptive DCT Filters

3.1.

DCT Filter Adaptive to Image Local Statistics

Here, we discuss the fourth question given in Section 1 (i.e., what methods of local estimation of noise statistics can be used). Some initial imagination concerning properties of local estimates of noise variance in blocks or scanning windows with fixed size can be got from earlier studies.^{49, 50, 51, 52} Generally speaking, noise local SD can be estimated in the spatial or spectral domain.^{49, 50, 51, 52, 53} These local estimates are then aggregated (processed) in a robust manner to provide blind estimation of noise characteristics based on a model of noise assumed known a priori. Pure additive or multiplicative noise variances can be blindly estimated.^{49, 50, 51, 52} More sophisticated methods allow estimating dependence of local variance on local mean for signal-dependent noise.^{36, 53}

The following are useful conclusions drawn from several papers: ^{36, 49, 50, 51, 52, 53}

• 1 Accuracy of local estimates of noise variance depends on several basic factors, namely, block size, image content, true values of noise variance, spatial correlation of noise, a used method for obtaining local estimates, and its parameters.

• 2 Better accuracy (in terms of smaller relative errors) is provided for simpler content (less textural) images if noise variance is rather large and noise is spatially uncorrelated.

• 3 The appropriate block (scanning window) size is from $7\times 7\text{to}9\times 9\text{pixels}$ ⁵ (i.e., $8\times 8\text{pixel}$ blocks used in DCT-based filtering conform well with this recommendation).

• 4 In any case, there are normal local estimates of noise variance obtained in homogeneous image blocks that are close to true values of noise variance and abnormal estimates obtained in heterogeneous image blocks that are considerably larger than the true values of noise variance.

These conclusions indicate that it is hard to expect good accuracy of noise variance evaluation for blocks of a limited size, especially if image content in a given block is quite complex and noise is not intensive. To prove this, consider first local estimates of noise variance in spatial domain

Eq. 6

$${\widehat{\sigma }}_c^2(n,m)=\sum _{q=n}^{n+7}\sum _{s=m}^{n+7}\frac{{(I_{qs}-{\widehat{\overline{I}}}_{nm})}^2}{63},$$

where $I_{qs}$ is $qs$ ’th pixel value within $nm$ ’th image block and ${\widehat{\overline{I}}}_{nm}$ denotes the $nm$ ’th block mean. It is easy to show that

Eq. 7

$${\widehat{\sigma }}_c^2(n,m)\approx \mathrm{Var}\left\{I^{\mathrm{tr}}(n,m)\right\}+{\sigma }_n^2(n,m),$$

where $\mathrm{Var}\left\{I^{\mathrm{tr}}(n,m)\right\}={\sum }_{q=n}^{n+7}{\sum }_{s=m}^{n+7}[{(I_{qs}^{\mathrm{tr}}-{\overline{I}}_{nm}^{\mathrm{tr}})}^2∕63]$ , ${\overline{I}}_{nm}^{\mathrm{tr}}={\sum }_{q=n}^{n+7}{\sum }_{s=m}^{n+7}{I}_{qs}^{\mathrm{tr}}∕64$ , ${\sigma }_n^2(n,m)$ denotes the noise variance in the $nm$ ’th block. Thus, conventional estimates of local variance [Eq. 6] are sensitive to image content in a block. It is an inherent property of such estimates, and the question is how to exploit this sensitivity in a reasonable manner. It can be a drawback or a positive feature for the considered situation of unknown noise variance. The reasons will be considered next.

Let us now concentrate on estimation of noise local SD in the spectral (DCT) domain. For a given block, it is possible to estimate $\widehat{\sigma }(n,m)$ in a spectral domain as^{39, 48}

Eq. 8

$$\widehat{\sigma }(n,m)=1.483\mathrm{med}\{\left|D(k,l,m,n)\right|:k=1,\dots ,8;l=1,\dots ,8;\text{except}k=l=1\},$$

where med{ } denotes the median of a data sample. This estimate is similar to the median of absolute deviations used as a robust data scale estimator.⁹ Here [in Eq. 8], the noise SD (scale) estimate $\widehat{\sigma }(n,m)$ exploits properties of transform coefficients when mixed Gaussian–Laplacian or Gaussian scale mixture models for their description are used.⁵⁴ Two properties of orthogonal transforms and robust estimation are exploited in Eq. 8. First, noise after an orthogonal transform spreads between all coefficients, and, if noise is independent identically distributed (i.i.d), its power spreads uniformly where spectral coefficients occur to be Gaussian random variables. Then, the SD of spectral coefficients is proportional to the noise SD. Second, the use of the robust estimate (median) provides less sensitivity to outliers where the outliers in the considered case are induced by image content.

However, although the used estimator (sample median) is robust, the image local content (details, texture, edges) present in a given block leads to the positive bias of the estimate $\widehat{\sigma }(n,m)$ .⁵² Note that if noise is spatially correlated, considerable bias can be observed for estimates obtained according to Eq. 6⁵¹ (these estimates in homogeneous regions are, on average, smaller than the true value of noise SD).

Let us analyze how great the influence of image local content and other factors could be on accuracy of local estimates. Consider a test image [e.g., the standard test image Barbara, (Fig. 2 )] corrupted by white additive noise with a zero mean and constant variance ${\sigma }_n^2=100$ . The values of $\widehat{\sigma }(n,m)$ from Eq. 6 are visualized in Fig. 2; larger values are lighter (magnification by 7 is used for better visualization). It is seen that the values $\widehat{\sigma }(n,m)$ for blocks that correspond to homogeneous image regions vary a little. They are mostly rather small and approximately equal to ${\sigma }_n$ of additive noise. The values $\widehat{\sigma }(n,m)$ for the blocks that correspond to image heterogeneities such as edges, textures, and details, are commonly slightly larger due to influence of image content.

Fig. 2

(a) Noise-free test image Barbara, (b) visualized $\widehat{\sigma }(n,m)$ of test image Barbara (values have been multiplied by 7) and (c) ${\widehat{\sigma }}_c(n,m)$ (values have been multiplied by 3.5).

Behavior of ${\widehat{\sigma }}_c(n,m)$ is similar [see Fig. 2 magnification is by 3.5], but its sensitivity to image heterogeneities in blocks is considerably larger. This can be also seen from the comparison of histograms of $\widehat{\sigma }(n,m)$ and ${\widehat{\sigma }}_c(n,m)$ presented in Figs. 3 and 3 respectively. In both cases, normal local estimates (which are obtained in homogeneous image blocks) are grouped near the true value equal to 10. But for the histogram in Fig. 3 there are more abnormal local estimates and they are, on average, larger.

Fig. 3

Histograms of (a) $\widehat{\sigma }(n,m)$ and (b) ${\widehat{\sigma }}_c(n,m)$ for the image Barbara, ${\sigma }_n^2=100$ .

This analysis allows proposing a locally adaptive DCT-based filter that will be further referred as LA DCT-1. This filter algorithm is as follows:

• 1 For each given block, estimate $\widehat{\sigma }(n,m)$ according to Eq. 8 and set a local threshold as $T(n,m)=\beta \widehat{\sigma }(n,m)$ .

• 2 All other operations (DCT in blocks, thresholding, inverse DCT, and averaging for overlapping blocks) are performed as described in Section 2.

We prefer to use $\widehat{\sigma }(n,m)$ , but not ${\widehat{\sigma }}_c(n,m)$ , because $\widehat{\sigma }(n,m)$ is less sensitive to heterogeneities. Note that the use of a larger threshold in DCT-based filtering leads to oversmoothing. As follows from the algorithm description, the proposed filter LA DCT-1 adapts to noise local characteristics. This is adaptation mechanism 1.

3.2.

Improved Locally Adaptive DCT Filter

In this section, we consider one more way (mechanism 2) to further improve performance of the locally adaptive DCT-based filter. Our approach is based on the following general idea and assumptions. Suppose that we have succeeded in discriminating homogeneous and heterogeneous regions of an image. Because in homogeneous regions, local estimates $\widehat{\sigma }(n,m)$ are quite close to true values of local SD of noise, then it is a correct decision to set the local threshold equal to $2.6\cdot \widehat{\sigma }(n,m)$ . On the contrary, if a given block corresponds to an image heterogeneous region, then, most probably, $\widehat{\sigma }(n,m)$ is larger than the local SD of noise. Then, it is reasonable to set the local threshold less than $2.6\cdot \widehat{\sigma }(n,m)$ as $\beta (n,m)\widehat{\sigma }(n,m)$ , where $\beta (n,m)<2.6$ . In general, there are many possible ways to set $\beta (n,m)$ . It can be fixed and equal to some ${\beta }_{\mathrm{het}}⩽2.6$ in a simplest case or it can be determined in a more complicated manner.

Thus, a primordial task is to design some discriminator for homogeneous and heterogeneous blocks. It is easy to resolve this task if one knows the local SD of noise in advance.⁷ However, it becomes more complicated if the local SD of nonstationary noise is unknown. To get around this shortcoming, let us exploit the properties of local estimates $\widehat{\sigma }(n,m)$ (in spectral domain) and ${\widehat{\sigma }}_c(n,m)$ (in spatial domain) established in Section 3.1.

Note that usually⁷ ${\widehat{\sigma }}_c(n,m)$ is defined in square-shaped scanning windows, where window side size is odd (e.g., 5 or 7).^{7, 23} Because here we deal with DCT-based filtering in $8\times 8$ blocks, it is possible to calculate ${\widehat{\sigma }}_c(n,m)$ either according to Eq. 6 or as follows:

Eq. 9

$${\widehat{\sigma }}_c^2(n,m)=\sum _{q=n+1}^{n+Q}\sum _{s=m+1}^{n+S}\frac{{(I_{qs}-{\widehat{\overline{I}}}_{qs})}^2}{QS-1},{\widehat{\overline{I}}}_{qs}=\sum _{q=n+1}^{n+Q}\sum _{s=m+1}^{n+S}\frac{I_{qs}}{QS},Q=S=6.$$

The map of the estimates ${\widehat{\sigma }}_c(n,m)$ for $Q=S=6$ is visualized in Fig. 2 for the test image Barbara corrupted by zero mean additive noise with the constant variance ${\sigma }^2=100$ . Because the estimate ${\widehat{\sigma }}_c(n,m)$ is much more sensitive to image heterogeneities in a block than the estimate $\widehat{\sigma }(n,m)$ , let us exploit this difference for discriminating homogeneous and heterogeneous image blocks. For this purpose, consider a ratio $R(n,m)=[{\widehat{\sigma }}_c(n,m)∕\sigma (n,m)]$ . Its example for the test image Barbara and additive Gaussian noise with ${\sigma }_n^2=100$ is visualized in Fig. 4 . We represent the ratio map as an image $R_{\mathrm{vis}}(n,m)=\left[40R(n,m)\right]$ , $n=1,\dots ,N_{\mathrm{im}}-7$ , $m=1,\dots ,M_{\mathrm{im}}-7$ , where [•] means rounding to the nearest positive integer not larger than 255, $N_{\mathrm{im}}\times M_{\mathrm{im}}$ is the image size, 40 is a magnification coefficient used to visualize better the maps $R(n,m)$ .

Fig. 4

Magnified ratio image $R_{\mathrm{vis}}(n,m)=\left[40R(n,m)\right]$ for the test image Barbara (a) corrupted by additive noise and (b) corrupted by Poisson noise.

Joint visual analysis of the test image Barbara [Fig. 2] and the ratio map [Fig. 4] shows that quite large ratios (essentially larger than unity, indicated by brighter color pixels) are observed for blocks located in the heterogeneous image regions. Figure 4 presents the ratio map for the same test image but corrupted by Poisson noise. The obtained maps are rather similar. This indicates that the method of analyzing local activity (heterogeneity) based on $R(n,m)$ is applicable to different types of noise. This observation allows expectation that the proposed principle will also work well enough for any nonstationary noise.

There are several possible ways to exploit the aforementioned property. The simplest one is to adaptively set the local threshold

Eq. 10

$$T(n,m)=\{\begin{array}{cc}\beta \widehat{\sigma }(n,m),& \text{if}R(n,m)<T_R\\ {\beta }_{\mathrm{het}}\widehat{\sigma }(n,m),& \text{if}R(n,m)⩾T_R\end{array}\phantom{\}},$$

where $\beta =2.6$ as in LA DCT-1, $T_R$ is a preset threshold, ${\beta }_{\mathrm{het}}<2.6$ is a factor that determines the hard threshold for DCT coefficients in edge/detail neighborhoods and in textural regions. The DCT-based filter described by Eqs. 5, 9, 10 (further referred as to LA DCT-2) belongs to the class of locally adaptive hard-switching filters,⁷ where hard switching relates to the parameter $\beta$ . This filter implies both adaptions to noise characteristics (mechanism 1) and to image content (mechanism 2).

LA DCT-2 uses two new parameters, $T_R$ and ${\beta }_{\mathrm{het}}$ . Consider distributions of the ratio $R(n,m)$ , which is a random variable. Suppose that the estimates ${\widehat{\sigma }}_c(n,m)$ and $\widehat{\sigma }(n,m)$ are obtained in a homogeneous image region. Then, their means under condition of i.i.d. Gaussian noise are the same. Thus, it is possible to expect that the distribution of the ratio $R(n,m)$ should have a mode in the neighborhood of unity. As follows from the ratio maps in Figs. 4 and 4, ${\widehat{\sigma }}_c(n,m)$ is usually larger than $\widehat{\sigma }(n,m)$ in image heterogeneous regions. Therefore, the distribution of the ratio $R(n,m)$ may have a heavy “right-hand” tail. The histogram of the obtained $R(n,m)$ values for the test image Barbara (pure additive Gaussian noise, ${\sigma }_n^2=100$ ) is presented in Fig. 5 . The histogram for the image Lenna corrupted by Poisson noise is shown in Fig. 5. Both distributions, as expected, have modes in the neighborhoods of unity, and they possess a heavy right-hand tail. Similar shapes of $R(n,m)$ value distributions have been observed for other conventional test images as Peppers, Goldhill, etc. The detailed analysis of values of these histograms shows that the distribution mode, which corresponds to homogeneous image regions, is about unity. Being random, the values $R(n,m)$ form the quasi-Gaussian part of this distribution. As in any discrimination or detection task, a larger preset discrimination threshold provides better correct detection of homogeneous regions but larger probability of recognizing heterogeneous region as homogeneous. Thus, threshold setting is a compromise. Analysis of histogram data shows that a trade-off in discriminating blocks that correspond to homogeneous and heterogeneous regions can be provided by setting $T_R\approx 1.35$ .

Fig. 5

Histograms of the ratio $R(n,m)$ for the test images: (a) Barbara corrupted by zero mean pure additive Gaussian noise with ${\sigma }_n^2=100$ and (b) Lenna corrupted by Poisson noise.

4.

Local Adaptive DCT Filter Performance Analysis

4.2.

Filter Performance in Spatially Correlated Noise Environment

The case of spatially uncorrelated (i.i.d.) noise has been studied. However, in practice, it is possible that noise can be spatially correlated. Furthermore, 2-D autocorrelation function of spatially correlated noise can be unknown in advance.^{7, 51}

We have analyzed how the proposed LA DCT-2 filter performs in the situation of spatially correlated noise. For this purpose, we have simulated additive Gaussian zero mean spatially correlated noise with variance 100. Spatially correlated noise has been modeled by applying the $3\times 3$ window mean filter to originally i.i.d. zero mean Gaussian 2-D noise and adjusting a desired variance. The results are presented in Table 4 for three values of $T_R$ and several values of ${\beta }_{\mathrm{het}}$ . The test images are the same as earlier.

Table 4

Values of MSEout for several βhet and TR for the considered test images corrupted by spatially correlated noise.

Image		Lenna	Barbara	Baboon	Peppers	Goldhill
TR	βhet
1.37	2.6	54.2	74.7	189.5	74.0	82.3
	2.3	76.8	76.5	178.5	75.6	83.3
	2.0	79.1	78.7	167.8	77.0	84.2
	1.7	81.4	81.3	157.6	78.2	85.2
	1.4	83.8	84.1	147.7	79.3	86.2
	1.1	86.1	87.1	138.7	80.2	87.1
	0.8	88.1	89.7	131.4	81.0	88.1
1.33	2.3	76.9	76.6	177.3	75.8	83.3
	2.0	79.4	78.6	165.7	77.7	84.4
	1.7	81.9	81.6	154.6	78.8	85.5
	1.4	84.4	84.6	143.9	80.1	86.6
	1.1	86.9	87.6	134.3	81.1	87.8
	0.8	89.1	90.5	126.6	82.2	88.9
1.29	2.3	77.1	76.6	176.3	75.9	83.4
	2.0	79.7	79.1	163.6	77.7	84.6
	1.7	82.3	81.9	151.6	79.4	85.8
	1.4	85.1	85.0	140.2	80.8	87.1
	1.1	87.7	88.2	130.0	82.2	88.4
	0.8	90.1	91.2	122.0	83.4	89.7

Analysis of the data presented in Table 4 shows the following:

• 1 For rather simple images (Lenna, Peppers, Barbara, Goldhill), there is some reduction of noise, but it is considerably smaller than in the earlier considered case of i.i.d. noise. For the image Baboon, it is still difficult to provide appropriate performance of the filter.

• 2 For the case of spatially correlated noise, the recommendation to set $T_R\approx 1.3$ , ${\beta }_{\mathrm{het}}\approx 1.5$ is incorrect. The results obtained for $T_R=1.37$ and ${\beta }_{\mathrm{het}}=2.6$ occurred to be better than for $T_R\approx 1.3$ , ${\beta }_{\mathrm{het}}\approx 1.5$ for four test images.

Thus, additional analysis of why it happens is needed. Some reasons are clear; they come from studies done recently.^{44, 56} In particular, it has been demonstrated that frequency-dependent thresholds should be used in DCT-based filter if noise is spatially correlated. Denoising with frequency-dependent thresholding can be carried out if either spatial correlation properties of noise are known in advance or preestimated. Thus, it becomes necessary to adapt not only to local statistics of noise and image local content, but also to spatial correlation characteristics of noise.

Another reason is that statistical characteristics of ${\widehat{\sigma }}_c(n,m)$ and $\widehat{\sigma }(n,m)$ as well as their ratio $R(n,m)$ change if noise is spatially correlated with respect to the i.i.d. noise case. Let us analyze histograms of $R(n,m)$ for spatially correlated noise. Two typical histograms are presented in Fig. 6 . Their shapes are similar to the shapes of the histograms represented in Fig. 5. However, there is an obvious difference. The modes of the histograms in Fig. 6 occur not in the neighborhood of unity but for larger values.

Fig. 6

Histograms of the ratio $R(n,m)$ for the test images: (a) Lenna and (b) Barbara corrupted by zero mean spatially correlated additive Gaussian noise with ${\sigma }_n^2=100$ .

This phenomenon has its explanation. For spatially correlated noise, both estimates ${\widehat{\sigma }}_c(n,m)$ and $\widehat{\sigma }(n,m)$ in homogeneous blocks (that form histogram maximum) become biased, smaller than the true value of noise SD.⁵¹ But the bias ${\Delta }_c$ of the estimates ${\widehat{\sigma }}_c(n,m)$ can be about $(0.1\dots 0.2){\sigma }_n^2$ , whereas for the estimates $\widehat{\sigma }(n,m)$ the bias $\Delta$ can reach $(0.7\dots 0.8){\sigma }_n^2$ .⁵¹ Thus, the estimates ${\widehat{\sigma }}_c(n,m)$ in image homogeneous blocks are, on average, larger than the estimates $\widehat{\sigma }(n,m)$ . This leads to mode shifting to values larger than unity.

Let us analyze the mode estimates $\widehat{r}$ for distributions of $R(n,m)$ . The obtained data are presented in Table 5 . The method⁴⁹ has been applied for mode estimation.

Table 5

Estimates r̂ for the test images corrupted by different types of noise.

Noise type	Test image
	Lenna	Barbara	Baboon	peppers	Goldhill
Additive white	1.006	1.034	1.141	1.000	1.047
Poisson	1.005	1.034	1.122	0.997	1.049
Additive spatially correlated	1.646	1.779	1.432	1.434	1.602

As is seen, for spatially uncorrelated noise the estimates $\widehat{r}$ are quite close to unity for all test images for both additive i.i.d. Gaussian and Poissonian noise. One interesting observation is that, for images with more complicated structure (e.g., Baboon), the estimates $\widehat{r}$ are larger. Similarly, if noise is spatially correlated, the estimates $\widehat{r}$ are considerably larger than 1.0 for all five considered test images.

Two conclusions follow from this analysis:

• 1 Practical cases when an image is either corrupted by spatially correlated noise or is textural can be discriminated from cases when an image is rather simple and corrupted by i.i.d. noise. This can be done by comparing $\widehat{r}$ to some threshold $T_{\mathrm{spc}}$ (e.g., approximately equal to 1.15: if $\widehat{r}<T_{\mathrm{spc}}$ ), image filtering by LA DCT-2 can be efficient (see data in Section 4.1), whereas in the opposite case, a more complicated filtering should be applied.

• 2 If noise is spatially correlated, it seems reasonable to use $T_R>1.3$ ; adaptation to noise spatial correlation is needed.

Thus, we have determined practical situations for which the designed filter LA DCT-2 is worth applying. A decision of whether or not filtering has to be applied can be undertaken automatically by analyzing $\widehat{r}$ for an image at hand. Design of adaptive DCT-based filters for complex-structure images corrupted by spatially correlated noise is a topic of future research.

5.

Performance Analysis for Real-Life Images

Let us present two examples of applying the proposed filters to real-life data. One example is a set of polarimetric radar data, presented as real valued (floating point) data arrays. These data have been obtained by a maritime coastal radar. The images are presented in Fig. 7 . HH means that a horizontally polarized signal is transmitted and then received; VH relates to a vertically polarized emitted and horizontally polarized received signal.

Fig. 7

Data registered for (a) HH and (b) VH polarization of radar signal.

A sensed area mainly corresponds to sea surface (that occupies basic part of images and form background), one large rock (small rocky island placed in the image central part, very bright pixels), several small ones (left lower corner, groups of bright pixels), and a shadowed zone behind the large rock. The horizontal axis of images corresponds to range and vertical axis relates to azimuth of the radar. During data acquisition, internal gain control was used to partly compensate dependence of backscattered signal mean on distance.

Visual analysis of these images allows concluding that background intensity varies depending on range. More detailed analysis has been done to confirm this conclusion. First, histograms (Fig. 8 ) have been obtained for manually selected six “homogeneous” regions (marked by frames in images in Fig. 7) for two different sectors and three different mean distances. Analysis of sample histograms in Fig. 8 shows that the histogram shape changes with range. The distributions are asymmetric with respect to their means. There is a heavy tail to the right side from the distribution mode. Distribution modes are also different for different ranges. Outliers are seldom (occur with quite small probabilities), their values differ from the mode values by approximately ten times.

Fig. 8

Sample histograms for the image in Fig. 7 (a) for small distance fragment 1 and (b) for middle distance fragment 3.

To carry out more thorough analysis, we have determined the following parameters for each image homogeneous region:

The results obtained for the image in Fig. 7 are given in Table 6 . For more distant fragments (5 and 6), the values $m_i$ and ${\sigma }_i^2$ are larger than for less distant fragments, especially for the fragments 1 and 2. In turn, the standard and robust (interquantile) estimates of relative variance for more distant fragments are smaller. This shows that noise is neither pure additive nor pure multiplicative. A specific “distance-dependent” noise is observed. This conclusion has been confirmed by results of analysis carried out for the image in Fig. 7. The observed nonstationarity results from joint influence of several factors, namely, specific features of the receiver amplifier gain control used, varying incidence (grazing) angle of backscattering from wavy sea surface, different mutual geometry of sea wave direction and radar azimuth, etc. This is the case when it is difficult to separate the influence of these factors. We assume that nonstationary noise local statistics does not change abruptly [i.e., it is possible to consider it practically constant for fragments (blocks) of relatively small size, let us say $8\times 8$ blocks commonly used in DCT-based filtering].

Table 6

Statistical characteristics of the selected “homogeneous” regions for the image in Fig. 7.

Figure 9 presents the processed images for the maritime radar data carlier given in Fig. 8. Note that in this case, a two-stage procedure has been applied. At the first stage, the CW Median filter²³ with the $3\times 3$ scanning window and the center element weight $W_c=5$ has been used to remove outliers. Then, at the second stage, the LA DCT-1 has been used. Comparing images in Figs. 7 and 9, one can conclude that noise is well suppressed and useful information in images is preserved well enough, although some oversmoothing of fine details is observed.

Fig. 9

Output images for (a) HH and (b) VH polarization of radar signal.

Let us give another example of a real-life image processing. The 224th subband image of the Lunar Lake AVIRIS^{57, 58, 59} image is presented in Fig. 10 . Noise present in this image is visually seen. The output image for LA DCT-1 is shown in Fig. 10. Noise is well suppressed, although sharp edges are slightly smeared. The ratio image $R_{\mathrm{vis}}(n,m)=\left[40R(n,m)\right]$ is visualized in Fig. 10. The most sharp edges and details are marked by brighter color pixels in this map. This allows better preservation of sharp edges and details in the output image by the designed LA DCT-2 [see Fig. 10 $T_R=1.35$ , ${\beta }_{\mathrm{het}}=1.7$ , $\widehat{r}=1.03$ ].

Fig. 10

(a) Original one-channel image of AVIRIS hyperspectral data, (b) the output of LA DCT-1 filter, (c) the ratio $R_{\mathrm{vis}}(n,m)$ image, and (d) the output of the hard switching LA DCT-2 filter.

6.

Conclusions and Future Work

It is shown that there are practical situations where noise is nonstationary and limited a priori information on its statistics is available. To filter images effectively under these conditions, two DCT-based filtering techniques to suppress nonstationary spatially uncorrelated noise have been proposed and studied.

The first technique is a locally adaptive filter based on local estimation of noise variance in blocks and setting the corresponding threshold proportionally to the obtained estimate of noise SD. It performs well enough for rather simple images. The second locally adaptive filter employs, in addition, the analysis of the ratio $R(n,m)$ and adapts to image content in a block. This leads to decreasing ${\mathrm{MSE}}_{\mathrm{out}}$ and improving edge/detail/texture preservation in processed images. As a result, the performance of this filter is comparable to performance of the best state-of-the-art methods.

The recommendations concerning proper selection of the filter adaptation parameters are given. The designed filters have been applied to real-life images and have demonstrated excellent results. The spatially correlated noise case has been studied as well. The proposed methods do not perform well enough for this case and should be further modified. A way to discriminate the cases of spatially uncorrelated and correlated noise has been proposed. The design of locally adaptive filters for spatially correlated noise case may be a subject of our future work.