[数字图像学笔记] 5.噪音生成——泊松噪音(二)

OpenCV与数字图像处理 专栏收录该内容
14 篇文章 0 订阅

泊松噪音

Knuth算法

首先,回顾泊松分布的函数:

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

其中, λ \lambda λ是期望值,而 e − λ e^{-\lambda} eλ则是单调递减的指数函数,而我们所需要关心的函数区间是 ∈ [ 0 , ∞ ] \in[0, \infty] [0,], 而观察函数图像,等效于一半指数函数 a x a^x ax,其中 0 < a < 1 0 < a < 1 0<a<1
在这里插入图片描述
另一方面,根据之前的关于 “泊松等待” 里介绍的,对于已发生的事件A,在接下来的时间里,随着时间增加,事件发生概率呈指数级下降。

即有

P ( t w a i t > t e v e n t ) = e − σ ⋅ t P(t_{wait} > t_{event}) = e^{- \sigma \cdot t} P(twait>tevent)=eσt

其中有

P w a i t > P e v e n t P_{wait} > P_{event} Pwait>Pevent

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

以图来说明这个过程:

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

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

在这里插入图片描述
对于第n个 Δ t \Delta t Δt内,这个概率就会更小,几乎忽略不计了

在这里插入图片描述
而我们的解体思想,就是在不同的事件时间内,随机生成一个 [ 0 , 1 ] [0, 1] [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 L > p L>p p p p 代表着当前 k k k 次的事件概率,而 L L L 则是0次事件出现的概率,根据公式,令 k = 0 k=0 k=0 , P ( x = k ) = e − λ λ 0 0 ! = e − λ P(x=k) = \frac{e^{- \lambda} \lambda ^ 0}{0!} =e^{-\lambda} P(x=k)=0!eλλ0=eλ

其次, 我们已经知道了 P w a i t > P e v e n t P_{wait} > P_{event} Pwait>Pevent 也就是说后面事件发生的概率必然比0次事件的概率低,所以得到了 P 0 > P k P_0 > P_k P0>Pk 也就是 e − λ > P k e^{-\lambda} > P_k eλ>Pk

最后,随着第 k k k次循环的 P k = P k − 1 × P k − 2 × P k − 3 × . . . × 1 P_k = P_{k-1} \times P_{k-2} \times P_{k-3} \times ... \times 1 Pk=Pk1×Pk2×Pk3×...×1

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

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

散列生成算法

我这里直接摆上代码吧

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

在这里插入图片描述

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

  • 0
    点赞
  • 0
    评论
  • 2
    收藏
  • 一键三连
    一键三连
  • 扫一扫,分享海报

©️2020 CSDN 皮肤主题: 程序猿惹谁了 设计师:白松林 返回首页
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、C币套餐、付费专栏及课程。

余额充值