### 1. Ideas

There are some similarities between BM3D and NLM algorithms. NLM has written before: go to

They all use the features of other regions of the image and the current block features to fuse into denoised image blocks. The main differences are as follows:

1. Search only within a fixed radius

2. Stack a limited number of the most similar blocks into a three-dimensional array. First do the 2D orthogonal transformation in the horizontal direction for the image block, and then do the 1D orthogonal transformation in the vertical direction (called Collaborative Filtering)

The first point is to reduce the amount of calculation. Because the noise is spatially independent, it is still necessary to use similar images of other regions to help denoise the current region in the case of a single image

The second point is the key idea: noise still exists after orthogonal transformation. For example, after DCT transformation, the real information of the image is often only distributed in the lower frequency range, while most areas of the high frequency range are often only noise. The method of setting the high frequency noise to zero by threshold can be used to roughly denoise, but it is easy to damage the edge of the image Texture and other high-frequency information, especially in the case of strong noise.

In order to prevent the texture from being hurt by mistake, it is necessary to carry out longitudinal orthogonal transformation to extract the low-frequency components of the same frequency position of similar blocks again. (guess) because it is a similar block, there is a high probability that there are similar real values on the same frequency coordinate. The vertical distribution of these values is also concentrated in the low frequency. After the orthogonal transformation again, the real values are more concentrated, and the noise components still spread throughout the whole range, making the denoising results more accurate and better retaining the original high-frequency components

### 2. Actual operation

BM3D algorithm is actually an algorithm, which has been done twice, but the second calculation is more simplified and Wiener filtering is introduced. Sometimes it may not be necessary to do a second calculation. This article temporarily discusses the first calculation, and the second one is not finished. Here, for convenience, all orthogonal transforms adopt DCT

Description of algorithm steps:

1. Traverse the picture according to the step size, and for each current block, find a group of similar blocks closest to its pixel value [DCT + Euclidean distance] within the specified radius to form a block set

2. Longitudinal 1D DCT

3. Use a fixed threshold to set the frequency coefficient with lower amplitude to zero and calculate the total number of zero values of the set

4. Do two inverse transformations to restore the whole image block

5. The image block is fused into the denoised image block, and the overlapping part with the leading block outputs the final result by weighted summation

The first step is similar to NLM. Search for similar blocks within a certain range. The judgment conditions of similar blocks are as follows:

amongIs the current block,Is a reference block

The methods of steps 2, 3 and 4 are relatively simple and will not be described in detail here

Step 5 is mainly introduced to prevent block effect. The weight of each block is determined by the number of non-zero values in its corresponding block set(), the specific formula is:

From the perspective of similarity, if the image sets are exactly the same, then after 3D transformation, especially the last 1D transformation, its matrix is very sparse. Because the value of each coordinate is the same in the longitudinal direction, it has only one low-frequency component in the longitudinal direction, so the non-zero value is very small and the weight is higher. On the contrary, if the sets are not similar, the weight is lower

From the perspective of noise, if the more noise components are not set to zero, the more coefficients remain after 3D change and the smaller the weight

Finally, the weight should be multiplied by kaiser window to get the final reduction formula

Where stack represents all overlapping blocks, z represents longitudinal coordinates, and k represents two-dimensional kaiser window

Finally, the python implementation is attached:

# -*- coding: utf-8 -*- import cv2 import numpy as np import math def process_bar(percent, start_str='', end_str='', total_length=0): bar = ''.join(["\033[31m%s\033[0m"%' '] * int(percent * total_length)) + '' bar = '\r' + start_str + bar.ljust(total_length) + ' {:0>4.1f}%|'.format(percent*100) + end_str print(bar, end='', flush=True) image = cv2.imread('F:/tf_learn--------/impixiv/593.jpg') image = cv2.cvtColor(image,cv2.COLOR_BGR2GRAY) size = image.shape # h w c image = np.array(image/255, dtype=np.float32) noise = np.random.normal(0, 0.001 ** 0.5, size) imnoise = image + 1*noise #cv2.imshow('result',imnoise) #cv2.waitKey() sigma = 1 searchsize = 16 #Image search area radius blocksize = 4#Image block size blockstep = 2#Search step blockmaxnum = 8#Maximum number of similar blocks searchblocksize = searchsize/blocksize#Number of blocks in radius kai = np.kaiser(blocksize,5) kai = np.dot(kai[:,None],kai[None,:])#Two dimensional kaiser diffT = 100 #Maximum diff allowed to be included in the array coefT = 0.0005 diffArray = np.zeros(blockmaxnum)#Similar block similarBoxArray = np.zeros((blocksize,blocksize,blockmaxnum)) #jinbuqu csdn newland #Make the image size conform to multiple newh = math.floor((size[0] - blocksize)/blockstep)*blockstep - searchsize neww = math.floor((size[1] - blocksize)/blockstep)*blockstep - searchsize newh = math.floor((size[0] - newh - blocksize)/blockstep)*blockstep + newh neww = math.floor((size[1] - neww - blocksize)/blockstep)*blockstep + neww imnoise = imnoise[0:newh,0:neww] #Initialize numerator and denominator imnum = np.zeros(imnoise.shape) imden = np.zeros(imnoise.shape) #Use the upper left corner as the coordinates of the block Each block is denoised independently, as long as we pay attention to the calculation of one block for y in range(0,newh-blocksize,blockstep): #Check whether it can go to the end completely. Is the scope right systart = max(0,y - searchsize) syend = min(newh - blocksize,y + searchsize - 1) process_bar(y/(newh-blocksize), start_str='speed of progress', end_str="100", total_length=15) for x in range(0,neww-blocksize,blockstep): sxstart = max(0,x - searchsize) sxend = min(neww - blocksize,x + searchsize - 1) #Sorting matrix initialization similarBoxArray[:,:,0] = imnoise[y : y + blocksize,x : x + blocksize] hasboxnum = 1 diffArray[0] = 0 #Not yourself for sy in range(systart,syend,blockstep): for sx in range(sxstart,sxend,blockstep): if sy == y and sx == x: continue diff = np.sum(np.abs(imnoise[y : y + blocksize,x : x + blocksize] - imnoise[sy : sy + blocksize,sx : sx + blocksize])) if diff > diffT: continue #Insert changeid = 0 if hasboxnum < blockmaxnum - 1: changeid = hasboxnum hasboxnum = hasboxnum + 1 else: #sort for difid in range(1,blockmaxnum - 1): if diff < diffArray[difid]: changeid = difid if changeid !=0: similarBoxArray[:,:,changeid] = imnoise[sy : sy + blocksize,sx : sx + blocksize] diffArray[changeid] = diff #Start doing dct for difid in range(1,hasboxnum): similarBoxArray[:,:,difid] = cv2.dct(similarBoxArray[:,:,difid]) #Start 1d, threshold operation and calculate the non-zero number notzeronum = 0 for y1d in range(0,blocksize): for x1d in range(0,blocksize): temp3ddct = cv2.dct(similarBoxArray[y1d,x1d,:]) zeroidx = np.abs(temp3ddct)<coefT temp3ddct[zeroidx] = 0 notzeronum = notzeronum +temp3ddct[temp3ddct!=0].size similarBoxArray[y1d,x1d,:] = cv2.idct(temp3ddct)[:,0] #jinbuqu csdn newland if notzeronum < 1: notzeronum = 1 #The numerator and denominator used to finally restore the image for difid in range(1,hasboxnum): #Weight calculation: kaiser * w weight = kai /( (sigma ** 2) * notzeronum) imidct = cv2.idct(similarBoxArray[:,:,difid]) #Numerator imnum[y : y + blocksize,x : x + blocksize] += imidct * weight imden[y : y + blocksize,x : x + blocksize] += weight x = x + blocksize y = y + blocksize #Complete denoising imout = imnum/imden cv2.imshow('in',imnoise) cv2.waitKey(0) cv2.imshow('out',imout) cv2.waitKey(0)

effect: