大家好，又见面了，我是你们的朋友全栈君。

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

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

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

1 图像扭曲

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

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

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

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

$$A=\left[\begin{array}{ccc}{x}_1& y_1& 1\\ x_2& y_2& 1\\ x_3& y_3& 1\end{array}\right]$$

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

$$B=[\begin{array}{ccc}{x}_1^{{}^{\prime }}& x_2^{{}^{\prime }}& x_3^{{}^{\prime }}\\ y_1^{{}^{\prime }}& y_2^{{}^{\prime }}& y_3^{{}^{\prime }}\end{array}]$$

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

$$M=B{A}^{-1}$$

对于三角形中的点$$p\left(x,y\right)$$，其映射后坐标$$p^{{}^{\prime }}=M\left[\begin{array}{c}x\\ y\\ 1\end{array}\right]$$

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{d{D}_B}{d{D}_A}}$$

$$=\frac{H_A\left(D_A\right)}{f^{\prime }(D_A)}=\frac{H_A\left(f^{-1}\left(D_B\right)\right)}{f^{\prime }(f^{-1}(D_B))}$$

即：

寻找函数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 滤波

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

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

$$\left[\begin{array}{ccccc}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{array}\right]$$

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

8 形态学扩张和腐蚀

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

f Θ b ( x , y ) = m i n { f ( x − x ′ ,   y − y ′ ) − b ( x ′ , y ′ ) ∣ ( x ′ , y ′ ∈ D b ) } f\Theta b\left( x,y \right) = min\{ f\left( x – x^{‘},\ y – y^{‘} \right) – b(x^{‘},y’)|(x^{‘},y^{‘} \in D_{b})\} fΘb(x,y)=min{
f(x−x′, y−y′)−b(x′,y′)∣(x′,y′∈Db​)}

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

f ⨁ b ( x , y ) = m a x { f ( x − x ′ ,   y − y ′ ) − b ( x ′ , y ′ ) ∣ ( x ′ , y ′ ∈ D b ) } f\bigoplus b\left( x,y \right) = max\{ f\left( x – x^{‘},\ y – y^{‘} \right) – b(x^{‘},y’)|(x^{‘},y^{‘} \in D_{b})\} f⨁b(x,y)=max{
f(x−x′, y−y′)−b(x′,y′)∣(x′,y′∈Db​)}

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

通过扩张和腐蚀的结合可实现结构开运算（$$AoB=\left(\text{AΘB}\right)⨁B$$）和结构闭运算（$$AoB=\left(A⨁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{array}{c}0\le i\le M-1\\ 0<k<N-1\end{array}}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{array}{c}0\le m\le M-1\\ 0<n<N-1\end{array}}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\ast 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}}^{\ast }\right)}^T={\mathcal{W}}^{\ast }$$，$$G={\mathcal{W}}^{\ast }g\mathcal{W}$$，其中$$F^{\ast }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}}^{\ast }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)$$

$$\{\begin{array}{c}{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{array}$$

同理，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)\ast \frac{256}{\log \left({\mathcal{W}}_{00}+1\right)}$$可以做进一步归一化。

算法代码可见github