noise estimation : Fast Method for Noise Level Estimation and Integrated Noise Reduction

1.该篇论文介绍了噪声估计和去除方法。

这里只简单介绍其噪声估计方法，想深入了解的可以查看原文。

2.首先高斯分布的标准差原理（高中知识学过）

3.这篇论文的噪声估计方法是朴素，简单的。

首先要明白，如果是没有纹理细节的平坦区域，比如白纸，那么这块区域本身的方差就是噪声水平。对于一个图像，并没有放置白纸，但是我们可以通过平坦区域来近似。 因此一个噪声估计方法就是，首先找到平坦区域，然后计算方差。

4.如何计算平坦区域

设置参数sigma_max。 根据噪声水平，该参数应该设置为实际噪声的3倍及以上。实际选择的时候这个参数很重要，作者说主要针对noise level < 12的场景，因此作者选取的sigma max 为 30， 35 左右。
然后计算邻域距离

如果 di 大于 sigma_max, 则该像素不属于平坦区域。
这样检测到的平坦区域如下图：
a是原图，b是上述方法检测得到的区域（黑色属于平坦区域）， c是另一篇论文的检测方法（Adaptive Temporal Filtering for CFA Video Sequences），可以看出两种方法得到的结果类似，因此证明该方法的有效和简答。

5.检测flat pixel后，计算 d的直方图

如果某个像素满足了flat pixel， 那么和邻域的8个差值就是 di, 所有flat pixel的di 的直方图，
然后可以拟合得到高斯分布的方差，即噪声水平。或者直方图的68分位

最后，直方图的68分位作为噪声水平的估计

6. code

# Fast Method for Noise Level Estimation and Integrated Noise Reduction
import glob
import os
import sys

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

import time

from scipy.optimize import least_squares
from scipy.stats import norm


def noise_estimate_fast_method(img):
    '''
    img: uint 8 dtype
    return: mu, sigma, [0-255]
    '''
    im = img.astype(np.float32)
    h, w = im.shape[:2]

    g = img
    kernel = 1
    g1 = np.pad(g, ((kernel, kernel), (kernel, kernel)), 'reflect')

    lap = []
    for i in range(2 * kernel + 1):
        for j in range(2 * kernel + 1):
            if not ((i == kernel) and (j == kernel)):
                lap.append(np.abs(g1[i:i + h, j:j + w] - g))

    lap = np.dstack(lap)

    sigma_max = 35
    mask = (lap > sigma_max).astype(np.uint8)
    mask = (np.sum(mask, -1) == 0)

    data = lap[mask]

    hist = cal_hist(data)[:sigma_max + 1]
    hist = hist / hist.sum()

    hist2 = np.zeros(2 * sigma_max + 1)
    a = hist[1:] / 2
    hist2[:sigma_max] = a[::-1]
    hist2[sigma_max] = hist[0]
    hist2[sigma_max + 1:] = a

    cumhist = cal_cum_hist(hist2)

    xx = np.arange(-sigma_max, sigma_max + 1, 1)
    show_fig = 0
    if show_fig:
        plt.figure()
        plt.plot(xx, hist2, 'g-')
        plt.plot(xx[sigma_max:], hist2[sigma_max:] * 2, 'r+')
        plt.plot(xx, cumhist, 'b--')
        plt.show()

    mu, sigma = fit_func_gauss(xx, cumhist)
    return mu, sigma
def demo():
    file = r'D:\code\denoise\sidd_python\mean-fusion-to-denoise-image-master\origin_gray.png'
    im = cv2.imread(file)[..., ::-1]
    #im = cv2.imread('./lena.png')
    im = im.astype(np.float32)

    noise_level = [0, 0.5, 1, 2, 5, 10, 15, 20, 30, 45, 50, 70]

    noise_level = [0, 0.5, 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15]
    for level in noise_level:
        sigma = level

        im_noise = im + np.random.randn(*im.shape) * sigma

        start = time.time()
        mu, est_level = noise_estimate_fast_method(im_noise[..., 1])
        end = time.time()
        time_elapsed = end -start

        str_p = "Time: {0:.4f}, Ture Level: {1:6.4f}, Estimated Level: {2:6.4f}"
        print(str_p.format(time_elapsed, level, est_level), 'mean:', mu)



def cal_hist(img):
    '''

    :param img: single channel image
    :return: hist
    '''
    hist, hist_edge = np.histogram(img, bins=256, range=[0, 256])
    # for i in range(256):
    #     hist[i] = np.count_nonzero(img==i)
    # print(hist_edge)
    return hist
def cal_cum_hist(hist):
    cum_hist = np.zeros_like(hist)
    for i in np.arange(1, len(hist)):
        cum_hist[i] = cum_hist[i-1] + hist[i]
    return cum_hist

def fit_func_gauss(x, y):
    '''
    指定函数拟合 ，然后插值xx
    高斯分布累计函数拟合
    '''
    def func(coef, x, y):
        w = 1 #/ (np.sqrt((y * (1 - y))) + 1e-2)
        t = y - norm.cdf(x, loc=coef[0], scale=coef[1])
        t = t * w
        return t.reshape(-1)

    coef = np.array([0, 1])
    res = least_squares(func, coef, args=(x, y), verbose=0)
    if res.success:
        mu, nu = res.x
    else:
        print('cal fail !')
        mu, nu = 0, 1
    return mu, nu

if __name__ == "__main__":
    # file = r'data\\cap_frame_100_gamma.png'
    # im = cv2.imread(file)[..., ::-1]
    # im = cv2.imread('./lena.png')
    demo()