一些基本数字图像处理算法
一些基本數字圖像處理算法
版權聲明:本文為原創文章,未經博主允許不得用于商業用途。
所有的圖像算法都在DIPAlgorithm類中,并且所有算法都為抽象成員函數。我已經按照java注釋規范為所有方法添加使用說明注釋,具體實現可見于DIPAlgorithm.java,這里只做算法說明。
1 圖像扭曲
模仿PS的扭曲功能,通過建立一個三角形映射網格實現對圖像的扭曲。
如上圖,一共設置了45個控制點圍成74個三角形網格
扭曲即形變處理其實是尋找一個函數,以所有網格頂點原始坐標為輸入,扭曲后所有網格頂點坐標為輸出。為了簡化計算任務,采用控制柵格插值方法,對每個三角網格獨立計算映射關系,如下圖:
即求解矩陣MMM滿足MA=BMA = BMA=B,其中AAA為原頂點的齊次矩陣:
A=[x1y11x2y21x3y31]A = \begin{bmatrix} x_{1} & y_{1} & 1 \\ x_{2} & y_{2} & 1 \\ x_{3} & y_{3} & 1 \\ \end{bmatrix} A=???x1?x2?x3??y1?y2?y3??111????
B為形變后頂點的其次矩陣:
B=[x1′x2′x3′y1′y2′y3′]B = \lbrack\begin{matrix} x_{1}^{'} & x_{2}^{'} & x_{3}^{'} \\ y_{1}^{'} & y_{2}^{'} & y_{3}^{'} \\ \end{matrix}\rbrack B=[x1′?y1′??x2′?y2′??x3′?y3′??]
M即為2×32 \times 32×3的映射矩陣,且由于三角形三點不共線,因此A為可逆陣,
M=BA?1M = BA^{- 1} M=BA?1
對于三角形中的點p(x,y)p\left( x,\ y \right)p(x,?y),其映射后坐標p′=M[xy1]p^{'} = M\begin{bmatrix} x \\ y \\ 1 \\ \end{bmatrix}p′=M???xy1????
2 直方圖計算
直方圖計算實際上即求圖像的概率密度函數PDF,只需遍歷一次所有像素點即可獲得。
3 直方圖均衡化算法
對于連續圖像直方圖均衡化其實是種點運算f,
對不同灰度值做映射,使得所有像素頻率相等。
對于點運算f,有如下性質:
DB=f(DA),HB(DB)ΔDB=HA(DA)ΔDAD_{B} = f\left( D_{A} \right),\ H_{B}\left( D_{B} \right)\Delta D_{B} = H_{A}\left( D_{A} \right)\Delta D_{A} DB?=f(DA?),?HB?(DB?)ΔDB?=HA?(DA?)ΔDA?
其中D為灰度值,H即為灰度值在圖像中的頻數,整理可得
HB(DB)=HA(DA)ΔDAΔDB=HA(DA)ΔDBΔDA=HA(DA)dDBdDAH_{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}}} HB?(DB?)=ΔDB?HA?(DA?)ΔDA??=ΔDA?ΔDB??HA?(DA?)?=dDA?dDB??HA?(DA?)?
=HA(DA)f′(DA)=HA(f?1(DB))f′(f?1(DB))= \frac{H_{A}\left( D_{A} \right)}{f'(D_{A})} = \frac{H_{A}\left( f^{- 1}\left( D_{B} \right) \right)}{f'(f^{- 1}(D_{B}))} =f′(DA?)HA?(DA?)?=f′(f?1(DB?))HA?(f?1(DB?))?
即:
尋找函數f使得HB(D)H_{B}(D)HB?(D)為常數A0Dm,A0,Dm\frac{A_{0}}{D_{m}},A_{0},D_{m}Dm?A0??,A0?,Dm?。
由(1)可知,KaTeX parse error: Undefined control sequence: \ at position 58: …\right)}{f'(D)}\? ?\Rightarrow f^{…
即f(D)=DmCDF(D)f\left( D \right) = D_{m}CDF(D)f(D)=Dm?CDF(D),CDF即累積分布函數
因此只需求得直方圖的前序和即可獲得映射關系。
4 圖像灰度化
目前比較符合人眼的灰度化權重為0.299、0.578和0.114,為了加速計算使用近似公式D=(3r+g+6b)/10D = (3r + g + 6b)/10D=(3r+g+6b)/10
5 圖像二值化
我使用的二值化算法為OSTU大律二值化算法。二值化操作即利用分割閾值u,將圖片分為前景后景兩部分。OSTU大律法認為使得前景像素和背景像素灰度方差g最大的閾值即為最佳分割閾值。
g=w0w1(u0?u1)2g = w_{0}w_{1}\left( u_{0} - u_{1} \right)^{2} g=w0?w1?(u0??u1?)2
其中w0,w1w_{0},\ w_{1}w0?,?w1?為前景、后景在圖像中的比例,KaTeX parse error: Undefined control sequence: \ at position 7: u_{0},\? ?u_{1}為前景、后景的平均灰度。
在實現時,只需遍歷所有灰度,利用CDF求出每種灰度的方差,取最大者作為閾值即可。
6 前景分離
目前主流的前景分離為深度學習算法。這里只使用了最基本的閾值分離法,分別為RGB三個通道設置不同閾值,將小于閾值的像素作為背景,大于閾值的作為前景。
7 濾波
我使用的濾波方法是高斯濾波和中值濾波,高斯濾波即使用二維高斯函數作為濾波函數,中值濾波即使用鄰域的中位數作為濾波函數。
高斯濾波器為線性濾波器,可以有效消除高斯噪聲。由于高斯函數離中值越近權重越大,因此相對于均值濾波器更加柔和,對邊緣的保留效果更好。這里我使用的是如下矩陣做卷積:
[1232124642367632464212321]\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} ???????12321?24642?36763?24642?12321????????
中值濾波器為非線性濾波器,可以有效的去除椒鹽噪聲和斑點噪聲并且不會使圖像變模糊。
8 形態學擴張和腐蝕
形態學腐蝕可記為AΘB\text{AΘB}AΘB,其中A為輸入圖像,B為結構單元。對于二值圖像,當且僅當當前像素點滿足腐結構單元時才會被保留。對于灰度圖像,則可類比為最小值,即
fΘb(x,y)=min{f(x?x′,y?y′)?b(x′,y′)∣(x′,y′∈Db)}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?BA\bigoplus BA?B,對于二值圖像,將每個有效像素點的鄰域結構單元置1,對于灰度圖像則取最大值,即
f?b(x,y)=max{f(x?x′,y?y′)?b(x′,y′)∣(x′,y′∈Db)}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=(AΘB)?BAoB = \left( \text{AΘB} \right)\bigoplus BAoB=(AΘB)?B)和結構閉運算(AoB=(A?B)ΘBAoB = \left( A\bigoplus B \right)\text{ΘB}AoB=(A?B)ΘB)對圖像進行粗化、細化、濾波等處理
9 傅里葉變換和濾波
變換公式
傅里葉變換可以將信號從時域轉換到頻域,因此可以看出許多時域中不明顯的特征。二維傅里葉變換(CFT)公式如下:
F(u,v)=?f(x,y)e?2πj→(ux+vy)dxdyF\left( u,v \right) = \iint_{}^{}{f\left( x,y \right)e^{- 2\pi\overrightarrow{j}(ux + vy)}}\text{dxdy} F(u,v)=??f(x,y)e?2πj?(ux+vy)dxdy
其中j→2=?1,f,F{\overrightarrow{j}}^{2} = - 1,f,Fj?2=?1,f,F,同樣二維傅里葉逆變換公式如下:
f(x,y)=?F(u,v)e2πj→(ux+vy)dudvf\left( x,y \right) = \iint_{}^{}{F\left( u,v \right)e^{2\pi\overrightarrow{j}(ux + vy)}}\text{dudv} f(x,y)=??F(u,v)e2πj?(ux+vy)dudv
對于離散函數,可以定義離散二維傅里葉變換(DFT)和逆變換:
G(m,n)=1MN∑0≤i≤M?10<k<N?1g(i,k)e?2πj→(imM+jnN)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(m,n)=MN?1?0≤?i?≤?M?10<k<N?1??∑?g(i,k)e?2πj?(Mim?+Njn?)
g(i,k)=1MN∑0≤m≤M?10<n<N?1g(m,n)e2πj→(imM+jnN)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})}} g(i,k)=MN?1?0≤?m?≤?M?10<n<N?1??∑?g(m,n)e2πj?(Mim?+Njn?)
DFT可以理解為對連續二維信號進行了頻率為M,
N的采樣,之后通過計算其和頻域空間M*N個基向量的相關性(在該方向投影)將時域信號映射到頻域。iDFT可以理解為通過M*N個基向量合成原始時域信號。
矩陣表示
傅里葉變換實際上是一種線性變換,因此在實際計算中常常將ggg擴充為N?NN*NN?N方陣,此時DFT可以通過矩陣表示:G=W?1gW,Wik=1Ne2πj→ikNG = \mathcal{W}^{- 1}g\mathcal{W},\mathcal{W}_{\text{ik}} = \frac{1}{N}e^{2\pi\overrightarrow{j}\frac{\text{ik}}{N}}G=W?1gW,Wik?=N1?e2πj?Nik?。
易知Wik=Wki\mathcal{W}_{\text{ik}} = \mathcal{W}_{\text{ki}}Wik?=Wki?,且為正交矩陣,因此W\mathcal{W}W為酉矩陣,即W?1=(W?)T=W?\mathcal{W}^{- 1} = \left( \mathcal{W}^{*} \right)^{T} = \mathcal{W}^{*}W?1=(W?)T=W?,G=W?gWG = \mathcal{W}^{*}g\mathcal{W}G=W?gW,其中F?FF^{*}FF?F。
由于傅里葉變換為酉變換,即Wt=W?1\mathcal{W}^{t} = \mathcal{W}^{- 1}Wt=W?1
圖像的傅里葉變換
對于二維圖片可以看作二維矩陣,因此可以進行DFT。二維圖片經過DFT后獲得的復矩陣的模矩陣可以表示每個頻率信號的強度(也可看作先做自相關后再進行傅里葉變換),經過適當處理即可轉化為灰度能量譜圖片。
線性噪聲在頻域中通常為點或線,因此可以通過傅里葉變換后進行濾波再通過逆變換復原圖片。
算法實現
在實際實現時,根據歐拉公式,e?j→t=cost?j→sint,ej→t=cost+j→sinte^{- \overrightarrow{j}t} = cost - \overrightarrow{j}\text{sint},\ e^{\overrightarrow{j}t} = cost + \overrightarrow{j}\text{sint}e?j?t=cost?j?sint,?ej?t=cost+j?sint,因此傅里葉變換的核矩陣可以表示為Wik=cos?(2πik)?j→sin?(2πik)N\mathcal{W}_{\text{ik}} = \frac{\cos\left( 2\pi ik \right) - \overrightarrow{j}\sin\left( 2\pi ik \right)}{N}Wik?=Ncos(2πik)?j?sin(2πik)?,為方便運算將W\mathcal{W}W分解為虛部系數Wlm\mathcal{W}_{\text{lm}}Wlm?和實部系數Wre\mathcal{W}_{\text{re}}Wre?,其中則W=Wre+j→Wlm\mathcal{W} = \mathcal{W}_{\text{re}} + \overrightarrow{j}\mathcal{W}_{\text{lm}}W=Wre?+j?Wlm?。變換結果同樣分解為G=Gre+j→GlmG = G_{\text{re}} + \overrightarrow{j}G_{\text{lm}}G=Gre?+j?Glm?,則DFT可以表示為:
G=W?gW=(Wre?j→Wlm)g(Wre+j→Wlm)=WregWre+WlmgWlm?j→(WlmgWre+WregWlm)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) G=W?gW=(Wre??j?Wlm?)g(Wre?+j?Wlm?)=Wre?gWre?+Wlm?gWlm??j?(Wlm?gWre?+Wre?gWlm?)
{Gre=WregWre+WlmgWlmGlm=?WlmgWre?WregWlm\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.\ {Gre?=Wre?gWre?+Wlm?gWlm?Glm?=?Wlm?gWre??Wre?gWlm???
同理,iDFT可以表示為:
g=(Wre+j→Wlm)(Gre+j→Glm)(Wre?j→Wlm)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) g=(Wre?+j?Wlm?)(Gre?+j?Glm?)(Wre??j?Wlm?)
其中,為了將能量譜轉化為可見的灰度圖,為能量譜取對數值進行歸一化。且由于在頻域中兩個維度頻率都為0時(即W00\mathcal{W}_{00}W00?處)為圖像能量的總和,因此通過log(e+1)?256log?(W00+1)log(e + 1)*\frac{256}{\log\left( \mathcal{W}_{00} + 1 \right)}log(e+1)?log(W00?+1)256?可以做進一步歸一化。
算法代碼可見github
總結
以上是生活随笔為你收集整理的一些基本数字图像处理算法的全部內容,希望文章能夠幫你解決所遇到的問題。
- 上一篇: 《人工智能-一种现代的方法》阅读笔记
- 下一篇: 只有一百行的xss扫描工具——DSXS源