# 數位影像處理期末 Digital Image Process Final :::info + 課程名稱：數位影像處理 Digital Image Process + 授課老師：[name=鄭錫齊 Shyi-Chyi Cheng ] + 開課學期：1121 ::: [TOC] ## CH5 臨域處理 概念就是把指定的遮罩(奇數邊矩形或其他形狀覆蓋到影像的每個像素上 ### 空間旋積 空間濾波類似於空間旋積 旋積與濾波運算方法相同，差別在濾波器再將成將佳之前必須旋轉180度 $\sum^1_{s=-1}\sum^2_{t=-2}m(-s,-t)p(i+s,j+t)$ $\sum^1_{s=-1}\sum^2_{t=-2}m(s,t)p(i-s,j-t)$ 如果遮罩係數沒有乘以180度旋轉，則叫就空間相關性(correlation)運算(即空間濾波) ### 運算符號 平均濾波的運算結果可以寫成 $\cfrac{1}{9} a+\cfrac{1}{9} b+\cfrac{1}{9} c+\cfrac{1}{9} d+\cfrac{1}{9} e+\cfrac{1}{9} f+\cfrac{1}{9} g+\cfrac{1}{9} h+\cfrac{1}{9} i$ 該濾波可以用矩陣表示$\begin{bmatrix} \cfrac{1}{9} & \cfrac{1}{9} & \cfrac{1}{9}\\ \cfrac{1}{9} & \cfrac{1}{9} & \cfrac{1}{9}\\ \cfrac{1}{9} & \cfrac{1}{9} & \cfrac{1}{9}\\ \end{bmatrix}$ 範例$\begin{bmatrix} 1 & -2 & 1\\ -2 & 4 & -2\\ 1 & -2 & 1\\ \end{bmatrix}$運算過程 $a-2b+x-2d+4e-2f+g-2h+i$ ### 影像邊緣 * 忽略邊緣（ignore the edges） * 補零（pad with zeros） * 重複（repeating） ![image](https://hackmd.io/_uploads/B1sz2zDDT.png) * 鏡射（mirroring） ![image](https://hackmd.io/_uploads/B1ZN2Gwwp.png) ### 高通濾波器：Laplacian濾波器 係數和為0且奇數為係數和與偶數位係數和為0的卷積 ### 頻率：低通與高通濾波器 數值超出0到255的範圍 * 將負值變為正值（make negative values positive ） * 數值裁剪（clip values） $y=\begin{cases} 0 \quad if \quad x < 0 \\ x \quad if \quad 0\leq x\leq 255\\ 255 \quad if \quad 255<x \end{cases}$ * 比例轉換（scaling transformation） ![image](https://hackmd.io/_uploads/SyrkQmPva.png) #### Python ```python=! import numpy as np import scipy.ndimage as ndi x = [[170, 240, 10, 80, 150], [230, 50, 70, 140, 160], [40, 60, 130, 200, 220], [100, 200, 190, 210, 30], [110, 180, 250, 20, 90]] x = np.array(x, dtype=np.uint8) mask = np.ones((3,3)) /9 result = ndi.convolve(x,mask,mode='constant') print(result) ``` ```python=! from skimage import io import scipy.ndimage as ndi img = io.imread('cameraman.tif') mask = np.ones((9,9)) /81 result = ndi.convolve(img,mask,mode='constant') io.imshow(result) ``` ```python=! from skimage import io import scipy.ndimage as ndi img = io.imread('cameraman.tif') result = ndi.uniform_filter(img,[9,9],mode='constant') io.imshow(result) ``` ### 高斯濾波器 $f(x)=e^{(-\cfrac{x^2}{2\sigma ^2})}$ ### 邊緣銳利化 ![image](https://hackmd.io/_uploads/Hy67E7vwp.png) #### Python ```python=! k = io.imread('koala.tif') #低通濾波 u = np.array([[-2.0,-2.0,-2.0],[-2.0,25.0,-2.0],[-2.0,-2.0,-2.0]]) /9 #濾波 result = ndi.convolve(k.astype(float),u) ``` ### 高增幅濾波 與去銳利化遮罩濾波類似的方法是高增幅（high-boost）濾波器 `高增幅影像=A(原始影像)-(低通濾波影像)` * 其中A為增幅系數 * 若A = 1，則高增幅濾波器就成為典型的高通濾波器。 ### 非線性濾波器 * 最大濾波器（maximum filter） * 最小濾波器（minimum filter） * 中值濾波器(median filter) * 排序濾波器(rank order filter) * 均方根濾波器(mean squared filter) $\cfrac{\sqrt{\sum X ^ 2}}{N}$ ### 邊緣保持模糊濾波器(Eedg-preserving blurring filter) 中值濾波器 #### Python ```python=! cm=ndi.median_filter(c,size=[3,3]) ``` ### Kawahara 濾波器 * 針對卷積核内部的所有像素点二次進行分割成四個區域 * 找出四個區域當中數值最平穩的區域，求出這個區域内的平均值 * 以這個平均值作為卷積核中心像素点的數值 ![image](https://hackmd.io/_uploads/rJTsvQwPa.png) ### 雙邊濾波器(Bilateral filter) 結合兩個函數～一個函數是由幾何空間距離決定濾波器係數，另一個由像素色差決定濾波器係數 $I_p=\cfrac{1}{W_p}\sum_{q\in S}G_{\sigma_s}(||p-q||)G_{\sigma_r}(|I_p-I_1|)I_q$ $W_p=\sum_{q\in S}G_{\sigma_s}(||p-q||)G_{\sigma_r}(|I_p-I_1|)I_q$ 其中 $p$：目標像素 $q$：目標像素周圍的一個像素 $I_p$：目標像素的顏色 $I_q$：目標像素周圍的一個像素的顏色 $S$：目標像素周圍的像素群 $G_s$：與標準偏差的高斯濾波氣(加權的像素根據距離) $G_r$：與標準偏差的高斯濾波氣(加權的像素根據像素色差) ## CH6 影像幾何 ### 數據(線性)內插法 除了第一點與最後一點，$x'_i$ 完全不會與原始的$x_j$ 相對應。 必須以已知的鄰近$f(x_j)$ 值來估算函數值$f(x'_i)$。 這種以周圍數值估算函數值的方法稱為內插法( interpolation） ![image](https://hackmd.io/_uploads/BkQaM2Uva.png) 解多項式$\begin{cases} 1 = a + b \\ 4 = 8a + b \end{cases}$ 得 $x=\cfrac{3}{7}x' +\cfrac{4}{7}$ 通式$F(新數值) = \lambda f(x_2)+(1-\lambda)f(x_1)$ #### 雙線性內插法 ![image](https://hackmd.io/_uploads/H1jP4h8Pp.png) $\begin{cases} f(x, y') = \mu f(x, y + 1)+(1 - \mu)f(x,y)\\ f(x', y') = \lambda f(x + 1, y') + (1 - \lambda)f(x,y')\\ f(x + 1, y') = \mu f(x + 1, y + 1) + (1 - \mu)f(x+1, y) \end{cases}$ ### 鄰近內插法 拿原始數據最靠近新位置的地方內插 ![image](https://hackmd.io/_uploads/ryei9hLPa.png) ### 立方內插 利用已知的四個點進行內插計算 四個已知點推出三次函數進行內插 $x_0, x_1, x_2, x_3$分別對應到內插函數$f(x'_i)$ $f(x'_0) = x_0、f(x'_1) = x_1、f(x'_2) = x_2、f(x'_3) = x_3$ 求出內插函式$f(x')$後代入要找的點 #### Python 使用skimage的Transform模組中的Rescale函數 ```Python=! from skimage import transform transform.rescale(A, k, order) #order = 0 鄰近內插 #order = 1 雙線性內插 #order = 2 立方內插 ``` ### 空間濾波器放大 $m = \begin{bmatrix} 1 & 2 \\ 3 & 4\\ \end{bmatrix}$ $m_{2}[i, j]=\begin{cases}m[i//2, j//2] \\ 0\end{cases} = \begin{bmatrix} 1 &0& 2&0 \\ 0&0&0&0&\\ 3 &0&4&0\\ 0&0&0&0&\\ \end{bmatrix}$ 再對$m_2$進行內插 #### Python ```python=! hz = np.zeros((2*r,2*c)) # 將原始圖像放入新的矩陣中 hz[::2,::2] = head # 定義內插核 ne = np.array([[0,0,0],[0,1,1],[0,1,1]]) bi = np.array([[1,2,1],[2,4,2],[1,2,1]])/4.0 bc = np.array([[1,4,6,4,1],[4,16,24,16,4],[6,24,36,24,6],[4,16,24,16,4],[1,4,6,4,1]])/64.0 # 進行內插運算 hz_ne = ndi.correlate(hz,ne) hz_bi = ndi.correlate(hz,bi) hz_bc = ndi.correlate(hz,bc) ``` ### 縮小 * 又稱影像最小化 * 次取樣 ### 旋轉 $r = \begin{bmatrix} cos\theta & -sin\theta \\ sin\theta & cos\\ \end{bmatrix}$ ![image](https://hackmd.io/_uploads/B1L7haUvT.png) 旋轉後部分點會脫離座標 找到旋轉後還在座標上的點 旋轉回原本圖形做內插 #### Python ```python=! c_ne=transform.rotate(c,60, order=0) c_bi=transform.rotate(c,60, order=1) c_bc=transform.rotate(c,60, order=2) ``` ### 透視變形矯正 透視變形是指一般在拍攝影片、照片的時候，你可以看到一個物體及其周圍區域與標準鏡頭中看到的相比完全不同 ![image](https://hackmd.io/_uploads/BJQu6TUDT.png) ![image](https://hackmd.io/_uploads/BJcqaaUvp.png) $x=ax'+by'+cx'y'+d$ $y=ex'+fy'+gx'y'+h$ 代入$x,y,x',y'$各四個解8條聯立 #### Python ```python=! # 使用 cv2.getPerspectiveTransform() 獲取變換矩陣 M M = cv2.getPerspectiveTransform(src, dst) # 使用 cv2.warpPerspective() 將圖像扭曲為頂視圖 warped = cv2.warpPerspective(img, M, (w, h), flags=cv2.INTER_LINEAR) ``` ## CH7 傅立葉轉換 ### 傅立葉轉換 $f(x)=\sum^\infty_{-\infty}F(\omega)e^{i\omega x}d\omega$ $F(\omega)=\cfrac{1}{2\pi}\sum^\infty_{-\infty}f(x)e^{-i\omega x}dx$ ### 一維離散傅立葉轉換 $f=[f_0, f_1, f_2,...,f_{N-1}]$ $F=[F_0, F_1, F_2,...,F_{N-1}]$ $F_u=\cfrac{1}{N}\sum^{N-1}_{x=0}exp[-2\pi i\cfrac{xu}{N}]f_x$ 矩陣乘積表示 $\mathbf{F}=\boldsymbol{Ff}$ 其中$\mathbf{F}$是一個$N\times N$的矩陣 $\boldsymbol{F}_{m,n}=\cfrac{1}{N}exp[-2\pi i\cfrac{mn}{N}] \rightarrow \boldsymbol{F}_{m,n}=\cfrac{1}{N}\omega^{mn}$ N給定後可以定義$\omega=exp[\cfrac{-2\pi i}{N}]$ 假設$f=[1,2,3,4]$，因此$N=4$ $\omega=exp[\cfrac{-2\pi i}{4}]=exp[-\cfrac{\pi i}{2}]=\mathbf{cos}(-\cfrac{\pi}{2})+i\mathbf{sin}(-\cfrac{\pi}{2})=-i$ ### 反一維離散傅立葉轉換 $f_x=\sum^{N-1}_{u=0}exp[2\pi i\cfrac{xu}{N}]F_u$ 與DFT差別 1. 沒有縮放係數1/N 2. 指數函數中的符號改為正號 3. 總和索引變數為u，而非x ### 一維離散傅立葉轉換特性 * 線性 * 由DFT矩陣乘積可推出 * 假設$f$和$g$是相同長度的兩個向量，$p$和$q$是純量 另$h=pf+qg$ 若$F、G、H$分別為$f、g、h$的DFT 則$H=pF+qG$ * 平移 * 將向量x的各個元素$x_n$乘以(-1)^n，也就是每兩個變號 另新的向量為$x'$，$x'$的DFT$x'$若將左有兩邊互換，就和x的DFTx相同 **白話文：把一段頻率相反後結果相同** * 縮放 * $\boldsymbol{F}[f(kx)] = (1/k)\boldsymbol{F}(\omega/k)$ * 共軛對稱$F_k=\bar{F}_{N-k}$ * 捲積$z=x\otimes y, z_y=\cfrac{1}{N}\sum^{N-1}_{n=0}x_ny_{k-n}$ 要使用空間濾波器S對影像M進行旋轉計算 * 快速傅立葉轉換 ### 二維離散傅立葉轉換特性 * 相似性 * DFT作為空間濾波器 * 分離性 * 線性 * DC係數 $F(0,0)=\sum^{M-1}_{x=0}\sum^{N-1}_{y=0}f(x,y)$ * 捲積定理：$F(M\times S)=F(M)\cdot F(S')$ * 平移 所有元素f(x,y)以(-1)^{x+y}處理 * 共軛對稱 $F(u,v)=F^*(-u+pM,-v+qN)$ #### Python ```python=! cf = fftshift(fft2(c)) #傅立葉轉換後移位 d= 15 ar = np.arange(-128,128) x,y = np.meshgrid(ar,ar) z= np.sqrt(x**2+y**2) #低於d的設成0，其餘設成1 f= (z>15)*1 #高於15的設為1，低於15設成0 #大於改成小於為低通濾波 cf1 = cf*f #套用至cf上 ifft2(cf1) #反傅立葉轉換 ``` ### Butterworth 濾波函數 理想濾波器直接切除傅利葉轉換距中心某個距離外的部分。 這種截頻點的使用十分方便，但缺點是結果會產生不必要的瑕疵（波紋）。 要避免這種現象，可使用截頻點較不銳利的圓形當作濾波矩陣。 #### Python ```python=! cf = fftshift(fft2(c)) #傅立葉轉換後移位 d = 15 ar = np.arange(-128, 128) x,y = np.meshgrid(ar, ar) z = np.sqrt(x ** 2 + y ** 2) #低於d的設成0，其餘設成1 bh = 1 - 1.0 / (1.0 + ((x ** 2 + y ** 2) / d ** 2) ** 2) #高於15的設為1，低於15設成0 #去除1-為低通濾波 cfbh = cf*bh #套用至cf上 ifft2(cfbh) #反傅立葉轉換 ``` ### 高斯濾波 ```python=! cf = fftshift(fft2(c)) ar = np.arange(-128, 128) x,y = np.meshgrid(ar, ar) sigma = 30.0 gh = 1.0 - np.exp(-(x ** 2 + y ** 2) / sigma ** 2) #高於30的為1，低於30的設為0 #去除1-為低通濾波 cfgh = cf * gh ``` ### 同態濾波 $f(x,y)=i(x,y)r(x,y)$ $i(x,y)$為照明(illumination) $r(x,y)$為反射(reflectance) ## CH8 影像復原 ### 影像品質劣化模型 $g(x,y)=f(x,y)*h(x,y)$ 若影像發生隨機的錯誤，劣化的影像為 $g(x,y)=f(x,y)*h(x,y)+n(x,y)$ 也可以在頻率做劣化 $G(x,y)=G(x,y)G(x,y)+N(x,y)$ ### 雜訊 * 因為外來干擾而造成任何形式的影像訊號劣化。 * 錯誤顯示在輸出影像上的形式，會依照訊號所受干擾的類型而有所不同。 * 一般都可以分析出會發生何種錯誤 * 影像會出現哪一種類型的雜訊，可以選擇最適合的方法來去除雜訊效果。 ### 鹽和胡椒雜訊 * 又稱為脈衝雜訊（impulse noise）、散粒雜訊（shot noise）或二元雜訊（binary noise）。 * 鹽和胡椒雜訊（Salt and Pepper Noise）可能是因為影像訊號受到突如其來的強烈干擾而產 * 呈現方式是整個影像任意散布黑色或白色（或兩者皆有）的像素。 ```python=! gn = util.noise.random_noise(b,mode = 's&p') ``` ### 高斯雜訊 * 高斯雜訊（Gaussian noise）是白雜訊（white noise）的理想形態，是由訊號中隨機的擾動而造成。 * 若設影像為I，高斯雜訊為N，兩者相加便可得到有雜訊的影像 ```python=! gn = util.noise.random_noise(b,mode = 'gaussian') ``` ### 斑點雜訊 * 像素值乘上亂數數值，則會產生斑點（speckle）雜訊（或簡稱斑點），因此又稱為乘積雜訊（multiplicative noise）。 ```python=! gn = util.noise.random_noise(b,mode = 'speckle') ``` ### 週期性雜訊 ```python=! x,y = np.mgrid[0:r,0:c].astype('float32') p = np.sin(x / 3 + y / 3) + 1.0 gp = (2 * util.img_as_float(b) + p / 2) / 3 ``` ### 去除鹽和胡椒雜訊 * 低通濾波器：平均濾波器 ```python=! # bn為經雜訊處理後的影像 brn = ndi.uniform_filter(bn,3) ``` * 中位數濾波器 ```python=! brn = ndi.median_filter(bn,3) ``` ### 排序濾波 中位數濾波其實是排序濾波（rank-order filtering）的一種特殊形式。 ```python=! cross = np.array([[0,1,0],[1,1,1],[0,1,0]]) brn = ndi.median_filter(bn,footprint = ross) ``` ### 歧異點方法 * 使用中位數濾波器進行運算，基本上速度頗慢：計算每個像素值都至少得排序9 個數值 * 較快去除雜訊的方法如下： 1. 選擇閥值D。 2. 比較指定像素值p 與周圍8 個像素的平均數值m。 3. 若$| p – m |> D$，則認定該像素為雜訊；反之則否。 4. 若像素為雜訊，則用m 取代；反之，則保留原數值。 ```python=! ave = np.array([[1,1,1],[1,0,1],[1,1,1]])/8.0 bnave = ndi.convolve(bn,ave) D = 0.5 r = (abs(bn - bnave) > D) * 1.0 brn = r * bnave + (1 - r) * bn ``` ### 去除高斯雜訊 * 影像平均 假設共有100 張影像，每一張都有雜訊 因為$N_i$ 為常態分布且平均值為0，所以所有$N_i$ 的平均數可能相當接近0，$N_i$數目越大，則越接近0。 白話文：因為雜訊平均分布，所以加一堆雜訊讓他抵銷 ```python=! for i in range(num_img): bn[:,:,i] = util.noise.random_noise(b,mode = 'gaussian') brn = np.mean(bn,axis = 2) ``` * 平均濾波 ```python=! brn = ndi.uniform_filter(brn,3) ``` ### 可適性濾波 * 可適性濾波器（adaptive filters）是會隨著遮罩下的灰階值改變其特性的濾波器。 * 最小均方誤差濾波器（minimum mean-square error filter） $m_f+\cfrac{\sigma_f^2}{\sigma_f^2+\sigma_g^2}(g-m_f)$ 遮罩中的平均值不一定為0。 * Wiener 濾波器（Wiener filters） $m_f+\cfrac{max{0,\sigma_f^2-n}}{max{\sigma_f^2,n}}(g-m_f)$ ```python=! from scipy.signal import wiener brn = wiener(bn,[3,3]) ``` ### 去除週期性雜訊 * 帶拒濾波（band reject filtering） 用傅立葉轉換去除頻率 ```python=! gpf = fftshift(fft2(gp)) z = np.sqrt((x - 128) ** 2 + (y - 128) ** 2) d = np.sqrt(300) k = 5 br = (z < np.floor(d - k)) | (z > np.ceil(d + k)) pr = np.abs(ifft2(gpf * br)) ``` * 陷波濾波（notch filtering） ### 反轉濾波 $X(i, j) = \cfrac{Y(i,j)}{F(i,j)}L(i,j)$ ```python=! bworth = 1./(1+(np.sqrt(2)-1)*((x**2+y**2) /15**2)**2) hf = fftshift(fft2(h)) hw=hf*bworth hwa = np.abs(ifft2(hw)) ``` ### 除式設限反轉濾波 ```python=! bworth = 1./(1+(np.sqrt(2)-1)*((x**2+y**2) /15**2)**2) bw= np.copy(bworth) bw[np.where(bw<d)]=1 fbw= fftshift(fft2(blur))/bw ba = np.abs(ifft2(fbw)) ``` ### 去除動態模糊 要去除影像的模糊現象，必須將其傅利葉轉換結果除以模糊濾波器的傅利葉轉換 也就是先要產生一個模糊濾波器的傅利葉轉換矩陣 * 直接除法去除 ```python=! m= np.ones((1,7))/7 cm = ndi.correlate(cr,m) m2 = np.zeros((256,256)) m2[0,0:7]=m mf = fft2(m2) bmi = np.abs(ifft2(fft2(cm)/mf)) ``` * 對除法設限去除動態模糊 ```python=! m= np.ones((1,7))/7 cm = ndi.correlate(cr,m) m2 = np.zeros((256,256)) m2[0,0:7]=m mf = fft2(m2) d=0.2 mf[np.where(np.abs(mf)<d)]=1 bmi = np.abs(ifft2(fft2(cm)/mf)) ```