Python Implementation for BM3D Denoising of Color Image

BM3D denoising is a popular method used for removing image noise. The original BM3D paper [1] proposed the algorithm for gray scale image. Then in [2], the method is extended for color images.

Python script for grayscale image BM3D can be found here by liuhuang31. But Python script for color image BM3D can't be identified after searching in Internet. Note that this Python implementation supports BM3D for color images but the kernel part is written in C by Marc Lebrun. By modifying liuhuang31's script, we provide a Python implementation for color BM3D here. In this implementation, all codes including algorithm kernel part is written in Python.

Below is an example including noisy image and images after step 1 and step 2 of denoising processing of BM3D. Noisy image is generated by adding random noise to reference noise-free image. PSNR is peak SNR calculated between each individual image and the reference noise-free image.

Noisy image (PSNR = 22.41dB)

After step 1 of removing noise (PSNR = 29.36dB)

After step 2 of removing noise (PSNR = 30.13dB)

While our Python script is an extension of liuhuang31's work, we should note that we have made a few corrections:

Correction 1:
In the original code, delta of two images is calculated as below. But since img1 and img2 are both uint8 type, the calculation result is wrong.

D = numpy.array(img1 - img2, dtype=numpy.int64)

In the modified code, img1 and img2 are first typecast to int64 and then delta between them is found.

D = numpy.array(img2, dtype=numpy.int64) - numpy.array(img1, dtype=numpy.int64)

Correction 2:
In the original code, to find boundary, it uses the following code. But last line has a typo. shape[0] should be shape[1].

    if LX < 0:   LX = 0
    elif RX > _noisyImg.shape[0]:   LX = _noisyImg.shape[0]-_WindowSize
    if LY < 0:   LY = 0
    elif RY > _noisyImg.shape[0]:   LY = _noisyImg.shape[0]-_WindowSize

The modified code is:

    if LX < 0:   LX = 0
    elif RX > _noisyImg.shape[0]:   LX = _noisyImg.shape[0]-_WindowSize
    if LY < 0:   LY = 0
    elif RY > _noisyImg.shape[1]:   LY = _noisyImg.shape[1]-_WindowSize

Correction 3:
The last one is probably the trickiest one. To calculate Wiener filter in step 2, noise variance is sigma^2. To match with noise variance, signal power should be normalized by the count of similar blocks. This normalization is not in liuhuang31's original code. Other BM3D code such as Marc Lebrun's code also uses this normalization.

                Norm_2 = numpy.float64(tem_Vct_Trans.T * tem_Vct_Trans)
                m_weight = Norm_2/Count/(Norm_2/Count + sigma_color[ch]**2)

[1] K. Dabov et. al., Image denoising by sparse 3D transform-domain collaborative filtering, IEEE Trans. Image Proc., Vol 16-8, Aug 2007
[2] K. Dabov et. al, Color image denoising via sparse 3D collaborative filtering with grouping constraint in luminance-chrominance space, ICIP 2007