[Python图像处理] 基于曲率驱动算法的图像去噪

发布时间：2023年12月29日

基于曲率驱动算法的图像去噪

前言

在本节中，我们将学习如何使用曲率驱动的滤波器和各向异性扩散滤波器来从带有噪声的图像中删除噪声。

图像梯度

图像I的梯度κ可以定义如下：

$\begin{aligned} k(I)&=\nabla \cdot (\frac {\nabla I} {|\nabla I|}) \\ &=\frac \partial {\partial x}(\frac {I_x} {\sqrt {I_x^2+I_y^2}})+\frac \partial {\partial y}(\frac {I_y} {\sqrt {I_x^2+I_y^2}}) \\ &=\frac {{\sqrt {I_x^2+I_y^2}} \cdot I_{xx}-I_x\cdot \frac {2\cdot I_x \cdot I_{xx}+2\cdot I_y \cdot I_{xy}} {2\cdot \sqrt {I_x^2+I_y^2}}} {I_x^2+I_y^2} + \frac {{\sqrt {I_x^2+I_y^2}} \cdot I_{yy}-I_y\cdot \frac {2\cdot I_x \cdot I_{yx}+2\cdot I_y \cdot I_{yy}} {2\cdot \sqrt {I_x^2+I_y^2}}} {I_x^2+I_y^2} \\ &=\frac {{I_x^2\cdot I_{yy}+I_y^2\cdot I_{xx} -2\cdot I_y} \cdot I_{xy}} {(I_x^2+I_y^2)^{\frac 3 2}} \end{aligned}$

其中 $I_x=\frac {\partial I}{\partial x}$ ， $I_y=\frac {\partial I}{\partial y}$ ， $\nabla=\left( \begin{matrix} \frac \partial {\partial x} \\ \\ \frac \partial {\partial y} \end{matrix} \right)$ 。

SimpleITK 库的 CurvatureFlowImageFilter 和 MinMaxCurvatureFlowImageFilter 实现了曲率驱动的图像去噪算法。灰度输入图像中的等亮度轮廓被视为一个级别集，然后，使用基于曲率的速度函数来演化级别集。对于速度函数，其中 $κ$ 是点 $x$ 处等亮度轮廓的平均曲率：

$F_{minmax}=\left\{ \begin{array}{rcl} max(k,0) & & {Avgstencil(x)≤T_{threshold}}\\ min(k,0) & & {otherwise} \end{array} \right.$

曲率驱动的图像滤波器的优势在于，仅在一个区域内进行平滑，从而保留了清晰的边界。但是，应该注意的是，随着每个轮廓收缩至零并消失，该滤波器的连续应用将导致最终所有信息被删除。MinmaxCurvatureFlowMageFilter 使用最小/最大曲率流对图像进行降噪。

曲率驱动的滤波器

在本节中，我们将使用灰度输入图像，向其添加随机噪声，并使用曲率驱动的滤波器有效地从图像中删除噪声，并尽可能地保留边缘。

(1) 导入所需库，然后读取输入图像：

import SimpleITK as sitk
import matplotlib.pylab as plt

img = sitk.ReadImage('1.png', sitk.sitkFloat64)

(2) 使用 NormalizeImageFilter() 函数对输入图像执行归一化，将输入图像的平均值置为零，方差置为 1。使用 ShotNoiseImageFilter() 函数向图像添加散弹噪声 (shot noise)。

 normfilter = sitk.NormalizeImageFilter()
caster = sitk.CastImageFilter()
caster.SetOutputPixelType(sitk.sitkFloat64)
        
tkfilter = sitk.ShotNoiseImageFilter()
tkfilter.SetScale(0.2)
img_noisy = tkfilter.Execute (img)
img_noisy = sitk.RescaleIntensity(img_noisy)

散弹噪声遵循泊松分布，噪声取决于图像中的像素强度。从带有噪声的图像中删除噪声，过滤器 CurvatureFlowImageFilter() 具有两个参数，即要执行的更新迭代次数以及每次更新之间的时间步长。
时间步应该足够小，以确保数值稳定性。当时间步长满足 Courant-Friedrichs levy 条件时，稳定性可以得到保证。广义上讲，该条件确保每个轮廓在每个时间步长中移动的网格位置不超过一个。时间步长通常由用户指定，必须根据不同应用手动调整。

(3) 滤波器利用多线程有限差分求解器层次结构，使用 CurvatureFlowFunction 对象和零通量 Neumann 边界条件计算更新：

tkfilter = sitk.CurvatureFlowImageFilter() 
tkfilter.SetNumberOfIterations(50)
tkfilter.SetTimeStep(0.1)
img_res_TK = tkfilter.Execute(img_noisy)

(4) 接下来，使用 MinMaxCurvatureFlowImageFilter() 函数利用最小/最大曲率流对图像进行去噪。
在最小/最大曲率流中，根据要移除的噪波的比例启用或禁用移动 movement；切换 switching 取决于每个点周围半径为 R 的区域的平均图像值；R (模板半径)的选择控制要去除的噪声的比例；阈值 threshold 是在与局部邻域极值处的点x处的梯度垂直的方向上获得的平均强度：

tkfilter = sitk.MinMaxCurvatureFlowImageFilter() 
tkfilter.SetNumberOfIterations(50)
tkfilter.SetTimeStep(0.1)
tkfilter.SetStencilRadius(4)
img_res_TK1 = tkfilter.Execute(img_noisy)![请添加图片描述](https://img-blog.csdnimg.cn/b5078592e19743a6ada1125b8e03c25f.jpeg)

img_res_TK1 = sitk.RescaleIntensity(img_res_TK1)

各向异性扩散

各向异性扩散方法用于减少图像中的噪声(或不需要的细节)，同时保留特定的图像特征。在许多应用程序中，均关注明暗过渡部分(边缘)。标准各向同性扩散方法会移动和模糊明暗边界，而各向异性扩散方法可以保留图像边缘。
这些方法可以视为计算图像多尺度描述的工具。将图像 $U (x)$ 嵌入导出图像的较高维函数，该高维函数表示热扩散方程的解：

$\frac {dU(x)} {dt}=\nabla\cdot C\nabla U(x)$

在常数 $C$ 和初始条件下， $U (x, 0)$ 表示原始图像。将 $C$ 由常数扩展到函数的情况，可以得到以下方程：

$\frac {dU(x)} {dt}=C(x)\Delta U(x)+\nabla C(x)\nabla U(x)$

参数 C (conductance) 的选择会改变各向异性扩散的强度。通常，选择 C 作为图像特征的某些函数，以选择性地保留或删除这些特征。例如，边缘保留倾向于在更平滑的区域上进行，其中 C 根据梯度大小进行反向缩放：

$C(x)=e^{-(\frac {||\nabla U(x)||}{k})^2}$

在 ITK 中作为子类实现了以上方程，可以通过使用迭代有限正差分技术求解方程。

各向异性扩散滤波器

接下来，使用上述两个各向异性扩散滤波器实现图像去噪。

(1) 首先使用 CurvatureAnisotropicDiffusionImageFilter()，该滤波器使用修正曲率扩散方程对标量图像执行各向异性扩散，此滤波器的默认时间步长设置为最大理论稳定值，其中 N 是图像的维度：

tkfilter = sitk.CurvatureAnisotropicDiffusionImageFilter()
tkfilter.SetNumberOfIterations(100);
tkfilter.SetTimeStep(0.05);
tkfilter.SetConductanceParameter(3);
img_res_TK2 = tkfilter.Execute(img_noisy)
#img_res_TK1 = sitk.RescaleIntensity(img_res_TK1)

(2) 接下来，使用 GradientAnisotropicDiffusionMimageFilter() 函数从输入图像中移除噪声，获得输出图像。该滤波器使用基于梯度大小的方程在标量图像上执行各向异性扩散，在各向异性扩散滤波器层次结构中，时间步长将由用户明确设置，此处所指的时间步长完全对应于有限差分更新方程中的 $Δ T$ 。求解这类偏微分方程的适当时间步长取决于图像的维数和方程的阶数。该滤波器会自动尝试将其时间步长限制为稳定值，如果时间步长设置得太高，则会生成运行时警告。
参数 C 控制项 C 的灵敏度，值越低，扩散方程对图像特征的保留越强，较高的值将使过滤器更容易扩散图像特征，参数 C 通常设置在范围 0.5 至 2.0 之间：

tkfilter = sitk.GradientAnisotropicDiffusionImageFilter()
tkfilter.SetNumberOfIterations(100);
tkfilter.SetTimeStep(0.05);
tkfilter.SetConductanceParameter(3);
img_res_TK3 = tkfilter.Execute(img_noisy)

(3) 最后，绘制原始/噪声输入和使用不同算法获得的去噪输出图像：

plt.figure(figsize=(16,20))
plt.gray()
plt.subplots_adjust(0,0,1,1,0.01,0.05)
plt.subplot(321), plt.imshow(sitk.GetArrayFromImage(img)), plt.axis('off'), plt.title('Original', size=10)
plt.subplot(322), plt.imshow(sitk.GetArrayFromImage(img_noisy)), plt.axis('off'), plt.title('Noisy (with added Shot Noise)', size=10)
plt.subplot(323), plt.imshow(sitk.GetArrayFromImage(img_res_TK)), plt.axis('off'), plt.title('Denoised (with CurvatureFlowImageFilter)', size=10)
plt.subplot(324), plt.imshow(sitk.GetArrayFromImage(img_res_TK1)), plt.axis('off'), plt.title('Denoised (with MinMaxCurvatureFlowImageFilter)', size=10)
plt.subplot(325), plt.imshow(sitk.GetArrayFromImage(img_res_TK2)), plt.axis('off'), plt.title('Denoised (with CurvatureAnisotropicDiffusionImageFilter)', size=10)
plt.subplot(326), plt.imshow(sitk.GetArrayFromImage(img_res_TK3)), plt.axis('off'), plt.title('Denoised (with GradientAnisotropicDiffusionImageFilter)', size=10)
plt.show()

图像去噪