扩散模型在图像生成中的应用：从真实样例到逼真图像的奇妙转变

一、扩散模型

        扩散模型的起源可以追溯到热力学中的扩散过程。热力学中的扩散过程是指物质从高浓度往低浓度的地方流动，最终达到一种动态的平衡。这个过程就是一个扩散过程。 

        在深度学习领域中，扩散模型（diffusion models）是深度生成模型中新的SOTA。 扩散模型在图片生成任务中超越了原SOTA：GAN，并且在诸多应用领域都有出色的表现，如计算机视觉、NLP、波形信号处理、多模态建模、分子图建模、时间序列建模、对抗性净化等。

  1、扩散模型原理介绍

在深度学习中，生成模型的目标是根据给定的样本（训练数据）生成新样本。首先给定一批训练数据，假设其服从某种复杂的真实分布，则给定的训练数据可视为从该分布中采样的观测样本。如果能从这批观测样本中估计出训练数据的真实分布，就可以从该分布中源源不断地采样出新的样本。

生成模型作用是估计训练数据的真实分布，并将其假定为q(x)" role="presentation" style="position: relative;">q(x)$q(x)$，这个过程称为拟合网络。

如何确定估计的分布q(x)" role="presentation" style="position: relative;">q(x)$q(x)$和真实分布的差距？

 要求所有的训练数据样本采样自q(x)" role="presentation" style="position: relative;">q(x)$q(x)$的概率最大，主要思想就是统计学中的最大似然估计思想，这也是生成模型的基本思想之一，因此生成模型的学习目标就是对训练数据的分布进行建模。

        公认最早的扩散模型DDPM（Denoing Diffusion Probalilistic Model）的扩散模型原理就由此而来，扩散过程如下图所示，具体分为前向过程和反向过程两部分：

1.1、扩散模型前向过程

 前向过程是给数据中添加噪声的过程，假设给定一批训练数据，数据分布为  x0∼q(x0)" role="presentation" style="position: relative;">x0∼q(x0)$x_0\sim q(x_0)$   ，其中，0表示初始状态，即还没有开始扩散，即上图的原图（小狗）。

如上所述，将前向加噪声过程分为离散的多个时间步，在每一个时间步  ，给上一个时间步 的数据xt−1" role="presentation" style="position: relative;">xt−1$x_{t-1}$ 添加高斯噪声，从而生成带有噪声（简称“带噪”）的数据  xt" role="presentation" style="position: relative;">xt$x_t$ ，同时数据xt" role="presentation" style="position: relative;">xt$x_t$ 也会被送入下一个时间步  t+1" role="presentation" style="position: relative;">t+1$t+1$  以继续添加噪声。其中，噪声的方差由一个位于区间（0，1）的固定值   和当前时刻“带噪”的数据分布确定。在反复迭代和加噪（即添加噪声） 次之后，只要  足够大，根据马尔可夫链的性质就可以得到纯随机噪声分布的数据，即类似稳定墨水系统的状态。

 马尔可夫链：是马尔可夫过程的原始模型，也是一个表示状态转移的离散随机过程。该离散过程具有“无记忆”的性质，即下一状态的概率分布仅由当前状态表示，而与之前的所有状态无关，同时只要时间序列足够长，即状态转移的次数足够多，最终的概率分布将趋于稳定。 

 1.1.1、扩散模型公式

从时间步  到时间步   的单步扩散加噪声过程的数学表达式如下：

q(xt|xt−1)=N(xt;1−βtxt−1,βtI)" role="presentation" style="position: relative;">q(xt|xt−1)=N(xt;1−βt−−−−−√xt−1,βtI)$q(x_t|x_{t-1})=N(x_t;\sqrt{1-{\beta }_t}{x}_{t-1},{\beta }_{t}I)$                                                                                                                (1)

最终的噪声分布数学表达式如下：

q(x1:T|x0)=∏t=1Tq(xt|xt−1)" role="presentation" style="position: relative;">q(x1:T|x0)=∏Tt=1q(xt|xt−1)$q(x_{1:T}|x_0)=\prod _{t=1}^{T}q(x_t|x_{t-1})$                                                                                                                          (2)

1.2、扩散模型反向过程

         前向过程是将数据噪声化的过程（增加噪声），反向过程则是“去噪”的过程，即从随机噪声中迭代恢复出清晰数据的过程。从采样自高斯噪声xt∼N(0,1)" role="presentation" style="position: relative;">xt∼N(0,1)$x_t\sim N(0,1)$  的一个随机噪声中恢复出原始数据  x0" role="presentation" style="position: relative;">x0$x_0$  ，就需要知道反向过程中每一步的图像分布状态转移。DDPM也将这个方向过程定义为马尔可夫链，由一系列用神经网络参数化的高斯分布组成的，也就是需要训练的扩散模型。

        从时间步    到时间步    的单步反向去噪的过程，公式如下所示：

q(xt−1|xt,x0)=N(xt−1;μ~(xt,x0),βt~I)" role="presentation" style="position: relative;">q(xt−1|xt,x0)=N(xt−1;μ~(xt,x0),βt~I)$q(x_{t-1}|x_t,x_0)=N(x_{t-1};\tilde{\mu }(x_t,x_0),\widetilde{{\beta }_t}I)$

        由于反向过程的每一步都是参数化的高斯分布，因此可以分布求高斯分布的均值和方差。经过贝叶斯公式推导

 q(xt−1|xt,x0)=q(xt,xt−1,x0)q(xt,x0)=q(xt|xt−1,x0)⋅q(xt−1,x0)q(xt|x0)⋅q(x0)=q(xt|xt−1,x0)⋅q(xt−1|x0)⋅q(x0)q(xt|x0)⋅q(x0)=q(xt|xt−1,x0)⋅q(xt−1|x0)q(xt|x0)∝exp⁡(−12((xt−αtxt−1)2βt+(xt−1−α¯t−1x0)21−α¯t−1−(xt−α¯tx0)21−α¯t))) \textcolorblue;because (1) and (2)=exp⁡(−12(xt2−2αtxtxt−1+αtxt−12βt+xt−12−2α¯t−1x0xt−1+α¯t−1x021−α¯t−1−(xt−α¯tx0)21−α¯t)))=exp⁡(−12(\textcolorred(αtβt+11−α¯t−1))xt−12−\textcolorgreen(2αtβtxt+2α¯t−11−α¯t−1x0)xt−1+C(xt,x0))" role="presentation" style="position: relative;">q(xt−1|xt,x0)=q(xt,xt−1,x0)q(xt,x0)=q(xt|xt−1,x0)⋅q(xt−1,x0)q(xt|x0)⋅q(x0)=q(xt|xt−1,x0)⋅q(xt−1|x0)⋅q(x0)q(xt|x0)⋅q(x0)=q(xt|xt−1,x0)⋅q(xt−1|x0)q(xt|x0)∝exp(−12((xt−α−−√txt−1)2βt+(xt−1−α¯t−1−−−−√x0)21−α¯t−1−(xt−α¯t−−√x0)21−α¯t))) \textcolorblue;because (1) and (2)=exp(−12(x2t−2α−−√txtxt−1+αtx2t−1βt+x2t−1−2α¯t−1−−−−√x0xt−1+α¯t−1x201−α¯t−1−(xt−α¯t−−√x0)21−α¯t)))=exp(−12(\textcolorred(αtβt+11−α¯t−1))x2t−1−\textcolorgreen(2αt−−√βtxt+2α¯t−1−−−−√1−α¯t−1x0)xt−1+C(xt,x0))$\begin{array}{rl}q({\mathbf{x}}_{t-1}|{\mathbf{x}}_t,{\mathbf{x}}_0)& =\frac{q({\mathbf{x}}_t,{\mathbf{x}}_{t-1},{\mathbf{x}}_0)}{q({\mathbf{x}}_t,{\mathbf{x}}_0)}\\ & =\frac{q({\mathbf{x}}_t|{\mathbf{x}}_{t-1},{\mathbf{x}}_0)\cdot q({\mathbf{x}}_{t-1},{\mathbf{x}}_0)}{q({\mathbf{x}}_t|{\mathbf{x}}_0)\cdot q({\mathbf{x}}_0)}\\ & =\frac{q({\mathbf{x}}_t|{\mathbf{x}}_{t-1},{\mathbf{x}}_0)\cdot q({\mathbf{x}}_{t-1}|{\mathbf{x}}_0)\cdot q({\mathbf{x}}_0)}{q({\mathbf{x}}_t|{\mathbf{x}}_0)\cdot q({\mathbf{x}}_0)}\\ & =\frac{q({\mathbf{x}}_t|{\mathbf{x}}_{t-1},{\mathbf{x}}_0)\cdot q({\mathbf{x}}_{t-1}|{\mathbf{x}}_0)}{q({\mathbf{x}}_t|{\mathbf{x}}_0)}\\ & \propto \exp \left(-\frac{1}{2}(\frac{({\mathbf{x}}_t-{\sqrt{\alpha }}_t{\mathbf{x}}_{t-1}{)}^2}{{\beta }_t}+\frac{({\mathbf{x}}_{t-1}-\sqrt{{\overline{\alpha }}_{t-1}}{\mathbf{x}}_0{)}^2}{1-{\overline{\alpha }}_{t-1}}-\frac{({\mathbf{x}}_t-\sqrt{{\overline{\alpha }}_t}{\mathbf{x}}_0{)}^2}{1-{\overline{\alpha }}_t)})\right)\text{textcolor}blue\text{;because (1) and (2)}\\ & =\exp \left(-\frac{1}{2}(\frac{{\mathbf{x}}_t^2-2{\sqrt{\alpha }}_t{\mathbf{x}}_t{\mathbf{x}}_{t-1}+{\alpha }_t{\mathbf{x}}_{t-1}^2}{{\beta }_t}+\frac{{\mathbf{x}}_{t-1}^2-2\sqrt{{\overline{\alpha }}_{t-1}}{\mathbf{x}}_0{\mathbf{x}}_{t-1}+{\overline{\alpha }}_{t-1}{\mathbf{x}}_0^2}{1-{\overline{\alpha }}_{t-1}}-\frac{({\mathbf{x}}_t-\sqrt{{\overline{\alpha }}_t}{\mathbf{x}}_0{)}^2}{1-{\overline{\alpha }}_t)})\right)\\ & =\exp \left(-\frac{1}{2}(\text{textcolor}red(\frac{{\alpha }_t}{{\beta }_t}+\frac{1}{1-{\overline{\alpha }}_{t-1})}){\mathbf{x}}_{t-1}^2-\text{textcolor}green(\frac{2\sqrt{{\alpha }_t}}{{\beta }_t}{\mathbf{x}}_t+\frac{2\sqrt{{\overline{\alpha }}_{t-1}}}{1-{\overline{\alpha }}_{t-1}}{\mathbf{x}}_0){\mathbf{x}}_{t-1}+C({\mathbf{x}}_t,{\mathbf{x}}_0)\right)\end{array}$

上式巧妙地通过贝叶斯公式将逆向过程转换为前向过程，且最终得到的概率密度函数和高斯概率密度函数的指数部分，

exp⁡(−(x−μ)22σ2)=exp⁡(−12(\textcolorred1σ2x2−\textcolorgreen2μσ2x+μ2σ2))" role="presentation" style="position: relative;">exp(−(x−μ)22σ2)=exp(−12(\textcolorred1σ2x2−\textcolorgreen2μσ2x+μ2σ2))$\exp (-\frac{(\mathbf{x}-\mu {)}^2}{2{\sigma }^2})=\exp \left(-\frac{1}{2}(\text{textcolor}red\frac{1}{{\sigma }^2}{\mathbf{x}}^2-\text{textcolor}green\frac{2\mu }{{\sigma }^2}\mathbf{x}+\frac{{\mu }^2}{{\sigma }^2})\right)$

相对应。其中C(xt,x0)" role="presentation" style="position: relative;">C(xt,x0)$C({\mathbf{x}}_t,{\mathbf{x}}_0)$，是与xt−1" role="presentation" style="position: relative;">xt−1$x_{t-1}$无关的常数项。令

" role="presentation" style="position: relative;">$$

这个公式表示的是高斯分布（正态分布）的概率密度函数，其中：

- $xt−1$" role="presentation" style="position: relative;">$xt−1$$\${\mathbf{x}}_{t-1}\$$是当前时刻的观测值；
- $xt$" role="presentation" style="position: relative;">$xt$$\${\mathbf{x}}_t\$$ 和 $x0$" role="presentation" style="position: relative;">$x0$$\${\mathbf{x}}_0\$$分别是前一次和前几次的观测值；
- 是预测值，即使用贝叶斯方法预测出来的下一个时刻的值；
- $\textcolorredβ~t$" role="presentation" style="position: relative;">$\textcolorredβ~t$$\$\text{textcolor}red{\tilde{\beta }}_t\$$是预测误差的方差，也就是预测值与实际值之间的差异程度。

这个公式的意义在于，根据已知的前几次观测值和当前时刻的观测值，使用贝叶斯方法预测出下一个时刻的观测值，并计算出预测误差的方差。这个预测值和预测误差的方差可以用来评估模型的准确性和稳定性。

 将绿色部分与红色部分一一对应，一是计算平滑后的数据 $\textcolorredβ~t$" role="presentation" style="position: relative;">$\textcolorredβ~t$$\$\text{textcolor}red{\tilde{\beta }}_t\$$，二是计算平滑后的条件期望 。

  1. 平滑后的数据 $\textcolorredβ~t$" role="presentation" style="position: relative;">$\textcolorredβ~t$$\$\text{textcolor}red{\tilde{\beta }}_t\$$：

$\textcolorredβ~t=1αtβt+11−α¯t−1=βt⋅(1−α¯t−1)αt−α¯t+βt=βt⋅1−α¯t−11−α¯t$" role="presentation" style="position: relative;">$\textcolorredβ~t=1αtβt+11−α¯t−1=βt⋅(1−α¯t−1)αt−α¯t+βt=βt⋅1−α¯t−11−α¯t$$\$\text{textcolor}red{\tilde{\beta }}_t=\frac{1}{\frac{{\alpha }_t}{{\beta }_t}+\frac{1}{1-{\overline{\alpha }}_{t-1}}}=\frac{{\beta }_t\cdot (1-{\overline{\alpha }}_{t-1})}{{\alpha }_t-{\overline{\alpha }}_t+{\beta }_t}={\beta }_t\cdot \frac{1-{\overline{\alpha }}_{t-1}}{1-{\overline{\alpha }}_t}\$$

这里，$αt$" role="presentation" style="position: relative;">$αt$$\${\alpha }_t\$$和 $α¯t−1$" role="presentation" style="position: relative;">$α¯t−1$$\${\overline{\alpha }}_{t-1}\$$分别表示第 $t$" role="presentation" style="position: relative;">$t$$\$t\$$个观测值和前一个观测值的条件期望。$βt$" role="presentation" style="position: relative;">$βt$$\${\beta }_t\$$ 是平滑因子，用于调整数据的平滑程度。

  2. 平滑后的条件期望 ：

$\textcolorgreenμ~(xt,x0)=(2αtβtxt+2α¯t−11−α¯t−1x0)/(αtβt+11−α¯t−1)=(2αtβtxt+2α¯t−11−α¯t−1x0)⋅βt⋅1−α¯t−11−α¯t$" role="presentation" style="position: relative;">$\textcolorgreenμ~(xt,x0)=(2αt√βtxt+2α¯t−1√1−α¯t−1x0)/(αtβt+11−α¯t−1)=(2αt√βtxt+2α¯t−1√1−α¯t−1x0)⋅βt⋅1−α¯t−11−α¯t$$\$\text{textcolor}green\tilde{\mu }({\mathbf{x}}_t,{\mathbf{x}}_0)=(\frac{2\sqrt{{\alpha }_t}}{{\beta }_t}{\mathbf{x}}_t+\frac{2\sqrt{{\overline{\alpha }}_{t-1}}}{1-{\overline{\alpha }}_{t-1}}{\mathbf{x}}_0)/(\frac{{\alpha }_t}{{\beta }_t}+\frac{1}{1-{\overline{\alpha }}_{t-1}})=(\frac{2\sqrt{{\alpha }_t}}{{\beta }_t}{\mathbf{x}}_t+\frac{2\sqrt{{\overline{\alpha }}_{t-1}}}{1-{\overline{\alpha }}_{t-1}}{\mathbf{x}}_0)\cdot {\beta }_t\cdot \frac{1-{\overline{\alpha }}_{t-1}}{1-{\overline{\alpha }}_t}\$$

这里，使用了一个加权平均的方法来计算条件期望。具体来说，我们将当前观测值的预测值乘以平滑因子，然后加上前一个观测值的预测值乘以平滑因子。最后，将结果除以平滑因子的倒数。

这两个公式可以帮助我们在贝叶斯推倒扩散模型中对数据进行平滑处理，从而得到更加稳定和可靠的预测结果。

从上述可以看出，方差是一个定量（扩散过程参数固定），而均值是一个依赖于x0" role="presentation" style="position: relative;">x0$x_0$和xt" role="presentation" style="position: relative;">xt$x_t$的函数，因此需要使用扩散模型来优化参数。（此部分推导过程参考 扩散模型Diffusion Model原理）

贝叶斯分类器基本原理：

当事件B发生时，事件A发生的条件概率可以表示为：

其中，$P(AB)$表示事件A和事件B同时发生的概率；$P(B)$表示事件B发生的概率。根据全概率公式，有：

其中，表示在事件A发生的情况下，事件B发生的概率；$P(B| A)$" role="presentation" style="position: relative;">$P(B| A)$$\$P(B|A)\$$表示在事件A不发生的情况下，事件B发生的概率；表示事件A不发生的概率。将上述两个式子代入上式，得到：

因此，当事件B发生时，事件A发生的概率为1，即 。

二、极大似然估计原理

        极大似然估计的原理目的是利用已知的样本结果，反推最有可能（最大概率）导致这种结果的参数值。 极大似然估计是建立在极大似然原理的基础上的一个统计方法，极大似然估计提供了一种给定观测数据来评估模型参数的方法，即“模型已定，参数未知”。通过若干次试验，观察其结果，利用试验结果得到某个参数值能够使样本出现的概率最大，则称为极大似然估计。

        由于样本集中的样本都是独立同分布的，可以只考虑一类样本集，来估计参数向量θ" role="presentation" style="position: relative;">θ$\theta$，假设已知的样本集为：

D=x1,x2,...,xn" role="presentation" style="position: relative;">D=x1,x2,...,xn$D=x_1,x_2,...,x_n$

         联合概率密度函数 p(D|θ)" role="presentation" style="position: relative;">p(D|θ)$p(D|\theta )$  称为相对于 x1,x2,...,xn" role="presentation" style="position: relative;">x1,x2,...,xn$x_1,x_2,...,x_n$ 的  θ" role="presentation" style="position: relative;">θ$\theta$  似然函数。

l(θ)=p(D|θ)=p(x1,x2,...,xn)=∏i=1np(xi|θ)" role="presentation" style="position: relative;">l(θ)=p(D|θ)=p(x1,x2,...,xn)=∏ni=1p(xi|θ)$l(\theta )=p(D|\theta )=p(x_1,x_2,...,x_n)=\prod _{i=1}^{n}p(x_i|\theta )$

         如果  θ^" role="presentation" style="position: relative;">θ^$\hat{\theta }$  是参数空间中能够使似然函数    最大的值，则 θ^" role="presentation" style="position: relative;">θ^$\hat{\theta }$  应该是 "最可能" 的参数值，那么  θ^" role="presentation" style="position: relative;">θ^$\hat{\theta }$  就是 θ" role="presentation" style="position: relative;">θ$\theta$  的极大似然估计量。它是样本集的函数，记作：

   称为极大似然估计函数估计值。

三、极大似然估计例子

 假设有两个外形完全相同的箱子，1号箱子中有99只白球，1只黑球；2号箱子中只有1只白球，99只黑球。在一次实验中取出黑球，请问是从哪个箱子中取出的？

人们通常会猜测这只黑球最像是从2号箱取出来的，此时描述的“最像”就有“极大似然”的意思，这种想法常称为“极大似然原理”。

最大似然估计是一种参数估计的方法，它是概率论在统计学的应用之一。

最大似然估计是建立在这样的思想上：已知某个参数能使这个样本出现的概率最大，我们当然不会再去选择其他小概率的样本，所以干脆就把这个参数作为估计的真实值。

举个例子，假设我们有一组数据，其中包含了一些人的年龄和收入信息。我们想要估计这些人的平均年龄和平均收入。根据最大似然估计的思想，我们应该尽量让这组数据符合真实情况。也就是说，如果我们假设这些人的平均年龄是30岁，那么当这组数据中有人年龄确实为30岁时，我们就认为这个假设是正确的；如果有人年龄不是30岁，但其他方面都符合条件，我们也认为这个假设是正确的。同样地，如果我们假设这些人的平均收入是5000元/月，那么当这组数据中有人收入确实为5000元/月时，我们就认为这个假设是正确的；如果有人收入不是5000元/月，但其他方面都符合条件，我们也认为这个假设是正确的。

因此，在最大似然估计中，我们需要根据实际情况来选择合适的参数值。

四、优化目标

        扩散模型预测的噪声残差，即要求后向过程中预测的噪声分布与前向过程中施加的噪声分布之间的“距离”最小。扩散模型的最终优化目标的数学表达式如下：

其中，可以看出，在训练DDPM时，只要用一个简单的MSE(均方误差)损失最小化前向过程施加的噪声分布和后向过程预测的噪声分布，就能实现最终的优化目标。

五、从零开始搭建扩散模型（代码实现）

5.1、导入相关的包

导入所需的相关库：采用diffusers提供的DDPM模型接口

注：可能会遇到一些安装包问题，可采用下面的pip源

pip install -i https://pypi.tuna.tsinghua.edu.cn/simple [包名]

pip install -i https://pypi.doubanio.com/simple [包名]

pip install -i http://mirrors.aliyun.com/pypi/simple/ [包名]

import torch
import torchvision
from torch import nn
from torch.nn import functional as F
from torch.utils.data import DataLoader
from diffusers import DDPMScheduler,UNet2DModel
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings("ignore")
 
device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
print(device)

注：本文的代码都是在colab中实现，本地也可以运行。

5.2、数据集测试

本文采用官方提供的一个小型测试集FashMNIST：torchvision.datasets.FashionMNIST（包含10个类别），当然也可以使用其他测试集，这里主要是过一遍流程

dataset = torchvision.datasets.FashionMNIST(root="./FashionMNIST",train=True,download=True,
            transform=torchvision.transforms.ToTensor())
train_dataloader = DataLoader(dataset,batch_size=8,shuffle=True)
image,label = next(iter(train_dataloader))
print("image shape:",image.shape)
print("Label:",label)
 
# 查看数据的类别：10个类
classes_list = dataset.classes
print(classes_list)
plt.imshow(torchvision.utils.make_grid(image)[0],cmap="Greys")

输出：

image shape: torch.Size([8, 1, 28, 28])
Label: tensor([4, 3, 5, 8, 8, 1, 4, 3])
['T-shirt/top', 'Trouser', 'Pullover', 'Dress', 'Coat', 'Sandal', 'Shirt', 'Sneaker', 'Bag', 'Ankle boot']

5.3、扩散模型之退化过程

在扩散过程中需要为内容加入噪声，如何通过一个简单的方法控制图像的损坏程度，可以引入一个参数控制输入的“噪声量”。如下图绿色箭头方向（给图像中加入噪声）。

主要方法：将噪声张量与输入图像相乘，并将结果与1减去噪声的比例相加，以实现损坏效果

代码实现：

# 通过引入一个参数来控制输入的“噪声”来控制内容的损坏的程度
noise = torch.rand_like(image)
 
def corrupt(image,amount):
    """根据p为输入image添加噪声,这就是退化过程"""
    noise = torch.rand_like(image)
    amount = amount.view(-1,1,1,1)
    print(amount.shape)
    return image*(1-amount)+noise*amount
 
 
# 对输出结果进行可视化
fig,axs = plt.subplots(2,1,figsize=(10,5))
axs[0].set_title("Input data")
axs[0].imshow(torchvision.utils.make_grid(image)[0],cmap="Greys")
 
# 加入噪声:获得一个一维张量，包含从0-1之间的均匀分布的数值
p = torch.linspace(0,1,image.shape[0])
print(p.shape)
noised_image = corrupt(image,p)
print(noised_image.shape)
# 绘制加噪声的图像
axs[1].set_title("corrupt data")
axs[1].imshow(torchvision.utils.make_grid(noised_image)[0],cmap="Greys")

输出：

torch.Size([8])
torch.Size([8, 1, 1, 1])
torch.Size([8, 1, 28, 28])

5.3.1、其他噪声

图像中常见的噪声有以下四种：高斯噪声，泊松噪声，乘性噪声，椒盐噪声。

高斯噪声：是指其概率密度函数服从高斯分布（即正态分布）的一类噪声。如果一个噪声，它的幅度分布服从高斯分布，而它的功率谱密度又是均匀分布的，则称它为高斯白噪声。高斯白噪声的二阶矩不相关，一阶矩为常数，是指先后信号在时间上的相关性。产生原因：1）图像传感器在拍摄时视场不够明亮、亮度不够均匀；2）电路各元器件自身噪声和相互影响；3）图像传感器长期工作，温度过高。

泊松噪声：是指像素点出现的概率不是固定的，而是在某个范围内随机波动的一种噪声。产生原因：1）图像传感器在拍摄时视场不够明亮、亮度不够均匀；2）电路各元器件自身噪声和相互影响；3）图像传感器长期工作，温度过高。

乘性噪声：是指将两个或多个低频信号相乘后得到高频信号再进行处理时产生的伪像。产生原因：1）图像传感器在拍摄时视场不够明亮、亮度不够均匀；2）电路各元器件自身噪声和相互影响；3）图像传感器长期工作，温度过高。

椒盐噪声：是指由于图像切割引起的黑白相间的亮暗点噪声，与图像信号的关系是相乘。产生原因：1）图像传感器在拍摄时视场不够明亮、亮度不够均匀；2）电路各元器件自身噪声和相互影响；3）图像传感器长期工作，温度过高。

（1）高斯噪声

        高斯噪声是一种通过向输入数据添加均值为零和标准差 (σ)的正态分布随机值而产生的噪声。 正态分布，也称为高斯分布，是一种连续概率分布，由其概率密度函数 (PDF) 定义：

其中是随机变量， 是均值，$σ$" role="presentation" style="position: relative;">$σ$$\$\sigma \$$ 是标准差。高斯噪声在数字信号处理中被广泛应用，例如在通信、图像处理、语音识别等领域。 

# 读取一张图像并将图像和高斯噪声合并
import numpy as np
import matplotlib.pyplot as plt
import cv2
 
# 读取图像
image = cv2.imread('./image/2008_000536.jpg', cv2.IMREAD_GRAYSCALE)
 
# 设置参数
mu = 0
sigma = 1
size = image.shape
 
# 生成高斯噪声
noise = np.random.normal(mu, sigma, size)
 
# 将高斯噪声添加到图像上
noisy_image = image + noise
 
# 显示原始图像和带噪声的图像
plt.subplot(131), plt.imshow(image,cmap="gray"), plt.title('Original Image')
plt.subplot(132), plt.imshow(noise), plt.title('Noise')
plt.subplot(133), plt.imshow(noisy_image,cmap="gray"), plt.title('Noisy Image')
plt.show()

（2）椒盐噪声

        椒盐噪声是一种图像噪声，通常是由于图像切割引起的黑白相间的亮暗点噪声，与图像信号的关系是相乘。椒盐噪声的成因可能是影像讯号受到突如其来的强烈干扰而产生、类比数位转换器或位元传输错误等。例如失效的感应器导致像素值为最小值，饱和的感应器导致像素值为最大值 。

        在数字图像中，椒盐噪声是一种因为信号脉冲强度引起的噪声，可以用中值滤波来消除。给一副数字图像加上椒盐噪声的步骤如下：指定信噪比 SNR （其取值范围在 [0, 1]之间）计算总像素数目 SP，得到要加噪的像素数目 NP = SP * (1-SNR)随机获得要加噪的像素位置。

计算步骤如下：

1. 指定信噪比SNR（其取值范围在[0, 1]之间）。
2. 计算总像素数目SP。
3. 根据信噪比和总像素数目计算要加噪的像素数目NP = SP * (1-SNR)。
4. 随机获得要加噪的像素位置，即在图像中随机选择NP个像素点。
5. 对于每个被选中的像素点，将其值设置为一个随机数，这个随机数的范围在图像的最大值和最小值之间。
6. 将加噪后的图像保存或显示出来。

代码实现：

# 添加椒盐噪声
import random
# 在原图中增加椒盐噪声
def salt_and_pepper_noise(image, prob):
    output = np.zeros(image.shape, np.uint8)
    thres = 1 - prob 
    for i in range(image.shape[0]):
        for j in range(image.shape[1]):
            rdn = random.random()
            if rdn < prob:
                output[i][j] = 0
            elif rdn > thres:
                output[i][j] = 255
            else:
                output[i][j] = image[i][j]
    return output
 
# 生产指定大小的椒盐噪声
def generate_salt_pepper_noise(size, salt_ratio=0.05, pepper_ratio=0.05):
    # 创建一个大小为size的空数组
    noise = np.zeros((size, size))
 
    # 遍历数组中的每个元素
    for i in range(size):
        for j in range(size):
            # 计算当前元素的盐和胡椒概率
            salt_prob = salt_ratio
            pepper_prob = pepper_ratio
 
            # 如果当前元素位于边缘，则增加胡椒概率
            if i == 0 or i == size - 1 or j == 0 or j == size - 1:
                pepper_prob += 0.05
 
            # 随机选择盐或胡椒
            if np.random.rand() < salt_prob:
                noise[i, j] = 1
            elif np.random.rand() < pepper_prob:
                noise[i, j] = 2
 
    return noise
 
 
# 读取图像
image = cv2.imread('./image/2008_000536.jpg',cv2.COLOR_BGR2GRAY)
pepper_image = generate_salt_pepper_noise(28,0.1)
salt_pepper = salt_and_pepper_noise(image,0.2)
plt.subplot(131),plt.imshow(image),plt.title("Original Image")
plt.subplot(132),plt.imshow(salt_pepper),plt.title("Salt Image")
plt.subplot(133),plt.imshow(pepper_image),plt.title("Pepper Image")

5.4、扩散模型之网络模型

训练模型之前，需要一个能够输入28×28像素的噪声图像，然后输出相同大小图片的预测结果。本文采用经典的Unet网络，由一条“压缩路径”和一条“扩展路径”组成。“压缩路径”主要是通过压缩数据的维度，而“扩展路径”则将数据扩展回原始维度（类似于自动编码器）。

1、自动编码器：

是一种无监督学习模型，主要用于数据压缩和特征提取。它由两部分组成：编码器和解码器。

编码器将输入数据压缩成一个低维表示，解码器则将这个低维表示还原成原始数据。在这个过程中，自动编码器学习如何有效地将原始数据转换为压缩表示，以及如何从压缩表示中恢复原始数据。这种学习过程使得自动编码器能够捕捉到数据中的有用信息，同时去除冗余和噪声。

自动编码器的基本思想是通过训练一个神经网络来最小化输入数据和解码器输出之间的差异（即重构误差）。在训练过程中，网络的权重会不断地更新，以便更好地捕捉数据的特征和结构。

自动编码器在许多领域都有广泛的应用，如图像识别、语音识别、自然语言处理、推荐系统等。它们可以用于降维、特征提取、异常检测、生成模型等任务。

这段代码定义了一个名为BasicNet的类，用于实现Unet网络。Unet网络是一种用于图像分割的卷积神经网络，它包括一个编码器（下采样）和一个解码器（上采样）。这个简单的Unet网络包含三个卷积层和相应的激活函数、池化层和上采样层。

1. `__init__`方法：定义了网络的结构。
   - `self.down_layers`：定义了下采样层，包括三个卷积层，每个卷积层后面跟着一个ReLU激活函数和一个最大池化层。
   - `self.up_layers`：定义了上采样层，包括三个卷积层，每个卷积层后面跟着一个ReLU激活函数。
   - `self.act`：定义了激活函数，这里使用ReLU激活函数。
   - `self.downscale`：定义了最大池化层，用于降低特征图的空间尺寸。
   - `self.upscale`：定义了上采样层，用于增加特征图的空间尺寸。

2. `forward`方法：实现了网络的前向传播过程。
   - 对于下采样层，首先将输入x通过第一个卷积层、激活函数和池化层，然后将结果添加到h列表中，并将x通过最大池化层降维。重复这个过程，直到处理完所有的下采样层。
   - 对于上采样层，首先将输入x通过上采样层和激活函数，然后将结果与h列表中的最后一个元素相加，再通过第二个卷积层、激活函数和池化层。重复这个过程，直到处理完所有的上采样层。
   - 最后返回x作为网络的输出。

# Unet网络
class BasicNet(nn.Module):
    """一个十分简单的Unet网络部署"""
    def __init__(self,in_channels=1,out_channels=1):
        super().__init__()
        self.down_layers = torch.nn.ModuleList(
            [
                nn.Conv2d(in_channels,32,kernel_size=5,padding=2),
                nn.Conv2d(32,64,kernel_size=5,padding=2),
                nn.Conv2d(64,64,kernel_size=5,padding=2),
            ]
        )
        self.up_layers = torch.nn.ModuleList(
            [
                nn.Conv2d(64,64,kernel_size=5,padding=2),
                nn.Conv2d(64,32,kernel_size=5,padding=2),
                nn.Conv2d(32,out_channels,kernel_size=5,padding=2),
            ]
        )
        self.act = nn.ReLU()
        self.downscale = nn.MaxPool2d(2)
        self.upscale = nn.Upsample(scale_factor=2)
 
    def forward(self,x):
        h = []
        for i,l in enumerate(self.down_layers):
            x = self.act(l(x))
            if i < 2:
                h.append(x)
                x = self.downscale(x)
        
        for i,l in enumerate(self.up_layers):
            if i > 0:
                x = self.upscale(x)
                x += h.pop()
            x = self.act(l(x))
        return x
    
net = BasicNet()
x = torch.rand(8,1,28,28)
print(net(x).shape)
# 模型参数
sum([param.numel() for param in net.parameters()])
print(net)