BM3D算法半解，带python代码

1.思路

BM3D和NLM算法有一些相似之处，NLM的文章之前有写过：前往

他们都是用图像其他区域的特征与当前块特征融合成去噪后的图像块，主要的不同之处在以下几点：

        1.只在固定半径内搜索
        2.把有限个最相似的块叠成三维数组，先对图像块做水平方向上的2D的正交变换，再做纵向1D的正交变换（称为协同滤波（Collaborative Filtering））

第一点很好理解，就是减少了计算量，因为噪声在空间上是独立的，所以在单幅图片的情况下，仍然需要使用其它区域的相似图像来帮助当前区域进行去噪

而第二点则是最关键的思想：噪声在经过正交变换之后仍然存在，如经过DCT变换之后，图像真实信息往往只分布在较低的频率范围，而高频范围的大部分区域往往只有噪声，可以使用阈值将高频噪声置零的方法来粗略去噪，但是这样容易误伤图像的边缘、纹理等高频信息，尤其是在强噪声的情况下。

为了防止纹理被误伤，就需要进行纵向的正交变换，再次提取出相似块的相同频率位置的低频分量。（猜测）由于是相似块，所以同一个频率坐标上大概率有相似的真实值，这些值在纵向上的分布也是集中于低频的，再次正交变换之后，真实值更加集中，而噪声分量仍然遍布整个范围，使得去噪结果更加精确，能够更好地保留原有的高频分量

2.实际操作

BM3D算法其实是一个算法做了两次，但是第二次计算更加简化并引入维纳滤波。有时可能不必要做第二次计算。本文暂时讨论第一次的计算，第二次的还没看完。这里为了方便，正交变换全部采用DCT

---

算法步骤简述：

1.根据步长遍历图片，针对每个当前块，寻找指定半径内与其像素值的【DCT+欧式距离】最接近的一组相似块，组成块集合

2.纵向的1维DCT

3.使用固定阈值将幅度较低的频率系数置零，计算集合的零值总数

4.做两次反变换恢复整个图像块

5.将图像块融合为去噪后的图像块，与领块重叠的部分通过加权求和的方式，输出最终结果

---

第一步与NLM类似，搜索一定范围内的相似块，相似块的判断条件为：

其中是当前块，是参考块

第2，3，4步的方法比较简单，这里不详细说明

第5步主要是防止出现块效应而引入的，各块的权重由它对应的块集合中，非0值的数量决定()，具体公式为：

从相似的角度理解，如果图像集合完全同，那么经过3D变换，尤其是最后的1D变换之后，它的矩阵是非常稀疏的，因为纵向上每个坐标的值相同，对于纵向来说它就只有一个低频分量，因此非0值就非常少，权重更高，反之集合不相似，则权重更低

从噪声的角度理解，如果没有被置零的噪声分量越多，则3D变化后残留的系数越多，权重就越小

最后权重还要乘上kaiser窗，得到最后的还原公式

其中stack表示所有重叠块，z表示纵向坐标，k表示二维的kaiser窗

最后附上python实现：

# -*- 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 #图像搜索区域半径
blocksize = 4#图像块尺寸
blockstep = 2#搜索步长
blockmaxnum = 8#相似块最大数量
searchblocksize = searchsize/blocksize#半径内块个数
kai = np.kaiser(blocksize,5)
kai = np.dot(kai[:,None],kai[None,:])#二维kaiser

diffT = 100 #允许纳入数组的最大diff
coefT = 0.0005

diffArray = np.zeros(blockmaxnum)#相似块
similarBoxArray = np.zeros((blocksize,blocksize,blockmaxnum))
#jinbuqu csdn newland


#让图像尺寸符合倍数
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]

#初始化分子分母
imnum = np.zeros(imnoise.shape)
imden = np.zeros(imnoise.shape)

#将左上角作为块的坐标.每个块独立做一次去噪，只要关注一块的计算就行
for y in range(0,newh-blocksize,blockstep): #检查能不能完整走到底,范围对不对
    systart = max(0,y - searchsize)
    syend = min(newh - blocksize,y + searchsize - 1)
    process_bar(y/(newh-blocksize), start_str='进度', 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)
        
        #排序矩阵初始化
        similarBoxArray[:,:,0] = imnoise[y : y + blocksize,x : x + blocksize]
        hasboxnum = 1
        diffArray[0] = 0
        #不算自己
        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
                #塞入
                
                changeid = 0
                if hasboxnum < blockmaxnum - 1:
                    changeid = hasboxnum
                    hasboxnum = hasboxnum + 1
                else:
                    #排序
                    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
        

        #开始做dct
        for difid in range(1,hasboxnum):
            similarBoxArray[:,:,difid] = cv2.dct(similarBoxArray[:,:,difid])
        
        #开始做1d,阈值操作，计算非零个数
        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
            
        #最终恢复图像所用到的分子分母
        for difid in range(1,hasboxnum):
            #求权重：kaiser * w
            weight = kai /( (sigma ** 2) * notzeronum)
            imidct = cv2.idct(similarBoxArray[:,:,difid])
            
            #分子，分母
            imnum[y : y + blocksize,x : x + blocksize] += imidct * weight
            imden[y : y + blocksize,x : x + blocksize] +=  weight
            

        
        x = x + blocksize
    
    y = y + blocksize


#完成去噪
imout = imnum/imden



cv2.imshow('in',imnoise)
cv2.waitKey(0)
cv2.imshow('out',imout)
cv2.waitKey(0)

效果：