泊松噪音

Knuth算法

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

P(x=k)=e−λλkk!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,∞], 而观察函数图像，等效于一半指数函数axa^xax，其中0<a<10 < a < 10<a<1

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

即有

P(twait>tevent)=e−σ⋅tP(t_{wait} > t_{event}) = e^{- \sigma \cdot t}P(twait>tevent)=e−σ⋅t

其中有

Pwait>PeventP_{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>pL > pL>p，ppp 代表着当前 kkk 次的事件概率，而 LLL 则是0次事件出现的概率，根据公式，令 k=0k=0k=0 , P(x=k)=e−λλ00!=e−λP(x=k) = \frac{e^{- \lambda} \lambda ^ 0}{0!} =e^{-\lambda}P(x=k)=0!e−λλ0=e−λ。

其次，我们已经知道了 Pwait>PeventP_{wait} > P_{event}Pwait>Pevent 也就是说后面事件发生的概率必然比0次事件的概率低，所以得到了P0>PkP_0 > P_kP0>Pk 也就是 e−λ>Pke^{-\lambda} > P_ke−λ>Pk

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

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

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

散列生成算法

我这里直接摆上代码吧

def poisson_distribution(lam: float, limits: int):distributions = [] e_lam = math.pow(math.e, -lam)                # e^{-lambda}sum = 0for k in range(0, 100):lambda_k = math.pow(lam, k)                # lambda^kk_factorial = math.factorial(k)              # k!prob = (lambda_k / k_factorial) * e_lam    # poisson distribution when x=ksum = sum + probif limits - 1 <= k:others = 1 - sum                          # long tail incidentsdistributions.append(others)breakif prob > 0:distributions.append(prob)k = k + 1else:breakreturn distributionsdef generate_poisson_list(distribution, size=100):pos = 0poisson_list = []for prob in distribution:max = round(size * prob)     # create a list has size,# convert the probability to a certain length# of arrayp_list = []if max < 1:                  # size is not enough to init a list, skipbreakfor i in range(max):      # assign k to the certain length arrayp_list.append(pos)pos = pos + 1poisson_list.append(p_list)  # append the probability converted arraypoisson_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 = 0p = 1while p >= L:k = k + 1# Generate uniform random number u in [0, 1] and let p ← p × u.p = random.random() * preturn 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 valuenoise_value = 0for k in range(dts):noise_value = noise_value + poisson_value1(pvals)# add noise to original imagetemp = image[i][j] + noise_value - dts * lamb  # compensating if temp > 255:temp = 255if temp < 0:temp = 0# assign noised image to outputoutput[i][j] = tempreturn 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 valuenoise_value = 0for k in range(image[i][j]):noise_value = noise_value + poisson_value2(lamb)# add noise to original imagetemp = noise_valueif temp > 255:temp = 255# assign noised image to outputoutput[i][j] = tempreturn output

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

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

数字图像学笔记 —— 16. 图像退化与复原（自适应滤波之「最小均方差滤波」）
文章目录图像恢复的一般运算过程什么是「最小均方差滤波」实现步骤实现代码最后的结果图像恢复的一般运算过程我们从前几章的基本理论出发,退化信号恢复成原始信号的步骤,可以概括成两步基本公式.对 ...
数字图像学笔记——13. 图像退化与复原（退化函数的评估方法：观察法、实验法、数学建模法与湍流导致的退化）
在对受到多种原因影响的图像进行复原时,我们经常需要先行评估对图像质量产生影响的退化函数,有时甚至需要尝试建模.通过这些手段,能够最大程度上恢复图像上的噪音,并重建高清的图像细节. 文章目录线性位置不 ...
数字图像学笔记——14. 图像退化与复原（线性退化）
文章目录运动导致的退化(线性退化) 水平运动导致的退化垂直运动导致的退化运动导致的退化(线性退化) 在上一章 <数字图像学笔记--13. 图像退化与复原(退化函数的评估方法:观察法.实验法 ...
数字图像学笔记——10. 频域与傅里叶分析方法
频域滤波技术,目前主要使用的有两种类型,一种是傅立叶变换技术,还有一种是小波分析.基本逻辑就是把原始信号映射到频率空间中,使得在时域空间无法处理的信号,得以在另外一种空间体系下能够被较有效的处理. 目 ...
数字图像学笔记——3.彩色转黑白
文章目录一些说明关于示例代码关于依赖环境关于教材灰度图.亮度图(Gray Image) 彩色图转灰度图一般亮度转换(luminosity method) 亮度优先转换(luminosity ...
数字图像学笔记——4. 直方图计算、线性变换、对数变换、Gamma变换
文章目录灰度直方图(Gray Histogram) 直方图的计算方法简单的图像转换方法线性变换 / 图像翻转(Image Nagatives) 对数变换(Log Transformation) ...
iOS 使用CoreAudio生成白噪音频数据
iOS 使用CoreAudio生成白噪音频数据(空白音频) /// 生成一段白噪音频数据/// - Parameters:/// - startFrm: 开始frame/// - nFrames: 持 ...
动手学习深度学习-跟李沐学AI-自学笔记(1)
个人学习笔记,如有错误欢迎指正! 预备课课程必备网站:[课程主页][https://courses.d2l.ai/zh-v2],[教材][https://zh-v2.d2l.ai/],[课程论坛讨论 ...
李宏毅-2023春机器学习 ML2023 SPRING-学习笔记：3/24机器如何生成图像
目录 3/24 机器如何生成图像速览图像生成常见模型浅谈图像生成模型 Diffusion Model 原理 Stable Diffusion.DALL-E.Imagen 背后共同的套路 Varia ...

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

文章目录

泊松噪音

Knuth算法

散列生成算法

生成泊松噪音的图像

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

最新文章

热门文章