直方图处理

• 令 $r_k(k = 0,1,2...L-1)$ 表示一幅$L$级灰度数字图像 $f(x,y)$ 的灰度. $f$ 的非归一化直方图定义为: $h(r_{k}) = n_{k}, k = 0,1,2...L-1$
  • 式中,$n_k$是$f$中灰度为$r_k$的像素的数量,并且细分的灰度级称为直方图容器
• 归一化直方图定义为 $p(r_{k}) = \frac{h(r_{k})}{MN} = \frac{n_{k}}{MN}$
  • 式中, $M$和$N$分别是图像的行数和列数.对$k$的所有值, $p(r_{k})$的和总是1
• 像素占据整个灰度级范围并且均匀分布的图像, 将具有高对比度的外观和多种灰色调

直方图均衡化

• 由像素累计概率确定对应的新像素值
• 直方图均衡化的目的: 为了生成一幅具有均匀直方图的输出图像
• 变换函数: $s_k=T\left(r_k\right)=(L-1) \sum_{j=0}^k p_r\left(r_j\right), \quad k=0,1,2, \cdots, L-1$
• 直方图均衡化改变的只是对比度并非内容

直方图匹配(规范化)

• 寻找累计概率密度最相近的像素值作为新像素值
• 定义: 用于生成具有规定直方图的图像的方法

局部直方图处理

• 目的: 解决增强图像中几个小区域的细节
• 解决方法: 设计基于像素邻域的灰度分布的变换函数
• 局部直方图的处理过程:
  1. 定义一个邻域, 并将其中心在水平方向或垂直方向上从一个像素移动到另一个像素
  2. 在每个位置, 计算邻域中的各点的直方图, 得到直方图均衡化或直方图规定化变换函数, 将这个函数映射于邻域中心点像素的灰度
  3. 然后将邻域的中心移到一个相邻像素位置, 并重复上述过程

图像均值和方差的含义

• 均值是平均灰度的测度
• 方差(或标准差)是图像对比度的测度
• 简单理解：均值越大, 图像越亮; 方差越大, 图像对比度越高

MTF的理解

• MTF为调制传递函数, 为分辨率和对比度两者间的关系
  • 分辨率单位为每毫米线对(lp/mm), 对比度定义为最小与最大强度值从物平面传输到像平面的程度
  • MTF是衡量其在特定分辨率下将对比度从物体转移到图像的能力, 随着测试目标上的线间距减小(频率增加), 镜头越来越难以有效地传递对比度,结果MTF下降
  • MTF会随着空间分辨率地增加而降低
• 刃边法SFR求MTF
  • 对黑白斜边超采样得到ESF(边缘扩散函数)
  • 对ESF求导得到直线的变化率LSF(线扩散函数)
  • 对LSF进行FFT变换得到各个频率下MTF

PSF,LSF,ESF

• PSF(点扩散函数): 为点光源成像后的亮度分布函数, 点扩散函数是中心圆对称
• ESF(边缘扩散函数): 为一条由白到黑的线
• LSF(线扩散函数)
• 当获取点光源的亮度分布函数PSF(x,y)后, 对其进行二维傅里叶变换即可得到MTF
  • 但由于实际中点光源强度很弱, 一般较少采用
  • 相对于PSF而言, LSF能量得到一定加强, 因此用LSF更好

膨胀, 腐蚀, 开运算, 闭运算

• 膨胀
  • 利用卷积核B对图像A进行卷积操作, 卷积核可以是任意形状或大小
  • 卷积核B通常有一个自定义的参考点, 一般位于中心位置
  • 膨胀是求局部最大值的操作: 当卷积核B扫描图像A与其进行卷积操作时, 计算模板B覆盖的区域的最大值并将最大值赋给模板的参考点. 因为图像中亮点的灰度值大, 所以膨胀操作会使得图像中的高亮区域逐渐增长
• 腐蚀
  • 膨胀的反操作, 腐蚀计算的是局部区域的最小值
  • 将卷积核B与图像A进行卷积, 将B所覆盖区域的最小值赋给参考点
  • 腐蚀操作会使得图像中亮的区域变小, 暗的区域变大
• 开运算: 先腐蚀再膨胀
  • 开运算可以移除较小的明亮区域, 在较细的地方分离物体(假设小物体是亮色, 关注的前景是黑色, 可以使用开运算移出亮点)
• 闭操作: 先膨胀再腐蚀
  • 闭运算可以填充五体内的细小空洞, 连接邻近的明亮物体
• 顶帽操作: 原图与开操作的差
  • 局部亮度极大点被分割出来(可以分两步理解, 开操作移除了明亮的小区域, 当用原图减去开操作的结果之后, 之前被移除的明亮区域就会凸显出来, 因此看到的效果就是一些亮度较大的小区域)
• 黑帽操作: 闭操作与原图的差
  • 局部黑色的洞被分割出来

边缘检测算子及其对比

• 绝大多数的边缘检测方法可以分为两类: 基于查找的方法和基于零穿越的方法
  • 基于查找的方法: 通过寻找图像一阶导数中的最大和最小值来检测边缘, 通常是将边界定位在梯度最大的方向
  • 基于零穿越的方法: 通过寻找图像二阶导数零穿越来寻找边界, Laplacian过零点或非线性差分表示的过零点
  • Canny算子
• Roberts算子
  • 利用局部交叉差分寻找边缘的算子, 常用来处理具有陡峭的低噪声图像, 当图像边缘接近于正45度或负45度时, 该算法处理效果更理想
  • 优点: 从图像处理的实际效果看, 边缘定位较准, 对噪声敏感, 适用于边缘明显且噪声较小的图像分割
  • 缺点: 提取的边缘线条较粗
  • $g_x = \frac{\partial f}{\partial x}=P9-P5 \\ g_y = \frac{\partial f}{\partial y} =P8-P6$
  • 在代码实现方面, 便可以如同卷积一样构造两个滤波器矩阵对图像进行卷积, 假设使用第一个模板卷积后得到的结果为 $f_y$ , 使用第二个模板卷积后得到的结果为 $f_x$ , 那么最终的结果为这两个中间结果的加权平均, 例如 $0.5*f_x +0.5*f_y$
• Prewitt算子
  • 采用3*3模板对区域内的像素值进行计算, 在水平方向和垂直方向均比Roberts算子更加明显, 适合用来识别噪声较多灰度渐变的图像
  • 优点: Prewitt算子对噪声有抑制作用, 抑制噪声的原理是通过像素平均
  • 缺点: 该算子具有平滑的作用, 但是像素平均相当于对图像的低通滤波, 所以Prewitt算子对边缘的定位不如Roberts算子
  • $g_x = \frac{\partial f}{\partial x}=(P7+P8+P9)-(P1+P2+P3) \\ g_y = \frac{\partial f}{\partial y}=(P3+P6+P9)-(P1+P4+P7)$
• Sobel算子
  • Sobel算子在Prewitt算子的基础上增加了权重的概念, 认为相邻点的距离远近对当前像素点的影响是不同的, 距离越近的像素点对应当前像素的影响越大, 从而实现图像锐化并突出边缘轮廓
  • Sobel算子结合了高斯平滑和微分求导, 因此结果会具有更多的抗噪性, 当对精度要求不是很高时, Sobel算子是一种较为常用的边缘检测方法
  • 优点：由于该算子中引入了类似局部平均的运算, 因此对噪声具有平滑作用, 能很好的消除噪声的影响, 边缘定位效果不错. Sobel算子对于象素的位置的影响做了加权, 与Prewitt算子、Roberts算子相比因此效果更好
  • 缺点：但检测出的边缘容易出现多像素宽度
  • $g_x = \frac{\partial f}{\partial x}=(P7+2*P8+P9)-(P1+2*P2+P3) \\ g_y = \frac{\partial f}{\partial y}=(P3+2*P6+P9)-(P1+2*P4+P7)$
• Laplacian算子
  • 各向同性算子, 不能检测出边的方向, 对孤立像素的响应要比对边缘或线的相应更强烈, 只适用于无噪声图像; 存在噪声的情况下, 使用Laplace算子检测之前需低通滤波
  • Laplacian算子一般不用于边的检测, 常用来判断边缘像素位于图像的明区或暗区
  • 定义: $Laplacian(f)=\frac{\partial^2f}{\partial x^2}+\frac{\partial^2f}{\partial y^2}$
    • $\frac{\partial^2f}{\partial x^2} = \frac{f(x+1)-f(x)}{\partial x} = [f(x+1)-f(x)]-[f(x)-f(x-1)] \\ \frac{\partial^2f}{\partial y^2} = \frac{f(y+1)-f(y)}{\partial y} = [f(y+1)-f(y)]-[f(y)-f(y-1)]$
    • $\nabla^2(f)=\frac{\partial^2f}{\partial x^2}+\frac{\partial^2f}{\partial y^2} =[f(x+1,y)+f(x-1,y)+f(x,y+1)+f(x,y-1)]-4f(x,y)$
  • 该算子对孤立点或端点更为敏感, 因此特别适用于以突出图像中的孤立点、孤立线或线端点为目的的场合. 同梯度算子一样, 拉普拉斯算子也会增强图像中的噪声, 有时用拉普拉斯算子进行边缘检测时, 可将图像先进行平滑处理
  • Laplacian算子进行边缘检测并没有像Sobel或Prewitt那样的平滑过程, 所以它会对噪声产生较大的响应, 并且无法分别得到水平方向、垂直方向或者其他固定方向的的边缘. 但是它只有一个卷积核, 所以计算成本会更低
  • Laplacian锐化
    • 由于拉普拉斯是一种微分算子, 它的应用可增强图像中灰度突变的区域, 减弱灰度的缓慢变化区域. 因此, 锐化处理可选择拉普拉斯算子对原图像进行处理, 产生描述灰度突变的图像, 再将拉普拉斯图像与原始图像叠加而产生锐化图像
    • $g(x, y)= \begin{cases}f(x, y)-\nabla^2 f(x, y) & \text { 如果拉普拉斯掩模中心系数为负 } \\ f(x, y)+\nabla^2 f(x, y) & \text { 如果拉普拉斯掩模中心系数为正 }\end{cases}$
• Canny算子
  • 是传统图像里被用的最多的边缘检测算子
  • Canny提出了一个对于边缘检测算法的评价标准, 包括：
    1. 以低的错误率检测边缘, 也即意味着需要尽可能准确的捕获图像中尽可能多的边缘
    2. 检测到的边缘应精确定位在真实边缘的中心
    3. 图像中给定的边缘应只被标记一次, 并且在可能的情况下, 图像的噪声不应产生假的边缘
    • 检测算法要做到：**边缘要全, 位置要准, 抵抗噪声的能力要强. **
  • 该算子求边缘点的具体算法步骤如下
    1. 用高斯滤波器平滑图像：边缘检测算子受噪声的影响都很大. 那么, 我们第一步就是想到要先去除噪声, 因为噪声就是灰度变化很大的地方, 所以容易被识别为伪边缘
    2. 用一阶偏导有限差分计算梯度幅值和方向, 例如 Sobel
    3. 对梯度幅值进行非极大值抑制：sobel算子检测出来的边缘太粗了, 我们需要抑制那些梯度不够大的像素点, 只保留最大的梯度, 从而达到瘦边的目的. 通常灰度变化的地方都比较集中, 将局部范围内的梯度方向上, 灰度变化最大的保留下来, 其它的不保留, 这样可以剔除掉一大部分的点. 将有多个像素宽的边缘变成一个单像素宽的边缘. 即“胖边缘”变成“瘦边缘”
    4. 用双阈值算法检测和连接边缘：通过非极大值抑制后, 仍然有很多的可能边缘点, 进一步的设置一个双阈值, 即低阈值（low）, 高阈值（high）. 灰度变化大于high的, 设置为强边缘像素, 低于low的, 剔除. 在low和high之间的设置为弱边缘. 对每一个弱边缘进一步判断, 如果其领域内有强边缘像素, 保留, 如果没有, 剔除

2 空间滤波基础

• 基础知识的注意点:
  • 低通滤波器可以消除噪声，高通滤波器可以提取边缘
  • 为什么0表示黑色，255表示白色：可能是模拟人眼对光线的机理，黑色表示没有光线，白色相反

图像信息的高频信息

• 高频信息：图像中的边缘和其他急剧的灰度变化(如噪声)

基础滤波算法

• 平滑(低通)空间滤波器: 通过模糊图像来平滑图像, 使用积分运算实现
  • 用于降低灰度的急剧过度, 因为随机噪声通常是由灰度的急剧过渡组成
  • 均值滤波：实质是归一化之后地方框滤波, 不能很好地保护图像细节, 不能很好地去除噪声点, 但对周期性地干扰噪声有很好的抑制作用
  • 低通高斯滤波核: 高斯核是唯一可分离的圆对称核, 对于抑制服从正态分布的噪声非常有效, 可以产生更加均匀的平滑结果
• 锐化(高通)空间滤波器:突出灰度中的过渡, 使用微分运算实现
  • 拉普拉斯算子: $\Delta f= \frac{\partial ^{2}f}{\partial x^{2}}+ \frac{\partial ^{2}f}{\partial y^{2}}$
  • $\Delta f = f(x+1, y) + f(x-1, y) + f(x, y+1) + f(x, y-1)-4f(x,y)$
  • 使用拉普拉斯算子锐化图像的基本方法: $g(x, y)=f(x, y)+c\left[\nabla^2 f(x, y)\right]$
• 钝化遮掩: 从原图减去钝化(平滑后)图像, 可以用来锐化图像
  • 步骤：
    1. 模糊原图像
    2. 从原图像减去模糊后的图像(产生的差称为模板)
    3. 将模板与原图像相加
• 使用罗伯特交叉梯度算子和Sobel算子
  • 使用一阶梯度算子锐化图像咳用于工业缺陷检测
  • 使用拉普拉斯算子来突出图像细节, 使用梯度来增强突出边缘

双边滤波与导向滤波

• 双边滤波和引导(导向)滤波都是保边滤波器, 能在滤波过程中有效的保留图像中的边缘信息
• 双边滤波器
  • 同时考虑了像素空间的差异与强度差异的滤波器, 具有保持图像边缘的特征
  • 与高斯滤波器相比, 在原有高斯函数基础上增加一项使权重与图像边缘梯度相关, 不仅仅与空间距离相关, 且与灰度距离相关
  • $\begin{aligned} B F & =\frac{1}{W_q} \sum_{p \in S} G_s(p) G_r(p) * I_p \\ & =\frac{1}{W_q} \sum_{p \in S} \exp \left(-\frac{\|p-q\|^2}{2 \sigma_s^2}\right) \exp \left(-\frac{\left\|I_p-I_q\right\|^2}{2 \sigma_r^2}\right) * I_p \end{aligned}$
• 引导(导向)滤波器
  • 引导滤波的思想用一张引导图像产生权重, 从而对输入图像进行处理
  • 保边滤波后的结果和引导图像在滤波窗口内呈现线性关系: $q_i=a_k I_i+b_k,\forall i \in \omega_k$, $I$是引导图像, $q$是输出, $\omega$是以$k$为中心像素的窗口, $a$和$b$都为该窗口对应的线性系数
  • 期望$a$和$b$随着图像内容变化, 在边缘区域,$a$尽量大,保持梯度;在平滑区域$a$尽量小,以尽量平滑. 利用带有正则项的岭回归模型求解: $E(a_k, b_k)=\sum_{i \in \omega_k}((a_k I_i+b_k-p_i)^2+\epsilon a_k^2)$
    • 求解上述方程,可以得到: $a_k=\frac{\frac{1}{|\omega|} \sum_{i \in \omega_k} I_i p_i-\mu_k \bar{p}_k}{\sigma_k^2+\epsilon}$, $b_k=\bar{p}_k-a_k \mu_k$, $q_i=\frac{1}{|\omega|} \sum_{k \mid i \in \omega_k}(a_k I_i+b_k)$
    • 在公式中, $\mu_k$ 与 $\sigma_k$ 表示 $I_i$ 在窗口内的均值、标准差, $|w|$ 表示窗口内像素块的总数, $\overline{p_k}$ 表示窗口内输入图像 $p$ 的均值
  • 最终算法:
    1. $\begin{aligned} & \operatorname{mean}_I=f_{\text {mean }}(I) \\ & \operatorname{mean}_p=f_{\text {mean }}(p) \\ & \operatorname{corr}_I=f_{\text {mean }}(I . * I) \\ & \operatorname{corr}_{I p}=f_{\text {mean }}(I . * p)\end{aligned}$
    2. $\begin{aligned} & \operatorname{var}_I=\operatorname{corr}_I-\operatorname{mean}_I \cdot \operatorname{mean}_I \\ & \operatorname{cov}_{I p}=\operatorname{corr}_{I p}-\operatorname{mean}_I \cdot \operatorname{mean}_p\end{aligned}$
    3. $\begin{aligned} & a=\operatorname{cov}_{I p} \cdot /\left(\operatorname{var}_I+\epsilon\right) \\ & b=\operatorname{mean}_p-a . * \operatorname{mean}_I\end{aligned}$
    4. $\begin{aligned} \operatorname{mean}_a & =f_{\text {mean }}(a) \\ \operatorname{mean}_b & =f_{\text {mean }}(b)\end{aligned}$
    5. $q=\operatorname{mean}_a * I+\operatorname{mean}_b$
    • $f_{mean}$为窗口半径为$r$的均值滤波

图像的抗混叠滤波

因为无法对图像无限地取样，因此在数字图像中总会出现混叠

• 混叠：空间混叠和时间混叠
  • 空间混叠：由欠采样引起的，主要问题是引入伪影，即原始图像中未出现的线条锯齿、虚假高光和频率模式
  • 时间混叠：与动态图像序列中图像的时间间隔相关，如车轮效应
• 在重取样之前使用低通滤波器进行平滑，以衰减图像的高频分量，可以有效地抑制混叠，但图像清晰度下降

常见的三种图像插值方法

1. 最近邻插值: 在待求像素的四邻域中, 将距离待求像素最近的邻像素的灰度值赋给待求像素
2. 双线性插值: 利用待求像素四个相邻像素的灰度在两个方向上做线性插值
3. 三次内插法: 利用三次多项式求逼近理论上的最佳插值函数, 待求像素(x,y)的灰度值由其周围16个灰度值加权内插得到

图像降噪算法专题

空间域滤波降噪

该方法主要针对随机噪声
空间域滤波降噪的方法主要是通过分析在一定窗口大小内中心像素和其他相邻像素之间的直接联系，来获取新的中心像素值的方法

算术均值滤波

• 算术均值滤波：用像素邻域的平均值代替像素值，适用于脉冲噪声（灰度级一般与周围像素的灰度级不相关，且亮度高出其他像素值许多）；频域角度考虑，均值滤波器为低通滤波器，可以帮助消除图像尖锐噪声实现图像平滑模糊
  • 方法：取待处理像素的L邻域计算平均灰度以替代其像素值；均值滤波效果随L的增大而变得模糊，图像对比度随之降低
  • 优点：适用于具有一般随机干扰的信号（信号在平均值附近波动）
  • 缺点：耗费RAM资源，计算速度较慢

1

void  cv::blur(InputArray src, OutputArray dst, Size ksize, Point anchor = Point(-1,-1), int borderType = BORDER_DEFAULT)		

• 为提高均值滤波器的效果：
  • 增加对边缘方向的判断，以提高锐度，但做均值的像素减少去噪效果下降；
  • 提出NLM(none local mean)算法进行降噪

中值滤波

• 中值滤波器：使用像素邻域内的中值代替像素值。由于中值滤波不会处理最大值和最小值，故不会受到噪声的影响。也可用于处理边缘信息；但中值滤波器会清除掉某些区域的纹理信息
  • 方法：取待处理像素的L邻域计算中值以替代其像素值；对椒盐噪声的消除效果好于均值滤波
  • 优点：能克服偶然因素引起的波动干扰
  • 缺点：当噪声像素个数大于窗口像素总数的一半时，由于中值仍为噪声像素灰度值，滤波效果较差；若加大窗口尺寸，会使得原边缘像素被其他区域像素代替使图像变得模糊加大运算量
• 由于窗口大小固定，中值滤波不能兼顾去噪和保护图像细节
  • 自适应中值滤波(Adaptive Median Filter)，会自动调节窗口大小，并判断像素是否为噪声需要处理

1

void medianBlur(InputArray src, OutputArray dst, int ksize);

高斯滤波器

• 高斯滤波与均值滤波的区别：
  • 均值滤波中，滤波器的各个像素的权重相同，滤波器为线性滤波器；高斯滤波器中的像素的权重与其距中心的距离成比例，高斯滤波器为非线性滤波器
  • 高斯滤波考虑了空间距离对中心像素的影响, 添加了空间权重, 因此边缘较均值滤波更锐利
• 高斯滤波：高斯滤波矩阵的权值，随着与中心像素点的距离增加，而呈现高斯衰减的变换特性；离算子中心很远的像素点的作用很小，从而能在一定程度上保持图像的边缘特征
  • 方法：与均值滤波相似，生成L邻域的高斯滤波算子，计算L邻域内的像素与滤波算子的卷积结果代替像素值
    • $G B[I]_{\mathrm{p}}=\sum_{\mathbf{q} \in \mathcal{S}} G_\sigma(\|\mathbf{p}-\mathbf{q}\|) I_{\mathbf{q}}$
    • $G_\sigma(x)=\frac{1}{2 \pi \sigma^2} \exp \left(-\frac{x^2}{2 \sigma^2}\right)$
    • 其中$p$为中心像素位置, $q$为另一个像素位置, $I_q$为$q$像素像素值, 两者之间的距离使用$G_\sigma(||p-q||)$获取, $\sigma$决定了高斯分布的宽度, 大小决定了以$p$为中心涵盖了多少临近像素用于参与计算
  • 对高斯噪声的滤除效果较好

1

void cv::GaussianBlur(InputArray 	src, OutputArray 	dst, Size 	ksize, double 	sigmaX, double 	sigmaY, int 	borderType = BORDER_DEFAULT)

双边滤波

• 上述高斯滤波仅仅考虑空间位置关系, 并没有考虑图像本身像素值的变化情况(图像边缘像素值变化剧烈, 非边缘区域变化平坦), 因此需要能够衡量像素变化剧烈程度的量
• 双边滤波：非线性滤波，结合图像的空间邻近度和像素相似度的折中处理，同时考虑空间信息和灰度相似度(可解决图像边缘问题)，达到保边去噪的目的；
  • 方法：使用两个滤波器（一个函数由几何空间距离决定滤波器系数，另一个由像素差值决定滤波器系数）
    • $BF[I]_{\mathrm{p}}=\frac{1}{W_{\mathbf{p}}} \sum_{\mathbf{q} \in \mathcal{S}} G_{\sigma_s}(\|\mathbf{p}-\mathbf{q}\|) G_{\sigma_r}\left(I_{\mathbf{p}}-I_{\mathbf{q}}\right) I_{\mathbf{q}}$
    • $W_{\mathrm{p}}=\sum_{\mathrm{q} \in \mathcal{S}} G_{\sigma_s}(\|\mathrm{p}-\mathrm{q}\|) G_{\sigma_{\mathrm{r}}}\left(I_{\mathrm{p}}-I_{\mathrm{q}}\right)$
    • 简单来说双边滤波模板组合由两个模板生成：第一个为高斯模板(离中心点远的对中心点的贡献小)，第二个为以像素值的差值作为函数系数生成的模板(像素差值越大对中心点的贡献小)，两模板点乘即可得到最终滤波模板; 第一个高斯模板只要是进行滤波操作, 而第二个以像素值差值的模板主要是进行保边操作
• 参数讨论
  • 单单考虑值域上的滤波不考虑空间上高斯滤波, 值域滤波可以看作直方图变换, 使图像的值域范围朝向峰值方向压缩
  • $\sigma_s$与$\sigma_r$, 值越大对双边滤波器的贡献就越小, 相反值越小则贡献越大, 得到强调;
    • $\sigma_s$变小, 则更多的是利用邻域的值进行平滑, 说明空间上的信息更加重要;
    • $\sigma_r$变小, 则更多的是利用相似的值进行保边, 说明值域上的信息更加重要;
  • $\sigma_s$表示的是空间上的平滑, 对没有边缘的更加合适;而$\sigma_r$强调值域的差别, 对具有边缘的区域更加重要, 减小该值以更好的保边
    • $\sigma_s$变大, 图像每个区域的权重更多源于值域滤波的权重, 对空间邻域信息不敏感
    • $\sigma_r$变大, 图像每个区域的权重更多源于空间距离的权重, 相似于高斯滤波, 对空间邻域信息不敏感, 保边性能下降
    • 需要去除平滑区域的噪声, 应提高$\sigma_s$, 而需要保持边缘, 应减小$\sigma_r$
• 可以很好的保留图像边缘细节而滤除掉低频分量的噪声，但效率较低，花费时间较长

1

void bilateralFilter(InputArray src, OutputArray dst, int d, double sigmaColor, double sigmaSpace, int borderType=BORDER_DEFAULT )

引导滤波

• 引导滤波(导向滤波)是一种图像滤波技术，通过一张引导图对初始图像(输入图像)进行滤波处理，使得最后的输出图像大体上与初始图像相似，但是纹理部分与引导图相似
• 原理:
  • 假设引导图像$G$, 输出图像为$O$, 则对应图像的任何一个位置$k$, 滤波窗口为$w_k$, 则$O_i = a_k * G_i + b_k, i \in w_k$; 该方式就可以保证滤波输出与引导图尽量一致
  • 因此为使输入图像和输出图像在局部内容大致相同, 最优化目标为$min \sum_{i \in w_k} (O_i - I_i)^2 = min \sum_{i \in w_k}(a_k*G_i+b_k-I_i)^2 + \epsilon * a_k ^2$, $\epsilon * a_k ^2$为控制变量, 防止除0
  • 上述公式进行求导并使其等于0以求极值, 可计算得:
    • $b_k=mean(I)-a_k*mean(G)$
    • $a_k=\frac{cov(G, I)}{var(G)+\epsilon}$
• 去噪实验
  • 去噪时$O = \sum (mean(I) + a_k * (G - mean(G)))$, 考虑窗口内的I的均值滤波(低频分量)与引导图的高频细节的加权, 最终表现为I的平坦区域与G的梯度剧烈变化区域(I的低频分量与G的高频细节的结合)
  • 平坦区域, $var(I) \rightarrow 0, a_k \rightarrow 0, b_k \rightarrow mean(I), O \rightarrow mean(I)$更多的是均值滤波, 反之输出引导图本身, 梯度得到保留
  • 参数的影响
    • 滤波核大小: 直接涉及均值滤波核大小, 核越大周围像素影响越大, 平滑越明显
    • $\epsilon$的影响: 控制$a_k$, $\epsilon$越大, $a_k$越小, 则更接近平滑滤波, 高频分量保留更少, 保边效果差, 降噪效果更好; 否则高频分量保留更多, 保边效果更好
      • 若$\epsilon=0$, 则滤波器没有作用, 将输入原封输出
      • 若$\epsilon>0$, 在像素强度变化小的区域, 相当于加权平均滤波; 而变化大的区域, 滤波效果很弱, 有助于保持边缘
        • $\epsilon$的作用用于界定像素强度变化的大小, $\epsilon$越大滤波效果越明显
• 实际实现中, 对整幅图像同时求解(窗口内的值求和, 并不进行取平均操作), 随后进行一次均值滤波
• 典型应用：保边、图像平滑、抠图、图像增强、HDR压缩、图像抠图及图像去雾等场景
• 重要应用：去雾算法
  • 随后提出快速导向滤波算法：通过下采样减少像素点后通过上采样恢复到原有的尺寸大小

1

void cv::ximgproc::guidedFilter(InputArray 	guide, InputArray 	src, OutputArray 	dst, int 	radius, double 	eps, int 	dDepth = -1)

非局部均值去噪(NLM)

• NLM算法使用自然图像中普遍存在的冗余信息来去噪，利用整幅图像去噪，以图像块为单位在图像中寻找相似区域，再对区域求平均，能够较好的去掉图像中存在的高斯噪声
  • 与高斯滤波的差异: 使用相似度计算权重, 相似度越大贡献越大
  • 基本思想：当前像素的估计值由图像中与它具有相似邻域结构的像素加权平均得到
    • TV算法在平滑噪声的同时也把很多图像本身的纹理边缘细节去掉
    • 各向异性滤波AD算法在保持细节信息的同时也保留很多的噪声
    • NL算法在去除噪声和保持纹理细节方面都具有较好的效果
  • 算法执行过程：
    • 对于图像中每个像素点, 建立两个窗口, 大窗口的为搜索窗, 小窗口为邻域窗口;其中邻域窗口分别建立在该像素点和搜索窗中. 通过滑动搜索窗中的邻域窗口计算两个邻域窗口的相似度进行赋值
    • 权值计算
      • $\hat{u}(x_0, y_0) = \frac{1}{C(x_0, y_0)} \sum_{(x, y) \in B(x_0, y_0, r)} u(x, y) w((x_0, y_0),(x,y)), C(x_0, y_0) = \sum_{(x, y) \in B(x_0, y_0, r)} w((x_0, y_0),(x,y))$
      • $w((x_0, y_0),(x,y)) = exp {(-\frac{max(d^2-2\epsilon^2, 0.0)}{h^2})}$
      • $d = \frac{1}{(2f+1)^2} \sum_{(x,y) \in (x_0, y_0, f)}(u(x,y)-u(x_0,y_0))^2$
      • 其中$r$为搜索窗的大小, f为邻域大小; 即每个点的值为搜索窗每个点的加权和. 而该权重的计算是利用两个邻域对应像素点值的差值的平方(对应像素点的像素值欧式距离), 并将该值赋予高斯权重以形成权重
• 参数分析
  • 当两块区域很相似度, 其差值就越小, 则权重就越大;而相似度很小时, 差值较大, 权重较小. 完全考虑了图像像素值变化对滤波的影响
  • h控制保边降噪效果：h越大时, 高斯曲线越平缓, 去噪效果越好, 但保边效果下降图像模糊;
  • $\epsilon$表示噪声的标准差, 小于$2\epsilon$的块设置为1, 而d较大的块权重较小, 贡献较少. 该参数主要用于将相似的块平均到噪声水平
• 考察不同大小的相似窗和搜索窗，在不同的噪声程度下，最后取得的去噪效果的PSNR值
  • 当相似窗大小$(2f+1)^2$中的f取为4,即$9 \times 9$时，对不同的噪声程度，最后的去噪效果都有较大的PSNR值
  • 当搜索窗大小$(2s+1)^2$中的s取为6,即搜索窗大小为$13 \times 13$时，对不同的噪声程度，最后的去噪效果都有最大的PSNR值
• 指标
  • MSE
  • PSNR, 峰值信噪比,即峰值信号的能量与噪声的平均能量之比，通常表示的时候取对数变成分贝dB
    • 由于MSE为真实图像与含噪图像之差的能量均值，两者的差为噪声，因此PSNR即峰值信号能量与MSE之比，定义为$P S N R=10 \log _{10} \frac{\text { MaxValue }^2}{M S E}=10 \log _{10} \frac{2^{\text {bits }}-1}{M S E}$

1
2
3
4

void fastNlMeansDenoising( InputArray src, OutputArray dst, float h = 3, int templateWindowSize = 7, int searchWindowSize = 21);
void fastNIMeansDenoisingColored( InputArray src, OutputArray dst, float h = 3, float hColor = 3, int templateWindowSize = 7, int searchWindowSize = 21);
void fastNlMeansDenoisingMulti( InputArray Of Arrays srcImgs, OutputArray dst,int imgToDenoiseIndex, int temporalWindowSize,  float h = 3, int templateWindowSize = 7, int searchWindowSize = 21);
void fastNlMeansDenoisingColoredMulti( InputArray Of Arrays srcImgs, OutputArray dst, int imgToDenoiseIndex, int temporalWindowSize, float h = 3, float hColor = 3,int templateWindowSize = 7, int searchWindowSize = 21);

BM3D去噪

维纳滤波(最小均方误差滤波)

• 该方法建立在图像和噪声都是随机变量的基础上, 为了寻找为未污染图像$f$与一个估计$\hat{f}$, 使其均方误差最小, $e^2 = E{((x-\hat{f})^2)}$
• 频域推导过程
  • 原始信号为$X(\omega)$, 噪声信号为$V(\omega)$, 观测到的信号$Y(\omega)=X(\omega)+V(\omega)$$
  • 假定维纳滤波器为$H(\omega)$, 则滤波输出$Z(\omega) = H(\omega) Y(\omega)$
  • 频域的误差为$e(\omega) = X(\omega)-Z(\omega)$， 代价函数为$J(H) = E^2(e(\omega))=E^2(X)-H P_{yx} - H^{*}P_{xy}+H^2P_{yy}$, 其中P表示功率谱和互功率谱
  • 对上式求导可得维纳滤波器的最优解为$H(\omega) = \frac{P_{xy}(\omega)}{P_{yy}(\omega)}$
  • 由于纯净信号无法得到, 在纯净信号和噪声不相关的假设, 可得$P_{xy} = P_{xx}, P_{yy} = P_{xx}+P_{vv}$, 于是$H(\omega)=\frac{P_{xx}}{P_{xx}+P_{vv}}$

Kaiser窗

• 该窗可以同时调节主瓣和旁瓣的宽度, 具有更小的带外能量, 带外衰减大
• 可以低通滤波，同时可以降低高频

BM3D算法

• 2007年TIP的文章Image denoising by sparse 3D transform-domain collaborative ltering

  • 变体：彩色图CBM3D, 时域VBM3D, 沿袭BM3D的基于块处理的思想BM4D和VBM4D只能应用于离线处理

• 算法总体流程:

  • 一、基础估计
    1. 对于每个目标图块，在附近寻找最多MAXN1(超参数)个相似的图块,为了避免噪点的影响，将图块经过2D变换(代码中使用DCT变换)后再用欧式距离衡量相似程度。按距离从小到大排序后取最多前MAXN1个。叠成一个三维数组
       • $d\left(Z_{x_R}, Z_x\right)=\frac{\left\|\Upsilon^{\prime}\left(\mathcal{T}_{2 \mathrm{D}}^{\mathrm{ht}}\left(Z_{x_R}\right)\right)-\Upsilon^{\prime}\left(\mathcal{T}_{2 \mathrm{D}}^{\mathrm{ht}}\left(Z_x\right)\right)\right\|_2^2}{\left(N_1^{\mathrm{ht}}\right)^2}$
    2. 对3D数组的第三维,即图块叠起来后，每个图块同一个位置的像素点构成的数组，进行DCT变换后，采用硬阈值的方式将小于超参数的成分置为0.同时统计非零成分的数量作为后续权重的参考，后续将第三维进行逆变换
       • $\widehat{\mathbf{Y}}_{S_{x_R}^{\mathrm{ht}}}^{\mathrm{ht}}=\mathcal{T}_{3 \mathrm{D}}^{\mathrm{ht}{ }^{-1}}\left(\Upsilon\left(\mathcal{T}_{3 \mathrm{D}}^{\mathrm{ht}}\left(\mathbf{Z}_{S_{x_R}}^{\mathrm{ht}^{\mathrm{t}}}\right)\right)\right)$
       • 为什么这么做？
         • 传统方法由空域得到近似块，然后对近似块的每个像素——对应去平均，作为目标块每个像素的值。但是该策略对部分场景并不合适
           • 某些相似块拥有的噪声更小，相比其它相似块，该块的权重应更大，而不是简单取平均
           • 相似块图像信息冗余，从空域上看，两个有重叠的相似块，简单平均会造成目标块信息重复
           • 采用Collaborative ltering by shrinkage in transform domain能够加强相似块的稀疏性，同时降低相似块的噪声
    3. 将这些图块逆变换后放回原位，利用非零成分数量统计叠加权重，最后将叠放后的图除以每个点的权重就得到基础估计的图像，图像的噪声得到了较大的去除
       • $\hat{y}^{\text {basic }}(x)=\frac{\sum_{x_R \in X} \sum_{x_m \in S_{x_R}^{\mathrm{ht}}} w_{x_R}^{\mathrm{ht}} \widehat{Y}_{x_m}^{\mathrm{ht}, x_R}(x)}{\sum_{x_R \in X} \sum_{x_m \in S_{x_R}^{\mathrm{ht}} .} w_{x_R}^{\mathrm{ht}} \chi_{x_m}(x)}, \forall x \in X$
  • 二、最终估计：基于初步估计,进行改进的分组和协同维纳滤波
    1. 由于基础估计极大地消除了噪点，对于含噪原图的每个目标图块，可以直接用对应基础估计图块的欧氏距离衡量相似程度。按距离从小到大排序后取最多前MAXN1个。将基础估计图块、含噪原图图块分别叠成两个三维数组
       • $S_{x_R}^{\text {wie }}=\left\{x \in X: \frac{\left\|\widehat{Y}_{x_R}^{\text {basic }}-\widehat{Y}_x^{\text {basic }}\right\|_2^2<\tau_{\text {match }}^{\text {wie }}}{\left(N_1^{\text {wie }}\right)^2}\right\}$
    2. 对含基础估计3D数组的第三维，即图块叠起来后，每个图块同一个位置的像素点构成的数组，进行DCT变换，利用如下公式得到系数
       • $\mathbf{W}_{S_{x_R}^{\text {wie }}}=\frac{\left|\mathcal{T}_{3 \mathrm{D}}^{\text {wie }}\left(\widehat{\mathbf{Y}}_{S_{x_R}^{\text {basic }}}^{\text {bus }}\right)\right|^2}{\left|\mathcal{T}_{3 \mathrm{D}}^{\text {wie }}\left(\widehat{\mathbf{Y}}_{S_{x_R}^{\text {wic }}}^{\text {basic }}\right)\right|^2+\sigma^2}$
    3. 将系数与含噪3D图块相乘放回原处，最后做加权平均调整即可得到最终估计图。相对于基础估计图，还原了更多原图的细节

• 算法实现:(见下方代码)

  • 一、基础估计
    • 参数定义: BLK_Size相似窗大小, window_size搜索窗大小, search_step扫描步长, Threshold相似搜索的阈值, max_matched最大堆叠数量, Threshold_hard变换域内三维组系数的阈值
    1. 首先取中心像素的的BLK_Size相似窗大小的邻域图像进行DCT变换, 将变换结果和中心点的位置进行保存
    2. 相似块搜索: 在搜索窗区域进行搜索, 取每个搜索的相似块进行DCT变换, 变换后的结果与原始搜索窗的相似度$distance = \frac{(O-N)^2}{BLK\_Size^2}$, 保存大于阈值$distance > Threshold$相似块的变换结果和坐标, 在排序后取相似度最大的max_matched个结果进行后续操作
    3. 3维滤波: 对于堆叠后的相似块, 其形状定义为[max_matched, BLK_size, BLK_size]. 我们在对每个像素点对应的max_matched矩阵进行DCT变换(即[:, i, j]), 在每次的变换结果中, 将对应像素值小于阈值(Threshold_Hard)的数值置零. 统计这一维的非零元素个数并将置零后的结果iDCT变换后放后对应位置.
    4. 聚合操作: 若非零元素个数小于1, 则将个数置1. 利用非零元素个数和Kaiser窗计算权重, 将上述保存的3维相似块矩阵取每张图像进行iDCT变换([i,:,:]), 随后放回原位置. 每次处理使用累加操作。 在一整张图像的所有像素处理后, 将处理后的结果与权重相除得到基础估计的结果
  • 二、最终估计
    1. 相似块搜索: 与一相同, 但使用一处理后的基础估计结果进行搜索, 最终将基础估计结果的相似块(DCT结果)和含噪原始图像的相似块输出(DCT结果)
    2. 3D维纳协同滤波: 对基础估计图像, 取每个像素点进行DCT变换, 随后计算权重$m_weight = \frac{DCT^T \times DCT}{DCT^T \times DCT + \sigma ^ 2}$, 若该权重不为0, 则反变换的权重为$wiener = \frac{1}{m\_weight^2 \times \sigma^2}$, 随后对含噪原始图像的相似块进行滤波操作, 将其每个像素点进行DCT变换, 利用m_weight进行滤波, 随后经iDCT变换后返回
    3. 聚合操作: 过程类似一, 利用Kaiser窗和反变换的权重对处理后的含噪原始图像的相似块进行iDCT操作后返回, 最终可得到最终结果

• 参数的选择

  • 最大相似块的数量: 第一步硬阈值操作中设置为16, 第二步设置为32
  • 相似块匹配时的阈值: 噪声的$\sigma < 40, \tau_{hard}=\tau_{wien}=2500, else =4000$. 该值对结果具有较小影响, 设置过小导致滤波性能下降, 过大导致假细节
  • 相似块的大小: $k\_{hard}=8, k_{wien}=12$
  • 硬阈值降噪时变换域内三维组系数的阈值$\lambda_{3D}^{hard}$, 该值对结果具有较大的影响, 最优解为2.7

• 加速方法

  • 为了加快计算速度，在寻找相似块的步骤后，得到的块实际上已经进行了2D变换处理，然后再加上一个1D变换(使用1D-Haar离散小波变换),成为了3D变换, 使用2D+1D变换方法替代直接3D变换

• 难点

  • 2D变换与各种超参数并没有一个确定值，需要大量实验积累参数的设置
  • 若实际运用在图像降噪中，原始图像不会大量噪声，因此就不需要BM3D两步去噪
    • 可以将BM3D的两步拆开，采用前步的硬阈值、2D变换寻找相似块、1D变换升至3D域再加权平均，或后步直接使用维纳滤波，或许就已经有很好的效果了

• 参考文章：

  • 图像去噪：BM3D 去噪算法原理解析与代码实现
  • BM3D图像去噪算法原理及代码详解

代码实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343

# -*- coding: utf-8 -*-
import cv2
import numpy


cv2.setUseOptimized(True)

# Parameters initialization
sigma = 25
Threshold_Hard3D = 2.7*sigma           # Threshold for Hard Thresholding
First_Match_threshold = 2500             # 用于计算block之间相似度的阈值
Step1_max_matched_cnt = 16              # 组最大匹配的块数
Step1_Blk_Size = 8                     # block_Size即块的大小，8*8
Step1_Blk_Step = 3                      # Rather than sliding by one pixel to every next reference block, use a step of Nstep pixels in both horizontal and vertical directions.
Step1_Search_Step = 3                   # 块的搜索step
Step1_Search_Window = 39                # Search for candidate matching blocks in a local neighborhood of restricted size NS*NS centered

Second_Match_threshold = 400           # 用于计算block之间相似度的阈值
Step2_max_matched_cnt = 32
Step2_Blk_Size = 8
Step2_Blk_Step = 3
Step2_Search_Step = 3
Step2_Search_Window = 39

Beta_Kaiser = 2.0


def init(img, _blk_size, _Beta_Kaiser):
    """该函数用于初始化，返回用于记录过滤后图像以及权重的数组,还有构造凯撒窗"""
    m_shape = img.shape
    m_img = numpy.matrix(numpy.zeros(m_shape, dtype=float))
    m_wight = numpy.matrix(numpy.zeros(m_shape, dtype=float))
    K = numpy.matrix(numpy.kaiser(_blk_size, _Beta_Kaiser))
    m_Kaiser = numpy.array(K.T * K)            # 构造一个凯撒窗
    return m_img, m_wight, m_Kaiser


def Locate_blk(i, j, blk_step, block_Size, width, height):
    '''该函数用于保证当前的blk不超出图像范围'''
    if i*blk_step+block_Size < width:
        point_x = i*blk_step
    else:
        point_x = width - block_Size

    if j*blk_step+block_Size < height:
        point_y = j*blk_step
    else:
        point_y = height - block_Size

    m_blockPoint = numpy.array((point_x, point_y), dtype=int)  # 当前参考图像的顶点

    return m_blockPoint


def Define_SearchWindow(_noisyImg, _BlockPoint, _WindowSize, Blk_Size):
    """该函数返回一个二元组（x,y）,用以界定_Search_Window顶点坐标"""
    point_x = _BlockPoint[0]  # 当前坐标
    point_y = _BlockPoint[1]  # 当前坐标

    # 获得SearchWindow四个顶点的坐标
    LX = point_x+Blk_Size/2-_WindowSize/2     # 左上x
    LY = point_y+Blk_Size/2-_WindowSize/2     # 左上y
    RX = LX+_WindowSize                       # 右下x
    RY = LY+_WindowSize                       # 右下y

    # 判断一下是否越界
    if LX < 0:   LX = 0
    elif RX > _noisyImg.shape[0]:   LX = _noisyImg.shape[0]-_WindowSize
    if LY < 0:   LY = 0
    elif RY > _noisyImg.shape[0]:   LY = _noisyImg.shape[0]-_WindowSize

    return numpy.array((LX, LY), dtype=int)


def Step1_fast_match(_noisyImg, _BlockPoint):
    """快速匹配"""
    '''
    *返回邻域内寻找和当前_block相似度最高的几个block,返回的数组中包含本身
    *_noisyImg:噪声图像
    *_BlockPoint:当前block的坐标及大小
    '''
    (present_x, present_y) = _BlockPoint  # 当前坐标
    Blk_Size = Step1_Blk_Size
    Search_Step = Step1_Search_Step
    Threshold = First_Match_threshold
    max_matched = Step1_max_matched_cnt
    Window_size = Step1_Search_Window

    blk_positions = numpy.zeros((max_matched, 2), dtype=int)  # 用于记录相似blk的位置
    Final_similar_blocks = numpy.zeros((max_matched, Blk_Size, Blk_Size), dtype=float)

    img = _noisyImg[present_x: present_x+Blk_Size, present_y: present_y+Blk_Size]
    dct_img = cv2.dct(img.astype(numpy.float64))  # 对目标作block作二维余弦变换

    Final_similar_blocks[0, :, :] = dct_img
    blk_positions[0, :] = _BlockPoint

    Window_location = Define_SearchWindow(_noisyImg, _BlockPoint, Window_size, Blk_Size)
    blk_num = (Window_size-Blk_Size)/Search_Step  # 确定最多可以找到多少相似blk
    blk_num = int(blk_num)
    (present_x, present_y) = Window_location

    similar_blocks = numpy.zeros((blk_num**2, Blk_Size, Blk_Size), dtype=float)
    m_Blkpositions = numpy.zeros((blk_num**2, 2), dtype=int)
    Distances = numpy.zeros(blk_num**2, dtype=float)  # 记录各个blk与它的相似度

    # 开始在_Search_Window中搜索,初始版本先采用遍历搜索策略,这里返回最相似的几块
    matched_cnt = 0
    for i in range(blk_num):
        for j in range(blk_num):
            tem_img = _noisyImg[present_x: present_x+Blk_Size, present_y: present_y+Blk_Size]
            dct_Tem_img = cv2.dct(tem_img.astype(numpy.float64))
            m_Distance = numpy.linalg.norm((dct_img-dct_Tem_img))**2 / (Blk_Size**2)

            # 下面记录数据自动不考虑自身(因为已经记录)
            if m_Distance < Threshold and m_Distance > 0:  # 说明找到了一块符合要求的
                similar_blocks[matched_cnt, :, :] = dct_Tem_img
                m_Blkpositions[matched_cnt, :] = (present_x, present_y)
                Distances[matched_cnt] = m_Distance
                matched_cnt += 1
            present_y += Search_Step
        present_x += Search_Step
        present_y = Window_location[1]
    Distances = Distances[:matched_cnt]
    Sort = Distances.argsort()

    # 统计一下找到了多少相似的blk
    if matched_cnt < max_matched:
        Count = matched_cnt + 1
    else:
        Count = max_matched

    if Count > 0:
        for i in range(1, Count):
            Final_similar_blocks[i, :, :] = similar_blocks[Sort[i-1], :, :]
            blk_positions[i, :] = m_Blkpositions[Sort[i-1], :]
    return Final_similar_blocks, blk_positions, Count


def Step1_3DFiltering(_similar_blocks):
    '''
    *3D变换及滤波处理
    *_similar_blocks:相似的一组block,这里已经是频域的表示
    *要将_similar_blocks第三维依次取出,然在频域用阈值滤波之后,再作反变换
    '''
    statis_nonzero = 0  # 非零元素个数
    m_Shape = _similar_blocks.shape

    # 下面这一段代码很耗时
    for i in range(m_Shape[1]):
        for j in range(m_Shape[2]):
            tem_Vct_Trans = cv2.dct(_similar_blocks[:, i, j])
            tem_Vct_Trans[numpy.abs(tem_Vct_Trans[:]) < Threshold_Hard3D] = 0.
            statis_nonzero += tem_Vct_Trans.nonzero()[0].size
            _similar_blocks[:, i, j] = cv2.idct(tem_Vct_Trans)[0]
    return _similar_blocks, statis_nonzero


def Aggregation_hardthreshold(_similar_blocks, blk_positions, m_basic_img, m_wight_img, _nonzero_num, Count, Kaiser):
    '''
    *对3D变换及滤波后输出的stack进行加权累加,得到初步滤波的图片
    *_similar_blocks:相似的一组block,这里是频域的表示
    *对于最后的数组，乘以凯撒窗之后再输出
    '''
    _shape = _similar_blocks.shape
    if _nonzero_num < 1:
        _nonzero_num = 1
    block_wight = (1./_nonzero_num) * Kaiser
    for i in range(Count):
        point = blk_positions[i, :]
        tem_img = (1./_nonzero_num)*cv2.idct(_similar_blocks[i, :, :]) * Kaiser
        m_basic_img[point[0]:point[0]+_shape[1], point[1]:point[1]+_shape[2]] += tem_img
        m_wight_img[point[0]:point[0]+_shape[1], point[1]:point[1]+_shape[2]] += block_wight


def BM3D_1st_step(_noisyImg):
    """第一步,基本去噪"""
    # 初始化一些参数：
    (width, height) = _noisyImg.shape   # 得到图像的长宽
    block_Size = Step1_Blk_Size         # 块大小
    blk_step = Step1_Blk_Step           # N块步长滑动
    Width_num = (width - block_Size)/blk_step
    Height_num = (height - block_Size)/blk_step

    # 初始化几个数组
    Basic_img, m_Wight, m_Kaiser = init(_noisyImg, Step1_Blk_Size, Beta_Kaiser)

    # 开始逐block的处理,+2是为了避免边缘上不够
    for i in range(int(Width_num+2)):
        for j in range(int(Height_num+2)):
            # m_blockPoint当前参考图像的顶点
            m_blockPoint = Locate_blk(i, j, blk_step, block_Size, width, height)       # 该函数用于保证当前的blk不超出图像范围
            Similar_Blks, Positions, Count = Step1_fast_match(_noisyImg, m_blockPoint)
            Similar_Blks, statis_nonzero = Step1_3DFiltering(Similar_Blks)
            Aggregation_hardthreshold(Similar_Blks, Positions, Basic_img, m_Wight, statis_nonzero, Count, m_Kaiser)
    Basic_img[:, :] /= m_Wight[:, :]
    basic = numpy.matrix(Basic_img, dtype=int)
    basic.astype(numpy.uint8)

    return basic


def Step2_fast_match(_Basic_img, _noisyImg, _BlockPoint):
    '''
    *快速匹配算法,返回邻域内寻找和当前_block相似度最高的几个block,要同时返回basicImg和IMG
    *_Basic_img: 基础去噪之后的图像
    *_noisyImg:噪声图像
    *_BlockPoint:当前block的坐标及大小
    '''
    (present_x, present_y) = _BlockPoint  # 当前坐标
    Blk_Size = Step2_Blk_Size
    Threshold = Second_Match_threshold
    Search_Step = Step2_Search_Step
    max_matched = Step2_max_matched_cnt
    Window_size = Step2_Search_Window

    blk_positions = numpy.zeros((max_matched, 2), dtype=int)  # 用于记录相似blk的位置
    Final_similar_blocks = numpy.zeros((max_matched, Blk_Size, Blk_Size), dtype=float)
    Final_noisy_blocks = numpy.zeros((max_matched, Blk_Size, Blk_Size), dtype=float)

    img = _Basic_img[present_x: present_x+Blk_Size, present_y: present_y+Blk_Size]
    dct_img = cv2.dct(img.astype(numpy.float32))  # 对目标作block作二维余弦变换
    Final_similar_blocks[0, :, :] = dct_img

    n_img = _noisyImg[present_x: present_x+Blk_Size, present_y: present_y+Blk_Size]
    dct_n_img = cv2.dct(n_img.astype(numpy.float32))  # 对目标作block作二维余弦变换
    Final_noisy_blocks[0, :, :] = dct_n_img

    blk_positions[0, :] = _BlockPoint

    Window_location = Define_SearchWindow(_noisyImg, _BlockPoint, Window_size, Blk_Size)
    blk_num = (Window_size-Blk_Size)/Search_Step  # 确定最多可以找到多少相似blk
    blk_num = int(blk_num)
    (present_x, present_y) = Window_location

    similar_blocks = numpy.zeros((blk_num**2, Blk_Size, Blk_Size), dtype=float)
    m_Blkpositions = numpy.zeros((blk_num**2, 2), dtype=int)
    Distances = numpy.zeros(blk_num**2, dtype=float)  # 记录各个blk与它的相似度

    # 开始在_Search_Window中搜索,初始版本先采用遍历搜索策略,这里返回最相似的几块
    matched_cnt = 0
    for i in range(blk_num):
        for j in range(blk_num):
            tem_img = _Basic_img[present_x: present_x+Blk_Size, present_y: present_y+Blk_Size]
            dct_Tem_img = cv2.dct(tem_img.astype(numpy.float32))
            m_Distance = numpy.linalg.norm((dct_img-dct_Tem_img))**2 / (Blk_Size**2)

            # 下面记录数据自动不考虑自身(因为已经记录)
            if m_Distance < Threshold and m_Distance > 0:
                similar_blocks[matched_cnt, :, :] = dct_Tem_img
                m_Blkpositions[matched_cnt, :] = (present_x, present_y)
                Distances[matched_cnt] = m_Distance
                matched_cnt += 1
            present_y += Search_Step
        present_x += Search_Step
        present_y = Window_location[1]
    Distances = Distances[:matched_cnt]
    Sort = Distances.argsort()

    # 统计一下找到了多少相似的blk
    if matched_cnt < max_matched:
        Count = matched_cnt + 1
    else:
        Count = max_matched

    if Count > 0:
        for i in range(1, Count):
            Final_similar_blocks[i, :, :] = similar_blocks[Sort[i-1], :, :]
            blk_positions[i, :] = m_Blkpositions[Sort[i-1], :]

            (present_x, present_y) = m_Blkpositions[Sort[i-1], :]
            n_img = _noisyImg[present_x: present_x+Blk_Size, present_y: present_y+Blk_Size]
            Final_noisy_blocks[i, :, :] = cv2.dct(n_img.astype(numpy.float64))

    return Final_similar_blocks, Final_noisy_blocks, blk_positions, Count


def Step2_3DFiltering(_Similar_Bscs, _Similar_Imgs):
    '''
    *3D维纳变换的协同滤波
    *_similar_blocks:相似的一组block,这里是频域的表示
    *要将_similar_blocks第三维依次取出,然后作dct,在频域进行维纳滤波之后,再作反变换
    *返回的Wiener_wight用于后面Aggregation
    '''
    m_Shape = _Similar_Bscs.shape
    Wiener_wight = numpy.zeros((m_Shape[1], m_Shape[2]), dtype=float)

    for i in range(m_Shape[1]):
        for j in range(m_Shape[2]):
            tem_vector = _Similar_Bscs[:, i, j]
            tem_Vct_Trans = numpy.matrix(cv2.dct(tem_vector))
            Norm_2 = numpy.float64(tem_Vct_Trans.T * tem_Vct_Trans)
            m_weight = Norm_2/(Norm_2 + sigma**2)
            if m_weight != 0:
                Wiener_wight[i, j] = 1./(m_weight**2 * sigma**2)
            # else:
            #     Wiener_wight[i, j] = 10000
            tem_vector = _Similar_Imgs[:, i, j]
            tem_Vct_Trans = m_weight * cv2.dct(tem_vector)
            _Similar_Bscs[:, i, j] = cv2.idct(tem_Vct_Trans)[0]

    return _Similar_Bscs, Wiener_wight


def Aggregation_Wiener(_Similar_Blks, _Wiener_wight, blk_positions, m_basic_img, m_wight_img, Count, Kaiser):
    '''
    *对3D变换及滤波后输出的stack进行加权累加,得到初步滤波的图片
    *_similar_blocks:相似的一组block,这里是频域的表示
    *对于最后的数组，乘以凯撒窗之后再输出
    '''
    _shape = _Similar_Blks.shape
    block_wight = _Wiener_wight # * Kaiser

    for i in range(Count):
        point = blk_positions[i, :]
        tem_img = _Wiener_wight * cv2.idct(_Similar_Blks[i, :, :]) # * Kaiser
        m_basic_img[point[0]:point[0]+_shape[1], point[1]:point[1]+_shape[2]] += tem_img
        m_wight_img[point[0]:point[0]+_shape[1], point[1]:point[1]+_shape[2]] += block_wight


def BM3D_2nd_step(_basicImg, _noisyImg):
    '''Step 2. 最终的估计: 利用基本的估计，进行改进了的分组以及协同维纳滤波'''
    # 初始化一些参数：
    (width, height) = _noisyImg.shape
    block_Size = Step2_Blk_Size
    blk_step = Step2_Blk_Step
    Width_num = (width - block_Size)/blk_step
    Height_num = (height - block_Size)/blk_step

    # 初始化几个数组
    m_img, m_Wight, m_Kaiser = init(_noisyImg, block_Size, Beta_Kaiser)

    for i in range(int(Width_num+2)):
        for j in range(int(Height_num+2)):
            m_blockPoint = Locate_blk(i, j, blk_step, block_Size, width, height)
            Similar_Blks, Similar_Imgs, Positions, Count = Step2_fast_match(_basicImg, _noisyImg, m_blockPoint)
            Similar_Blks, Wiener_wight = Step2_3DFiltering(Similar_Blks, Similar_Imgs)
            Aggregation_Wiener(Similar_Blks, Wiener_wight, Positions, m_img, m_Wight, Count, m_Kaiser)
    m_img[:, :] /= m_Wight[:, :]
    Final = numpy.matrix(m_img, dtype=int)
    Final.astype(numpy.uint8)

    return Final

变换域去噪

• 基于频域的处理方法主要是用滤波器，把有用的信号和干扰信号分开，它在有用信号和干扰信号的频谱没有重叠的前提下，才能把有用信号和干扰信号完全分开，但实际上很难分开
  • 主要方法：傅里叶变换，小波变换等