噪声论文：An Efficient Statistical Method for Image Noise Level Estimation

4. 原文作者有matlab程序，也可以结合大佬的python程序来理解，python程序来源：https://github.com/zsyOAOA/noise_est_ICCV2015

#!/usr/bin/env python
# -*- coding:utf-8 -*-
# Power by Zongsheng Yue 2019-01-07 14:36:55
import sys

import cv2
import numpy as np
from matplotlib import pyplot as plt
from skimage import img_as_float
import time

def im2patch(im, pch_size, stride=8):
    '''
    Transform image to patches.
    Input:
        im: 3 x H x W or 1 X H x W image, numpy format
        pch_size: (int, int) tuple or integer
        stride: (int, int) tuple or integer
    '''
    if isinstance(pch_size, tuple):
        pch_H, pch_W = pch_size
    elif isinstance(pch_size, int):
        pch_H = pch_W = pch_size
    else:
        sys.exit('The input of pch_size must be a integer or a int tuple!')

    if isinstance(stride, tuple):
        stride_H, stride_W = stride
    elif isinstance(stride, int):
        stride_H = stride_W = stride
    else:
        sys.exit('The input of stride must be a integer or a int tuple!')


    C, H, W = im.shape
    num_H = len(range(0, H-pch_H+1, stride_H))
    num_W = len(range(0, W-pch_W+1, stride_W))
    num_pch = num_H * num_W
    pch = np.zeros((C, pch_H*pch_W, num_pch), dtype=im.dtype)
    kk = 0
    for ii in range(pch_H):
        for jj in range(pch_W):
            temp = im[:, ii:H-pch_H+ii+1:stride_H, jj:W-pch_W+jj+1:stride_W]
            pch[:, kk, :] = temp.reshape((C, num_pch))
            kk += 1

    return pch.reshape((C, pch_H, pch_W, num_pch))

def noise_estimate(im, pch_size=8):
    '''
    Implement of noise level estimation of the following paper:
    Chen G , Zhu F , Heng P A . An Efficient Statistical Method for Image Noise Level Estimation[C]// 2015 IEEE International Conference
    on Computer Vision (ICCV). IEEE Computer Society, 2015.
    Input:
        im: the noise image, H x W x 3 or H x W numpy tensor, range [0,1]
        pch_size: patch_size
    Output:
        noise_level: the estimated noise level
    '''

    if im.ndim == 3:
        im = im.transpose((2, 0, 1))
    else:
        im = np.expand_dims(im, axis=0)

    # image to patch
    pch = im2patch(im, pch_size, 3)  # C x pch_size x pch_size x num_pch tensor
    num_pch = pch.shape[3]
    pch = pch.reshape((-1, num_pch))  # d x num_pch matrix
    d = pch.shape[0]

    mu = pch.mean(axis=1, keepdims=True)  # d x 1
    X = pch - mu
    sigma_X = np.matmul(X, X.transpose()) / num_pch
    sig_value, _ = np.linalg.eigh(sigma_X)
    sig_value.sort()

    t = 0
    for ii in range(-1, -d-1, -1):
        tau = np.mean(sig_value[:ii])
        if np.sum(sig_value[:ii]>tau) == np.sum(sig_value[:ii] < tau):
            t = np.sqrt(tau)
            break
    return sig_value, t

if __name__ == '__main__':

    # file = r'D:\code\data\kodak_IDR\kodim24_gaussian_NOISY_25.0.png'
    #file = r'D:\code\data\kodak2\kodim12_gaussian_0.1.png'
    # file = r'D:\dataset\noise_img\cap_frame_2016_D65_3000lx_raw.jpg'
    # file = r'D:\dataset\noise_img\cap_frame_0001_raw.png'
    file = r'a.jpg'
    im = cv2.imread(file)
    #im = cv2.imread('./lena.png')
    im = img_as_float(im)

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

    for level in noise_level:
        sigma = level / 255

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

        start = time.time()
        eigval, est_level = noise_estimate(im_noise, 8)
        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*255))
        plt.plot(255* np.sqrt(eigval)[::-1][:])
    plt.grid()

    plt.show()

5. 作者提到 关于patch size的选择

patch_size比较大时，beta值越大，意味着可以允许比较多的主要特征比例。
但是另一方面， patch的数量越大的时候，冗余特征值的分布才符合

因此size 和 number 需要一个折衷。