图像处理 深度学习 实战 图像处理算法研究

一些基本数字图像处理算法

版权声明：本文为原创文章，未经博主允许不得用于商业用途。

所有的图像算法都在DIPAlgorithm类中，并且所有算法都为抽象成员函数。我已经按照java注释规范为所有方法添加使用说明注释，具体实现可见于DIPAlgorithm.java，这里只做算法说明。

1 图像扭曲

模仿PS的扭曲功能，通过建立一个三角形映射网格实现对图像的扭曲。

如上图，一共设置了45个控制点围成74个三角形网格

扭曲即形变处理其实是寻找一个函数，以所有网格顶点原始坐标为输入，扭曲后所有网格顶点坐标为输出。为了简化计算任务，采用控制栅格插值方法，对每个三角网格独立计算映射关系，如下图：

即求解矩阵$M$满足$MA = B$，其中$A$为原顶点的齐次矩阵：

$A = \begin{bmatrix} x_{1} & y_{1} & 1 \\ x_{2} & y_{2} & 1 \\ x_{3} & y_{3} & 1 \\ \end{bmatrix}$

B为形变后顶点的其次矩阵：

$B = \lbrack\begin{matrix} x_{1}^{$

M即为$2 \times 3$的映射矩阵，且由于三角形三点不共线，因此A为可逆阵，

$M = BA^{- 1}$

对于三角形中的点$p\left( x,\ y \right)$，其映射后坐标$p^{$

2 直方图计算

直方图计算实际上即求图像的概率密度函数PDF，只需遍历一次所有像素点即可获得。

3 直方图均衡化算法

对于连续图像直方图均衡化其实是种点运算f,
对不同灰度值做映射，使得所有像素频率相等。

对于点运算f，有如下性质：

$D_{B} = f\left( D_{A} \right),\ H_{B}\left( D_{B} \right)\Delta D_{B} = H_{A}\left( D_{A} \right)\Delta D_{A}$

其中D为灰度值，H即为灰度值在图像中的频数，整理可得

$H_{B}\left( D_{B} \right) = \frac{H_{A}\left( D_{A} \right)\Delta D_{A}}{\Delta D_{B}} = \frac{H_{A}\left( D_{A} \right)}{\frac{\Delta D_{B}}{\Delta D_{A}}} = \frac{H_{A}\left( D_{A} \right)}{\frac{dD_{B}}{dD_{A}}}$

$= \frac{H_{A}\left( D_{A} \right)}{f$

即：

寻找函数f使得$H_{B}(D)$为常数$\frac{A_{0}}{D_{m}},A_{0},D_{m}$。

由(1)可知，KaTeX parse error: Undefined control sequence: \ at position 58: …\right)}{f'(D)}\̲ ̲\Rightarrow f^{…

即$f\left( D \right) = D_{m}CDF(D)$，CDF即累积分布函数

因此只需求得直方图的前序和即可获得映射关系。

4 图像灰度化

目前比较符合人眼的灰度化权重为0.299、0.578和0.114，为了加速计算使用近似公式$D = (3r + g + 6b)/10$

5 图像二值化

我使用的二值化算法为OSTU大律二值化算法。二值化操作即利用分割阈值u，将图片分为前景后景两部分。OSTU大律法认为使得前景像素和背景像素灰度方差g最大的阈值即为最佳分割阈值。

$g = w_{0}w_{1}\left( u_{0} - u_{1} \right)^{2}$

其中$w_{0},\ w_{1}$为前景、后景在图像中的比例，KaTeX parse error: Undefined control sequence: \ at position 7: u_{0},\̲ ̲u_{1}为前景、后景的平均灰度。

在实现时，只需遍历所有灰度，利用CDF求出每种灰度的方差，取最大者作为阈值即可。

6 前景分离

目前主流的前景分离为深度学习算法。这里只使用了最基本的阈值分离法，分别为RGB三个通道设置不同阈值，将小于阈值的像素作为背景，大于阈值的作为前景。

7 滤波

我使用的滤波方法是高斯滤波和中值滤波，高斯滤波即使用二维高斯函数作为滤波函数，中值滤波即使用邻域的中位数作为滤波函数。

高斯滤波器为线性滤波器，可以有效消除高斯噪声。由于高斯函数离中值越近权重越大，因此相对于均值滤波器更加柔和，对边缘的保留效果更好。这里我使用的是如下矩阵做卷积：

$\begin{bmatrix} 1 & 2 & 3 & 2 & 1 \\ 2 & 4 & 6 & 4 & 2 \\ 3 & 6 & 7 & 6 & 3 \\ 2 & 4 & 6 & 4 & 2 \\ 1 & 2 & 3 & 2 & 1 \\ \end{bmatrix}$

中值滤波器为非线性滤波器，可以有效的去除椒盐噪声和斑点噪声并且不会使图像变模糊。

8 形态学扩张和腐蚀

形态学腐蚀可记为$\text{AΘB}$，其中A为输入图像，B为结构单元。对于二值图像，当且仅当当前像素点满足腐结构单元时才会被保留。对于灰度图像，则可类比为最小值，即

$f\Theta b\left( x,y \right) = min\{ f\left( x - x^{$

形态学扩张可看作腐蚀的逆操作,记作$A\bigoplus B$，对于二值图像，将每个有效像素点的邻域结构单元置1，对于灰度图像则取最大值，即

$f\bigoplus b\left( x,y \right) = max\{ f\left( x - x^{$

本程序将结构单元b统一设定为5*5矩形。

通过扩张和腐蚀的结合可实现结构开运算（$AoB = \left( \text{AΘB} \right)\bigoplus B$）和结构闭运算（$AoB = \left( A\bigoplus B \right)\text{ΘB}$）对图像进行粗化、细化、滤波等处理

9 傅里叶变换和滤波

变换公式

傅里叶变换可以将信号从时域转换到频域，因此可以看出许多时域中不明显的特征。二维傅里叶变换（CFT）公式如下：

$F\left( u,v \right) = \iint_{}^{}{f\left( x,y \right)e^{- 2\pi\overrightarrow{j}(ux + vy)}}\text{dxdy}$

其中${\overrightarrow{j}}^{2} = - 1,f,F$，同样二维傅里叶逆变换公式如下:

$f\left( x,y \right) = \iint_{}^{}{F\left( u,v \right)e^{2\pi\overrightarrow{j}(ux + vy)}}\text{dudv}$

对于离散函数，可以定义离散二维傅里叶变换（DFT）和逆变换：

$G\left( m,n \right) = \frac{1}{\sqrt{\text{MN}}}\sum_{\begin{matrix} 0 \leq \ i\ \leq \ M - 1 \\ 0 < k < N - 1\ \\ \end{matrix}}^{}{g\left( i,k \right)e^{- 2\pi\overrightarrow{j}(\frac{\text{im}}{M} + \frac{\text{jn}}{N})}}$

$g\left( i,k \right) = \frac{1}{\sqrt{\text{MN}}}\sum_{\begin{matrix} 0 \leq \ m\ \leq \ M - 1 \\ 0 < n < N - 1\ \\ \end{matrix}}^{}{g\left( m,n \right)e^{2\pi\overrightarrow{j}(\frac{\text{im}}{M} + \frac{\text{jn}}{N})}}$

DFT可以理解为对连续二维信号进行了频率为M,
N的采样，之后通过计算其和频域空间M*N个基向量的相关性（在该方向投影）将时域信号映射到频域。iDFT可以理解为通过M*N个基向量合成原始时域信号。

矩阵表示

傅里叶变换实际上是一种线性变换，因此在实际计算中常常将$g$扩充为$N*N$方阵，此时DFT可以通过矩阵表示：$G = \mathcal{W}^{- 1}g\mathcal{W},\mathcal{W}_{\text{ik}} = \frac{1}{N}e^{2\pi\overrightarrow{j}\frac{\text{ik}}{N}}$。

易知$\mathcal{W}_{\text{ik}} = \mathcal{W}_{\text{ki}}$，且为正交矩阵，因此$\mathcal{W}$为酉矩阵，即$\mathcal{W}^{- 1} = \left( \mathcal{W}^{*} \right)^{T} = \mathcal{W}^{*}$，$G = \mathcal{W}^{*}g\mathcal{W}$，其中$F^{*}F$。

由于傅里叶变换为酉变换，即$\mathcal{W}^{t} = \mathcal{W}^{- 1}$

图像的傅里叶变换

对于二维图片可以看作二维矩阵，因此可以进行DFT。二维图片经过DFT后获得的复矩阵的模矩阵可以表示每个频率信号的强度（也可看作先做自相关后再进行傅里叶变换），经过适当处理即可转化为灰度能量谱图片。

线性噪声在频域中通常为点或线，因此可以通过傅里叶变换后进行滤波再通过逆变换复原图片。

算法实现

在实际实现时，根据欧拉公式，$e^{- \overrightarrow{j}t} = cost - \overrightarrow{j}\text{sint},\ e^{\overrightarrow{j}t} = cost + \overrightarrow{j}\text{sint}$，因此傅里叶变换的核矩阵可以表示为$\mathcal{W}_{\text{ik}} = \frac{\cos\left( 2\pi ik \right) - \overrightarrow{j}\sin\left( 2\pi ik \right)}{N}$，为方便运算将$\mathcal{W}$分解为虚部系数$\mathcal{W}_{\text{lm}}$和实部系数$\mathcal{W}_{\text{re}}$，其中则$\mathcal{W} = \mathcal{W}_{\text{re}} + \overrightarrow{j}\mathcal{W}_{\text{lm}}$。变换结果同样分解为$G = G_{\text{re}} + \overrightarrow{j}G_{\text{lm}}$，则DFT可以表示为:

$G = \mathcal{W}^{*}g\mathcal{W =}\left( \mathcal{W}_{\text{re}} - \overrightarrow{j}\mathcal{W}_{\text{lm}} \right)g\left( \mathcal{W}_{\text{re}} + \overrightarrow{j}\mathcal{W}_{\text{lm}} \right) = \mathcal{W}_{\text{re}}g\mathcal{W}_{\text{re}} + \mathcal{W}_{\text{lm}}g\mathcal{W}_{\text{lm}} - \overrightarrow{j}\left( \mathcal{W}_{\text{lm}}g\mathcal{W}_{\text{re}} + \mathcal{W}_{\text{re}}g\mathcal{W}_{\text{lm}} \right)$

$\left\{ \begin{matrix} G_{\text{re}} = \mathcal{W}_{\text{re}}g\mathcal{W}_{\text{re}} + \mathcal{W}_{\text{lm}}g\mathcal{W}_{\text{lm}} \\ G_{\text{lm}} = - \mathcal{W}_{\text{lm}}g\mathcal{W}_{\text{re}} - \mathcal{W}_{\text{re}}g\mathcal{W}_{\text{lm}} \\ \end{matrix} \right.\$

同理，iDFT可以表示为：

$g = \left( \mathcal{W}_{\text{re}} + \overrightarrow{j}\mathcal{W}_{\text{lm}} \right)(G_{\text{re}} + {\overrightarrow{j}G}_{\text{lm}})\left( \mathcal{W}_{\text{re}} - \overrightarrow{j}\mathcal{W}_{\text{lm}} \right)$

其中，为了将能量谱转化为可见的灰度图，为能量谱取对数值进行归一化。且由于在频域中两个维度频率都为0时（即$\mathcal{W}_{00}$处）为图像能量的总和，因此通过$log(e + 1)*\frac{256}{\log\left( \mathcal{W}_{00} + 1 \right)}$可以做进一步归一化。

算法代码可见github