数字图像学笔记——7. 噪音生成（泊松噪音生成方法）

泊松噪音

Knuth算法

首先，回顾泊松分布的函数：

$$P(x=k)=\frac{e^{-\lambda }{\lambda }^k}{k!}$$

其中，$\lambda$是期望值，而$e^{-\lambda }$则是单调递减的指数函数，而我们所需要关心的函数区间是$\in [0,\infty ]$, 而观察函数图像，等效于一半指数函数$a^x$，其中$0<a<1$

另一方面，根据之前的关于 “泊松等待” 里介绍的，对于已发生的事件A，在接下来的时间里，随着时间增加，事件发生概率呈指数级下降。

即有

$$P(t_{wait}>t_{event})=e^{-\sigma \cdot t}$$

其中有

$$P_{wait}>P_{event}$$

这个限制条件存在。那么，假设打开快门的一瞬间，什么事也没有发生的概率必然是1， 而随着时间的流逝，比如经过3个 $\Delta t$ 之后，出现一个光子击中了像素传感器这种事，也可能出现所有的 $\Delta t$都结束后，一个光子都没有击中传感器的概率存在。

以图来说明这个过程：

对于刚打开快门时，发生事件的概率会非常大，那么我们生成的光子就有很大概率落在了绿色的事件时间窗口与红色的概率区间范围内，即第一个 $\Delta t$ 内，有一个点同时出现在红色区域。

如果没有击中像素传感器，那么在第二个时间窗口 $\Delta t$ 内，发生这件事的概率就会迅速下降。

对于第n个 $\Delta t$内，这个概率就会更小，几乎忽略不计了

而我们的解体思想，就是在不同的事件时间内，随机生成一个 $[0,1]$ 的随机数，用来表示这个点可能在第几次的时候击中传感器。

唐纳德·努斯祖师爷给了一个算法，可以用来模拟这个过程

algorithm poisson random number (Knuth):
    init:
         Let L ← exp(−λ), k ← 0 and p ← 1.
    do:
         k ← k + 1.
         Generate uniform random number u in [0,1] and let p ← p × u.
    while p > L.
    return k − 1.

我这里进行一下说明，这样你就明白了这个算法的伟大之处。

首先 ，这个算法的限制条件是 $L>p$，$p$ 代表着当前 $k$ 次的事件概率，而 $L$ 则是0次事件出现的概率，根据公式，令 $k=0$ , $P(x=k)=\frac{e^{-\lambda }{\lambda }^0}{0!}=e^{-\lambda }$。

其次， 我们已经知道了 $P_{wait}>P_{event}$ 也就是说后面事件发生的概率必然比0次事件的概率低，所以得到了$P_0>P_k$ 也就是 $e^{-\lambda }>P_k$

最后，随着第$k$次循环的$P_k=P_{k-1}\times P_{k-2}\times P_{k-3}\times ...\times 1$

而$p=p\times u$ 因上式得到了 $P_k=u^k$, 又因为$u\in (0,1)$ 所以这里其实得到的是单调递减的指数函数。所以，问题最后被简化成了每一次生成的$P_k$是否比第0次事件小，如果 $P_k>P_0$ 说明当前发生的事件可能是一个概率极低的事件。

说实话，初次看到这个算法确实不是很直观，需要发挥一些想象力，并且查阅了大量资料。所以我这里再提供一个我想到的散列生成算法，而且比较直观，你可以根据自己的需要实现。

散列生成算法

我这里直接摆上代码吧

def poisson_distribution(lam: float, limits: int):
    distributions = [] 
    e_lam = math.pow(math.e, -lam)  				# e^{-lambda}
    sum = 0

    for k in range(0, 100):
        lambda_k = math.pow(lam, k) 		 		# lambda^k
        k_factorial = math.factorial(k)  			# k!
        prob = (lambda_k / k_factorial) * e_lam  	# poisson distribution when x=k
        sum = sum + prob

        if limits - 1 <= k:
            others = 1 - sum  						# long tail incidents
            distributions.append(others)
            break

        if prob > 0:
            distributions.append(prob)
            k = k + 1
        else:
            break

    return distributions


def generate_poisson_list(distribution, size=100):
    pos = 0
    poisson_list = []

    for prob in distribution:
        max = round(size * prob)  	# create a list has size,
        							# convert the probability to a certain length
        							# of array
        p_list = []

        if max < 1:					# size is not enough to init a list, skip
            break

        for i in range(max):		# assign k to the certain length array
            p_list.append(pos)

        pos = pos + 1
        poisson_list.append(p_list)	# append the probability converted array

    poisson_list = [i for elem in poisson_list for i in elem] # flat the ragged arrays to an 1-d array!
    return poisson_list

其中，poisson_distribution 用来生成泊松分布，然后我们创建一个限定长度的散列，把概率问题通过generate_poisson_list转换为查表问题，接下来就只需要使用这个函数来随机生成一个给定范围的数，然后检查它落在了哪个事件区间里，嗯就是这么简单粗暴……

def poisson_value2(poisson_list):
    number = random.random() * (len(poisson_list) - 1)  # to avoid the situation when a[i]
    													# i = len(a)
    number = round(number)
	return poisson_list[number]  						# k = poisson_list[number]

knuth大神的代码实现则大概是这样的味道：

def poisson_value1(lamb):
    L = math.exp(-lamb)
    k = 0
    p = 1

    while p >= L:
        k = k + 1
        # Generate uniform random number u in [0, 1] and let p ← p × u.
        p = random.random() * p

    return k - 1

至于算法效率，我还没来得及对比，如果你感兴趣可以对比一下。接下来，我们试着生成包含泊松噪音的图片吧。

生成泊松噪音的图像

在获得了泊松噪音函数之后，为了获得泊松噪音添加后的图像，我们可以有两种不同的图像基本处理方法：

原始图像 + 噪音图像 - 溢出补偿:

def poisson_noise1(image, pvals, dts, lamb):
    output = np.zeros(image.shape, np.uint8)

    for i in range(image.shape[0]):
        for j in range(image.shape[1]):

            # get the poisson value
            noise_value = 0
            for k in range(dts):
                noise_value = noise_value + poisson_value1(pvals)

            # add noise to original image
            temp = image[i][j] + noise_value - dts * lamb  # compensating 
            if temp > 255:
                temp = 255
            if temp < 0:
            	temp = 0

            # assign noised image to output
            output[i][j] = temp

    return output

得出的结果是这样的：

根据光通量重新生成照片

def poisson_noise2(image, lamb):
    output = np.zeros(image.shape, np.uint8)

    for i in range(image.shape[0]):
        for j in range(image.shape[1]):

            # get the poisson value
            noise_value = 0
            for k in range(image[i][j]):
                noise_value = noise_value + poisson_value2(lamb)

            # add noise to original image
            temp = noise_value
            if temp > 255:
                temp = 255

            # assign noised image to output
            output[i][j] = temp

    return output

至此，图像学最常见的三种噪音算法演示完毕。