noise:时域降噪，图像平均的效果

VIDEO DENOISING BASED ON ADAPTIVE TEMPORAL AVERAGING
1.这篇论文时域去噪的论文很简单，就是对图像序列的当前帧查找附近的相似帧，然后进行平均。对于位移比较大的物体或者相机位移比较大时，都会出现重影的效果。

def denoise_adaptive_temp_average(imgs, thrA, thrB):
    imgs2 = []
    # 首先计算帧间差异
    N = len(imgs)
    print(N)
    d = np.zeros([N, N], dtype=np.float32)
    for i in range(N):
        for j in np.arange(i, N, 1):
            d[i][j] = np.mean(np.abs(imgs[i] - imgs[j]))
            d[j][i] = d[i][j]
    print("cal d finish! ")
    for i in np.arange(0, N):
        l = i
        r = i
        dl = 0
        dr = 0
        suml = 0
        sumr = 0
        while(dl < thrA and suml < thrB):
            if l > 0:
                l -= 1
                dl = d[l][i]
                suml += dl
            else:
                break
            #print(dl, suml, thrA, thrB)
        while (dr < thrA and sumr < thrB ):
            if r < N-1:
                r += 1
                dr = d[i][r]
                sumr += dr
            else:
                break
            #print(dr, sumr, thrA, thrB)
        # 计算 left - right image 的均值
        im_t = np.zeros(imgs[0].shape, dtype=np.float32)
        #print(im_t.dtype, imgs[0].dtype)
        for k in np.arange(l+1, r, 1):
            im_t += imgs[k]

        cnt = (r - l - 1)
        if cnt == 0:
            cnt = 1
            im_t = imgs[i]
        im_t = im_t // cnt
        imgs2.append(im_t.astype(np.uint8))

        #print(i, l, r, suml, sumr)
    return imgs2

2.即使只利用左右两帧的效果，也会有重影。可见判断运动时非常重要的。当r=1时，只利用左右两帧进行平均的code

def denoise_temp_average(imgs, r):
    imgs2 = []
    # 首先计算帧间差异
    N = len(imgs)
    s = imgs[0].shape

    for ii in np.arange(0, r, 1):
        imgs2.append(imgs[ii])

    for ii in np.arange(r, N-r):
        # cur = imgs[i]
        sum = np.zeros(s, dtype=np.float32)
        for k in np.arange(-r, r+1):
            sum += imgs[ii+k]

        sum = sum // (2*r+1)
        sum = np.clip(sum, 0, 255).astype(np.uint8)
        imgs2.append(sum)
    for ii in np.arange(N - r, N, 1):
        imgs2.append(imgs[ii])
    return imgs2

3.如果对每个像素进行判断呢，和当前像素在一定阈值内的相邻帧同一位置像素，作为候选点，然后进行平均。理论上应该更好，实际确实是感观更好一些

def denoise_adaptive_temp_pixel_average(imgs, r=3, thr1=10, thr2=0.2):
    '''
    前后各r帧中， 相同位置的pixel, 如果值相差在一定范围内，则作为候选帧
    '''
    s = imgs[0].shape
    imgs2 = []
    # 首先计算帧间差异
    N = len(imgs)
    if 2*r+1 > N:
        print("error: frame number is not enough, ", 2*r+1, "  ", N)
        return [-1]
    print("num is enough !")
    for ii in np.arange(0, r, 1):
        imgs2.append(imgs[ii])

    for i in np.arange(r, N - r):

        masks = np.zeros(imgs[0].shape, dtype=np.int32)
        sum = np.zeros(imgs[0].shape, dtype=np.float32)

        for k in np.arange(-r, r+1):
            mask1 = np.abs((imgs[i] - imgs[i+k])) < thr1
            mask2 = np.abs((imgs[i] - imgs[i+k]) / (imgs[i] + imgs[i+k] + 0.0001)) < thr2

            mask = np.logical_and(mask1, mask2).astype(np.int32)


            sum += imgs[i+k]*mask
            masks += mask
            #print(i, k, mask)
            # plt.figure()
            # plt.subplot(121)
            # plt.imshow((mask*255).astype(np.uint8), cmap='gray')
            # plt.subplot(122)
            # plt.imshow(imgs[i+k]*mask)

            #print(np.sum(masks<1)) # should be 0
            #print(masks)
        sum = sum // masks
        imgs2.append(np.clip(sum, 0, 255).astype(np.uint8))

    for ii in np.arange(N-r, N, 1):
        imgs2.append(imgs[ii])

    #plt.show()
    return imgs2

4.关于yuv文件转换的一些函数

"""
import os

import cv2
from matplotlib import pyplot as plt

"""
Threshold A is designed to react to abrupt changes in the input signal and
Threshold B is designed to react to continuous changes in the input signal

Threshold A and Threshold B were determined and used in further experiments. They are as follows:
ThresholdA = 5·σ, ThresholdB = 10·σ.
"""
import numpy as np

def split420(single_frame_data, height, width):
    t = height * width
    Y = single_frame_data[:t].reshape(height, width)
    U = single_frame_data[t:t//4*5].reshape(height//2, width//2)
    V = single_frame_data[t // 4 * 5:].reshape(height // 2, width // 2)
    return Y, U, V

def pixel_yuv2rgb(y, u, v):
    '''
    https://zhuanlan.zhihu.com/p/111069275
    像素的转换和转换公式
    '''
    r = y + (1.370705 * (v - 128))
    g = y - (0.698001 * (v - 128)) - (0.337633 * (u - 128))
    b = y + (1.732446 * (u - 128))
    return r, g, b
def cvt_yuv2rgb(Y, U, V):
    Y = Y.astype(np.float32)
    U = U.astype(np.float32)
    V = V.astype(np.float32)
    rgb = []
    h, w = Y.shape

    # for i in range(h):
    #     ii = i // 2
    #     for j in range(w):
    #         jj = j // 2
    #         y = Y[i, j]
    #         u = U[ii, jj]
    #         v = V[ii, jj]
    #         r, g, b = pixel_yuv2rgb(y, u, v)
    #         rgb.append([r,g,b])

    U2 = np.zeros([h, w], dtype=np.float32)
    U2[0::2, 0::2] = U
    U2[0::2, 1::2] = U
    U2[1::2, 0::2] = U
    U2[1::2, 1::2] = U
    V2 = np.zeros([h, w], dtype=np.float32)
    V2[0::2, 0::2] = V
    V2[0::2, 1::2] = V
    V2[1::2, 0::2] = V
    V2[1::2, 1::2] = V
    r, g, b = pixel_yuv2rgb(Y, U2, V2)
    rgb = cv2.merge([r,g,b])
    rgb = np.array(rgb).reshape(h, w, 3)
    rgb = np.clip(rgb, 0, 255).astype(np.uint8)
    return rgb

def cvt_single(single_frame_data, height, width, show_im=0):
    '''
    将一个YUV420 的数据解析为 Y， U， V ，并转换为rgb
    '''
    Y, U, V = split420(single_frame_data, height, width)
    if show_im:
        plt.figure()
        plt.subplot(131)
        plt.imshow(Y, cmap='gray')

        plt.subplot(132)
        plt.imshow(U)

        plt.subplot(133)
        plt.imshow(V)
        plt.show()

    rgb = cvt_yuv2rgb(Y, U, V)
    if show_im:
        plt.figure()
        plt.imshow(rgb)
        plt.show()
    return Y,U,V,rgb

5.主函数
测试图像：https://download.csdn.net/download/tywwwww/86264119

def cal_psnr(im1, im2):
    t = 10 * np.log10(255 ** 2 / np.mean(((im1 - im2) ** 2)))
    return t
if __name__ == "__main__":
    file = r'D:\dataset\video_test\bus_cif\bus_cif.yuv' # 90frame
    dirY = r'D:\dataset\video_test\bus_cif\Y'
    dirRGB = r'D:\dataset\video_test\bus_cif\RGB'
    if not os.path.exists(dirY):
        os.makedirs(dirY)
    if not os.path.exists(dirRGB):
        os.makedirs(dirRGB)

    data = np.fromfile(file, dtype=np.uint8)
    # print(len(data)/90)
    height = 288
    width = 352
    single_frame_size = height * width * 3 // 2
    frame_num = len(data) // single_frame_size

    fps = 30
    videoWriteRGB = cv2.VideoWriter(file[:-4] + 'RGB.avi', cv2.VideoWriter_fourcc('I', '4', '2', '0'), fps, (width, height))
    sigma = 20 # noise level

    Ys = []
    rgbs = []
    Yns = []
    rgbns = []
    start = 0
    for i in range(frame_num):
        print(i)
        single_frame_data = data[start:start+single_frame_size]
        start += single_frame_size
        Y, U, V, rgb = cvt_single(single_frame_data, height, width, show_im=0)
        # write
        videoWriteRGB.write(rgb[...,::-1])
        # add gauss noise
        Yn = Y + np.random.normal(0, sigma, Y.shape)
        rgbn = rgb + np.random.normal(0, sigma, rgb.shape)
        #Yn = np.clip(Yn, 0, 255).astype(np.uint8)         # 如果进行clip，降噪效果会变差
        #rgbn = np.clip(rgbn, 0, 255).astype(np.uint8)
        # save
        # cv2.imwrite(dirY + '\\'+str(i)+'.png', Y)
        # cv2.imwrite(dirRGB + '\\' + str(i) + '.png', rgb[...,::-1])
        # cv2.imwrite(dirY + '\\' + str(i) + 'n.png', Yn)
        # cv2.imwrite(dirRGB + '\\' + str(i) + 'n.png', rgbn[..., ::-1])
        # list
        Ys.append(Y)
        rgbs.append(rgb)
        Yns.append(Yn)
        rgbns.append(rgbn)

    thrA = 10 * sigma
    thrB = 30 * sigma

    USE_Y = 0
    USE_RGB = 1
    if USE_Y:
        denoised = denoise_adaptive_temp_average(Yns, thrA, thrB)
        denoised2 = denoise_temp_average(Yns, 2)
        denoised3 = denoise_adaptive_temp_pixel_average(Yns, r=5, thr1=4*sigma, thr2=0.2)

        NT = 99

        snr1 = cal_psnr(Ys[NT], Yns[NT])
        snr2 = cal_psnr(Ys[NT], denoised[NT])
        print(snr1, snr2)
        snr1 = cal_psnr(Ys[NT], Yns[NT])
        snr2 = cal_psnr(Ys[NT], denoised2[NT])
        print(snr1, snr2)
        snr1 = cal_psnr(Ys[NT], Yns[NT])
        snr2 = cal_psnr(Ys[NT], denoised3[NT])
        print(snr1, snr2)

        plt.figure()
        plt.subplot(131)
        plt.imshow(Ys[NT])
        plt.subplot(132)
        plt.imshow(Yns[NT])
        plt.subplot(133)
        plt.imshow(denoised3[NT])
        plt.show()

        for i in np.arange(0, frame_num):
            cv2.imwrite(dirY + '\\' + str(i) + 'dn.png', denoised[i])
        for i in np.arange(0, frame_num):
            cv2.imwrite(dirY + '\\' + str(i) + 'dn2.png', denoised2[i])
        for i in np.arange(0, frame_num):
            cv2.imwrite(dirY + '\\' + str(i) + 'dn3.png', denoised3[i])


    elif USE_RGB:
        denoised = denoise_adaptive_temp_average(rgbns, thrA, thrB)
        denoised2 = denoise_temp_average(rgbns, 2)
        denoised3 = denoise_adaptive_temp_pixel_average(rgbns, r=5, thr1=4*sigma, thr2=0.2)

        for i in np.arange(0, frame_num):
            cv2.imwrite(dirRGB + '\\' + str(i) + 'dn.png', denoised[i][..., ::-1])
        for i in np.arange(0, frame_num):
            cv2.imwrite(dirRGB + '\\' + str(i) + 'dn2.png', denoised2[i][..., ::-1])
        for i in np.arange(0, frame_num):
            cv2.imwrite(dirRGB + '\\' + str(i) + 'dn3.png', denoised3[i][..., ::-1])

方便具体分析:

"""
采用比较少的帧，更具体的分析下时域滤波的效果

QCIF Format
(176x144) = 25344

CIF Format
(352x288) = 101376
"""
import os

import cv2
import numpy as np
from matplotlib import pyplot as plt

from denoise_video.denoise_adaptive_temporal_averaging import denoise_temp_average, cal_psnr, \
    denoise_adaptive_temp_pixel_average

if __name__ == "__main__":
    # salesman_qcif
    dirY = r'D:\dataset\video_test\bus_cif\Y'
    dirRGB = r'D:\dataset\video_test\bus_cif\RGB'

    dir = dirRGB
    k = 100
    n = 20
    sigma = 20

    imgs = []
    imgns = []
    for i in range(k-n,k+n+1):
        file = os.path.join(dir, str(i)+'.png')
        img = cv2.imread(file)
        img = img[:, :, ::-1]
        imgn = img + np.random.normal(0, sigma, img.shape)
        imgs.append(img)
        imgns.append(imgn)
    r = 4
    imgs2 = denoise_temp_average(imgns, r)
    imgs3 = denoise_adaptive_temp_pixel_average(imgns, 2*r, thr1=2*sigma, thr2=0.2)


    psnr = cal_psnr(imgs[n], imgns[n])
    psnr2 = cal_psnr(imgs[n], imgs2[n])
    psnr3 = cal_psnr(imgs[n], imgs3[n])
    print('psnr :', psnr, psnr2, psnr3)

    plt.figure()
    plt.subplot(231)
    plt.imshow(imgs[n], cmap='gray')

    plt.subplot(232)
    plt.imshow(np.clip(imgns[n], 0, 255).astype(np.uint8))

    plt.subplot(233)
    plt.imshow(imgs2[n])
    plt.subplot(234)
    plt.imshow(imgs3[n])
    plt.show()