Noisy image block matching based on dissimilarity measure in discrete cosine transform domain

1. Introduction

The problem of image template location in a considered image is very common in different applications such as video compression, image clustering, vector quantization, nonlocal noise reduction [1–5]. To match templates, various measures of similarity/dissimilarity can be applied [1]. Among them, the most simple and popular measure is square L₂ norm or square Euclidean distance usually calculated as the sum of squared differences of corresponding pixels in two templates or image blocks. Despite its simplicity, the calculation of the square distances can cause significant computational load to find the best block matching in video compression [3] or modern noise suppression algorithm such as BM3D [5] and LPG-PCA [6], image inpainting [7]. The direct computation of all the matching requires O (KMNr²) operations [8], where M × N is the number of image pixels, r² denotes the patch size, K is the maximum offset.

In [3], an ample review of different existing block matching algorithms of reduced complexity are described: techniques based on a reduced set of candidates, techniques based on a reduced-complexity block distortion measure, techniques based on a subsampled block-motion field, hierarchical search techniques, fast full-search techniques. Among them, it was found that the fastest technique is the hierarchical search that usually performs on a mean pyramid calculated applying low-pass filtering and subsampling to the considered image(s) [3].

For many years, the patch matching based on the cross-correlation in Fourier transform domain [9, 12] was a state-of-the-art method [10]. This method requires approximately O (KMN  log ₂ (M)) operations [11] that is faster than the exhaustive search when r² > 30 log ₂ (M) [10]. However, when K is much larger than r, the direct method is faster than Fourier transform based method; the transform method is more efficient when r approaches K [11].

A fast method for image block dissimilarity computation of the sum of squared differences of pixel values within the patches is presented in [8]. It based on the integral image calculation [10] that permits to reduce the total time complexity independently from the patch size to O ((M - r) (N - r)) [8] operations at the expense of some flexibility and errors caused by limited numerical precision [8]. Besides, block matching using integral images has limited performance gain with respect to exhaustive block matching for small patch sizes [8].

From other hand, the matched blocks using the aforementioned methods become less similar when the image(s) are contaminated by noise. Obviously, the image data can be previously filtered to provide better matching but it can be inconvenient in some practical situations. For this reason, a similarity/dissimilarity measure robust to noise presence is required to find similar patches in noisy images.

In this paper, the feature vectors are formed as the patches in the discrete cosine transform domain that describe well the considered image blocks to perform fast block matching with the accuracy near to the full search technique or higher in the case of noise presence.

2. Patch dissimilarity measures

There exist a number of similarity/dissimilarity measures that can be applied to image processing. Here, only the simplest and fastest to calculate measures are described.

The standard dissimilarity measure between two image blocks of size r _x , r _y , sum of squared distances (SSD), can be formulated as d ( x , y , o x , o y ) = ∑ j = 0 r y - 1 ∑ i = 0 r x - 1 [ u ( x + i , y + j ) - u ( x + o x + i , y + o y + j ) ] 2 ,(1)

where u (x, y) is the image pixel at the coordinate x, y, o _x , o _y ∈ K denote the shift between the compared image blocks in the horizontal and vertical directions, respectively. The best matching is the minimum of the SSD found within the search window: u bestmatch ( x , y ) = arg min o x o y , ∈ K { d ( x , y , o x , o y ) } ,(2)

where the search area K = K _x · K _y is limited by offsets K _x , K _y in the horizontal and vertical directions, respectively. Often, the task is to find the set of the l best matched patches; this task is common in patch based filtering algorithms [5, 6]. Usually, it requires to calculate SSD (1) between all the blocks within the offset K, sort them and store the shortest block distances and their coordinates in the list p  (x, y) = [d₁, x₁, y₁; …; d _l , x _l , y _l ] ^T of the matched patches. The number of operations to find p may be large; for example, for the r² = 8 ×8 pixels block and the search window K = 20 × 20, at least O (Kr²) = 25600 sums and multiplications are required for each image pixel in (M - r) × (N - r) range. The simplest but not accurate way to speed up the calculations is find the distances between the means of the patches or even between the block central pixel values [7]. And the efficient block matching algorithm based on integral images requires O (4K) = 1600 operations for each pixel apart from the integral image calculation, but some errors caused by limited numerical precision might be introduced.

The number of operations to find p with a modification of a hierarchical search: firstly, the image is smoothed by some low-pass filter and then, distances are evaluated using only 4 pixels of the smoothed r² = 8 ×8 image blocks (without subsampling): d s ( x , y , o x , o y ) = [ u s ( x , y ) - u s ( x + o x , y + o y ) ] 2 + [ u s ( x + 4 , y ) - u s ( x + o x + 4 , y + o y ) ] 2 + [ u s ( x , y + 4 ) - u s ( x + o x , y + o y + 4 ) ] 2 + [ u s ( x + 4 , y + 4 ) - u s ( x + o x + 4 , y + o y + 4 ) ] 2(3)

where d _s denotes SSD calculated on the smoothed image u _s .

The image smoothing can be performed in spatial domain applying a low pass filter or in a transform domain. Such filtering apart from the complexity reduction introduces some robustness to the noise present in the contaminated images that improves the accuracy of block matching when searching for the minimum (2) or forming the list of patches p.

In this paper, for the low pass filtering is performed in discrete cosine transform (DCT) and consists in zeroing high frequency coefficients. Or similarly, one can take the first DCT coefficients that represent signal low frequencies and these coefficients can be efficiently used to calculate the dissimilarities between the image patches. Recall that 2D DCT for small r can be performed efficiently by matrix multiplications: u T ( x , y ) = A · u ( x , y ) · A T(4)

where u (x, y) is the image block at the coordinate x, y, u _T (x, y) is DCT transform of the image block at the coordinate x, y, and A is the DCT transformation matrix whose elements a_i,j are defined as [11] a i , j = β ( i ) 2 n cos ( π · i · ( 2 j + 1 ) 2 n ) , β ( i ) = { 1 2 , i = 1 1 otherwise(5)

The distances between the different patches in DCT domain may be calculated as d DCT ( x , y , o x , o y ) = ∑ ω y = 0 r ω ∑ ω x = 0 r ω [ u T ( x , y , ω x , ω y ) - u T ( x + o x , y + o y , ω x , ω y ) ] 2 ,(6)

where u _T (x, y, ω _x , ω _y ) is the ω _x , ω _y -th spectral coefficient of the image DCT patch u _T (x, y) at the coordinate x, y, o _x , o _y ∈ K denote the spatial shift between the compared image DCT patches in the horizontal and vertical directions, respectively, r _ω is the radius of the first spectral coefficients that participate in the distance calculation. Due to the energy compacting property of DCT, one can expect the distances calculated on the first DCT coefficients of the DCT patches, Equation (4), represent better the dissimilarities of the respective image blocks than d _h calculated on the smoothed image u _s . Obviously, when r = r _ω all spectral coefficient are participate in the distance calculation and the distance calculated with Equation (6) is equal to Equation (1)d _DCT (x, y, o _x , o _y ) = d (x, y, o _x , o _y ) because of energy normalization in Equation (5).

A simple modification of the distance Equation (6) is d absDCT ( x , y , o x , o y ) = ∑ ω y = 0 r ω ∑ ω x = 0 r ω [ | u T ( x , y , ω x , ω y ) | - | u T ( x + o x , y + o y , ω x , ω y ) | ] 2(7)

The distance Equation (7) may be useful, for example, to find the image clusters of homogeneous regions and edges, and is expected to be useful for deriving the spectral/statistical properties of the considered image blocks to improve the noise suppression of the modern noise suppression algorithms, such as BM3D or LG-PCA [5, 6].

In this paper, the distance Equation (6) is used to form the list p of matched patches when few DCT coefficients are involved in the distance calculation, r _ω = 2.7.

As shown below in Section 4, the properties of the measure Equation (6) are sometimes insufficient to determine with more precision the list p of matched patches, especially when noise is high and/or the considered image has texture areas; for these reasons, it must be improved.

To this end, the hierarchical approach may be useful: firstly, the algorithm searches for a rough approximation of the final list; secondly, a more exact search is performed over the provisional list found at the previous step.

In this Section the results of quantitative evaluation of the performance of different block matching algorithms described in Section 3 are presented.

To evaluate the results of dissimilarity calculation applied to block matching problem, the measure of similarity between two lists of the matched blocks p₁ and p₂ is introduced: the lists are more similar when more coordinates in p₁ are included in p₂ independently from the their positions in p₁ and p₂: C (k, m) = card (p _k ∩ p _m ). When both lists contain the same coordinates of the image patches, it means they are equal for non-local image filtering; as a result, the restored image fragments will not differ. And with the more different patches in the lists the image filtering result will be more different.

Additional measure of similarity/dissimilarity between two lists may be the distance between the patches with the same coordinates; and if the image restoration algorithm takes into account these distances for weighting, it affects the filtering results. Here, it is assumed that the non-local filtering algorithm does not use any weighting depending on the distances between the patches and does not take into account the values of these distances because the weighting may vary depending on the filtering algorithm.

Taking into account these observations, a measure for direct evaluation of the performance of different block matching algorithms is proposed. The algorithm ability to find similar block can be expressed in terms of a mean number of the coincided coordinates in the lists of the matched blocks p (x, y) and p _ref   (x, y) for all image pixels, normalized on the cardinality of the list p _ref and the total number of lists of the matched blocks: μ = ∑ y = 1 N - r ∑ x = 1 M - r card ( p ( x , y ) ∩ p ref ( x , y ) ) card ( p ref ) ( M - r ) ( N - r ) .(18)

where p _ref (x, y) denotes the reference list of patches, p (x, y) is the list of patches to compare.

Figures 1 and 2 illustrate the typical behavior of the accuracy of the block matching based on the distance Equation (6) in dependence of the radius r _ω and the noise variance. These results were obtained for the test image “Lenna” and “Baboon” for image block size r² = 8 ×8. As one can observe from these figures, the accuracy is decreased when the variance is growing for all radius lengths. However, the accuracy can be increased with the shorter r _ω for high noise levels. This increasing depends on the image data: for an image with plain areas, such as “Lenna” image, the results can be improved selecting r _ω = 3 for the noise variance σ² > 100. Contrary, for textural images, such as “Baboon”, the improvement is observed only for r _ω = 5, r _ω = 6. And only when the variance is greater 225, the results can be improved selecting r _ω = 4. This is due to the fact that the textural images contain high spatial frequencies, and suppressing these features with short r _ω results in a reduced accuracy of the block matching. In other words, the choice of the radius r _ω length depends on the local image block signal/noise ratio at the spatial frequencies that might be suppressed with a short r _ω .

OPEN IN VIEWER

Fig. 1.

Behavior of the accuracy Equation (18) of the simple block matching based on the distance Equation (6) depending on the radius r _ω and the noise variance for the test image “Lenna”, image block size is 8×8, search area K = 21×21.

OPEN IN VIEWER

Fig. 2.

Behavior of the accuracy Equation (18) of the simple block matching based on the distance Equation (6) depending on the radius r _ω and the noise variance for the test image “Baboon”, image block size is 8×8, search area K = 21×21.

A special case is when the radius r _ω = r and the size of the image block is equal to the size of DCT patch. In such conditions, the accuracy of the block matching based on the distance Equation (6) calculation is equal to the accuracy obtained with the Euclidean distance Equation (1), i.e., the performance is identical to the performance of the full search block matching, μ = 1.

To evaluate quantitatively and compare the accuracy of the of the proposed known image block matching algorithm to the known algorithms in noisy environment, their performance according to criterion Equation (18) was tested on the standard 512 × 512 test images “Lena”, “F-16”, “Peppers”, “Aerial”, “Baboon” and “Bridge” in the absence and presence of white additive Gaussian noise of different variance. The simulation results were obtained for image block size r² = 8 ×8 and the search area K = 21 × 21. As a reference, the lists of the matched blocks calculated by the full search algorithm on the noise free test images p _ref = p _FS were used. The obtained results are presented in the Table 1. Analyzing these results, one can conclude:

The proposed algorithms DCTMB, WDCTBM, ADCTBMandWADCTBM outperform the known algorithms (full search, integral image, hierarchical search) in the presence of noise in data;

New algorithms for similar image block matching in noisy environment have presented. The algorithms are based on the dissimilarity measure calculated as the distance between image patches in the discrete cosine transform domain. The algorithms perform the hierarchical search for the similar image blocks and hereby have a reduced complexity in comparison to the full search algorithm. Adjusting the radius of the distance calculation for spectral coefficient matching, the characteristics of the block matching algorithm can easily be adjusted to obtain a better accuracy of the matched block group. A higher accuracy is obtained using the local adaptation of the radius for the distance calculation outperforming the existing algorithms used to find groups of similar blocks. The performance of the different block matching algorithms were evaluated on the base of the proposed accuracy measure that uses as a reference the list of patches obtained with the full search algorithm in the absence of noise. The proposed locally adaptive algorithms ADCTBM, WADCTBM may be useful in different applications, such as image noise filtering, clustering and vector quantization. The evaluation of the proposed algorithm performance in such applications may be the subject of a future work.