IA10 - 機率方法在圖像去噪之應用Part 2: Color image

一張帶有高斯噪聲的Lena相片

這個專題將延續IA8的內容，但我們將跳脫黑白兩色並探討彩色照片的去噪方法，因為色彩充斥著我們日常的生活，這個世界上僅有非常少數的人因為先後天的因素不能感知色彩，或大腦對色彩的反應發生錯誤。所以色彩，併同聲音、觸覺等皆佔了人類日常生活中非常重要的一部份。有別於IA像素的毀損，噪聲代表的是該處的像素值 (被認為) 偏離了真正的值，這可能是肇因於多種因素，在這篇文章裡我們將著重於如何重建出真正的像素值。在這篇專題中對應的notebook和檔案皆儲存在IA10這個資料夾內請將其全部下載下來。

一、去噪: A naïve way

1.1 來自IA4的啟發

以一張灰階的圖片為例，我們可以將這張圖片表示為一個二維的陣列\(\mathbf{P}\)： \[ \mathbf{P}=\left(\begin{array}{cccc} p_{11} & p_{12} & \cdots & p_{1W}\\ p_{21} & p_{22}\\ \vdots & & \ddots\\ p_{H1} & & & p_{HW} \end{array}\right) \] 其中W和H表示圖片的寬和高，而\(p_{ij}\)表示在位置\((i,j)\)處的像素值，對灰階圖片而言則表示強度值，通常介於0~255的整數，因為共有\(2^8=256\)階，又稱8-bit或256色。補充一點，上式中的\(\mathbf{P}\)是灰階照片，所以陣列對應的shape是(H,W)，如果是彩色照片則含有3通道，對應的shape則變為(H,W,3)，最後一個3表示RGB三通道，而在python中的OpenCV套件cv2所讀入的彩色照片三通道的順序是BGR，如何轉換到matplotlib.pyplot的imshow可參考IA7或notebook的內容。之後為了方便我們將\(p_{ij}\)表示成\(p_i\)，表示在第\(i\)處的像素值，這是因為在正統圖像處理領域時常會將一個二維的圖像陣列重新reshape成\(WH\times 1\)的row vector，而這種排序又稱為lexicographical order，不過這裡我不探討這些技術細節。
　原則上\(p_i\)可能都含有來自噪聲的貢獻，假定噪聲的值是\(n_i\)，真正的像素值是\(x_i\)，則有 \[ p_i = x_i + n_i. \tag{1} \] 通常噪聲\(n_i\)可以視為是一個均值為0的擾動。假定我們有一台老舊的數位相機，其CCD和電路板由於劣化的關係會使得拍出來的相片產生雜訊，假設我們拍攝\(K\)張相同的照片，由於來自相同的場景，所以Eq. (1)中的\(x_i\)在這\(K\)張相片中應該都要一樣，我們會在每一張相片中量到不同的\(p_i\)是由於拍攝時的生成噪聲\(n_i\)的隨機性所導致的。所以從這\(K\)張相片中，我們可以推測在\(i\)處像素的期望值為 \[ \mathbb{E}(p_i) = \mathbb{E}(x_i + n_i) = \mathbb{E}(x_i) + \mathbb{E}(n_i).\tag{2} \] 而在上式中由於噪聲的均值為0，所以有\(\mathbb{E}(n_i)=0\)，這時候我們便將\(\mathbb{E}(p_i)\)表為這\(K\)張照片在\(i\)處所量到像素值\(p_i\)的平均，即 \[ \mathbb{E}(x_i)=\bar{p}_i =\frac{1}{K} \sum_{k=1}^K p_i^k. \tag{3} \] 透過Eq. (3)我們推論出了\(i\)處真實像素值\(x_i\)可以透過多張同場景拍攝的照片的像素值取平均\(\bar{p}_i\)而來，當照片愈多張的時候，\(\bar{p}_i\)愈接近\(x_i\)，這可方法和我們在IA4一次傳送多個一樣的訊號，最後在接收端用veto的方式來糾錯只差在我們是取平均數，而這個方法在早期天文學觀測所拍攝的相片中去噪也常常使用。

1.2 簡單的例子

讓我們來看一個簡單的例子，先用cv2套件讀入images這個資料夾下名為lena.jpeg的一張照片，並用以下方式對圖片生成高斯噪聲，其中變量scale對應的是噪聲的標準差 
# Define Gaussian noise function
def gauss_noise(im, mean = 0, sd = 10):
    """
    Function for generating Gaussian noise on an image.   
    ...
    """
    noise = np.random.normal(mean, sd, size = im.shape)
    noise_im = np.clip(im + noise, 0, 255) # chop values outside [0,255]
    return np.uint8(noise_im)
上面我們使用np.random.normal來實現隨機的高斯噪聲，np.clip是為了截斷超過\([0,255]\)之間的值，最後我們將回傳影像陣列的值轉換成8-bit unsigned integer以符合256色的規範。我們將讀入的Lena照用上面的函數生成8張雜訊相片，如圖1所示。

圖1. 8張帶有高斯噪聲的Lena照片

　那我們便可以依照Eq. (3)的方式將所有的照片在各個通到的像素值全部相加在取算術平均值，所以真正乾淨的圖片陣列\(\mathbf{X}\sim\sum_{k=1}^K{\mathbf{P}}/K\)，結果如圖2所示。不諱言的，去噪的結果其實蠻不錯的，可以看出比起圖1來，圖2的Lena乾淨非常多！

圖2. 透過Eq. (3)去噪後的Lena照片

　但有一個問題就是，圖2相比圖1而言雖然“看起來”好很多，但是這個看起來每個人有不同的基準，有人覺得這樣就很好，有人覺得差強人意，我們需要一個量化上的數字來判別通過Eq. (3)的方式去噪後，圖片劣化的程度確實有降低。一個常用的方法便是比較峰值訊噪比PSNR (peak signal-to-noise ratio)，其定義為 \[ {\rm PSNR}=10\log_{10} \frac{{\rm MAX}^2}{{\rm MSE}},\tag{4} \] 其中MAX指的是像素值的最大值，在這裡是255，而MSE則是目標圖片和真實圖片之間的方均根誤差。PSNR的單位是dB，值愈大代表劣化程度愈低，當然前題是我們得要有一個真實的圖片來和他做比較，而這裡不外乎的是在加上高斯噪聲前的lena.jpeg，那張圖片是我們真實且乾淨的照片，我們將會拿它來和去噪前和去噪後的圖片來做比較。實現PSNR的代碼如下 
# Define function for calculating PSNR index
def psnr(im, tr_im):
    """
    Function for calculating PSNR index
    ...
    """
    mse = (np.float32(im) - np.float32(tr_im))**2
    return 10*np.log10(255**2/np.mean(mse))
注意我們將輸入的影像陣列在計算時先轉成32-bit的浮點數，這是因為若使用uint8的話所有的值只能取在[0,255]之間的整數，超出255的部份會重新從0開始count，而小於0的部份會倒回去從255回推，這不符合一般我們數值計算的方式，所以要改成32-bit浮點數以擴大可處理的數字範圍，否則psnr函數輸出的值會出現問題。
　透過上面的代碼，我們計算出圖1中8張混有高斯噪聲的Lena相片的PSNR值都落在14.98 dB附近，而圖2去噪後PSNR則增為23.78 dB，所以Eq. (3)取多張照片算術平均像素值的去噪方法確實能有效降低雜訊並減少劣化。上述實現的過程請參照notebook的內容，另外在這裡我會時常使用雙斜線//來表示除法，這是由於8-bit色彩空間皆為整數，所以我們只取模數 (modulus) 並捨去餘數，這和Python 2的除法結果相當，當然若想用像np.uint8來強迫它變成8-bit無號整數也無不可，只是在輸出的像素值和我的方法上可能會有\(\pm 1\)的差距。
　事實上，我們並不總能一次擁有多張來自相同場景照片的奢侈，有時候我們僅僅只有唯一一張照片，所以利用Eq. (3)這種方法便不適合用在單張圖像的去噪，我們接著要討論的便是如何僅利用那僅有的一張照片來除去雜訊，不過在此之前，我們要先複習一下圖像陣列與馬可夫毯的概念。

二、複習圖像陣列與馬可夫毯的關係

在IA7的時候我們提到了圖像上的某處的像素與其鄰近的像素構成了一個馬可夫毯 (Markov blanket) 的概念，在這裡我們快速地複習一下，這裡唯一和IA7不同的是毀損處不存在像素灰度值，那裡填的是所謂的NaN，但含噪聲的像素則不同，該處是有灰度值的，但是該灰度值\(p_i\)並不是真正的值，而是已經如Eq. (1)一樣已經加上一個噪聲\(n_i\)干擾的貢獻，我們必須要重建真正的\(x_i\)。另外我們會使用灰度值和像素值來代表圖像陣列上某個位置\(i\)的像素所應填的數值，兩者在這篇專題中代表的是相同的意思，還有由於這篇專題和IA7之間的寫作時間有些間隔了，所以在符號的使用上會略有差異，況且這裡有些觀念在IA7裡也因為比較sloppy的關係而沒有提到或隨意帶過，我會在這裡較有系統性的重整這些內容，不熟這些觀念的讀者有時間的話建議不要跳過這一節。

2.1 圖像作為馬可夫毯的概述

如圖3所示，在中心的像素\(p_i\)及其鄰近的像素\(x_j\)所構成的結構為一馬可夫毯，在該圖上我們讓\(p_i\)與上下左右4個鄰居相連 (即相互影響)，所以這是一個四連通 (4-connected) 的結構，當然我們也可以讓\(p_i\)與更多的鄰居相連，構成一個多連通 (multi-connected) 的狀態，但事實上是與太遙遠的像素相關聯是沒有實際上的意義，這樣做反而會干擾並徒增計算上的困難，比如說一張照片上某處像素對應人臉的肉色，但其與遙遠處對應到天空藍色的像素值本應無關。

圖3. 圖像陣列上某個像素\(p_i\)與鄰近的像素\(x_i\)所構成的馬可夫毯

2.2 如何除去像素中的噪聲?

原則上，我們知道在\(i\)處真正的灰度值\(x_i\)和其周圍的鄰居\(x_j\)是有相關的，另為了方便起見我們有時候也會用\(\mathcal{N}\)來代表所有鄰居的集合，在這裡則是\(\mathcal{N}=\{x_1,x_2,x_3,x_4\}\)。這裡做一個假設，即\(x_i\)和其鄰居們\(\mathcal{N}\)相關聯的意義表雖然\(x_i\)和\(x_j\)之間必然存在差異，但它們間的差異應該遵從高斯分佈。基於這點，則當\(\mathcal{N}\)的狀態給定時，\(x_i\)最有可能的取值必然是使得它和\(\mathcal{N}\)之間的方均根誤差MSE最小，即 \[ \hat{x}_i={\rm argmin}~(x_i-x_j)^2. \tag{5} \] 上式中\(\hat{x}_i\)為真實像素值的估計，且由於鄰居不只一個，實際在計算時會對考慮所有的\(j\in\mathcal{N}\)，我們稍後再回來這個部分。
　另一方面如果雜號是高斯的，那麼同理我們也考慮觀測值\(p_i\)和真實值\(x_i\)之間的MSE為 \[ \hat{x}_i = {\rm argmin}~(x_i-p_i)^2.\tag{6} \] 還有計算時我們會同時考慮Eqs. (5,6)，所以trivally滿足Eq. (6)的情況\(x_i=p_i\)不見得滿足Eq. (5)，死了這條心吧！所以我們將Eqs. (5,6)整合在一起並改寫得到 \[ E(p_i,x_i,\mathcal{N}) = \Phi_i(p_i,x_i) + \sum_{j\in\mathcal{N}} \Psi_{ij} \tag{7} \] 其中\(\Phi_i\)即Eq. (6)且 \[ \Psi_{ij}=\lambda \min [(x_i-x_j)^2,\Psi_{\rm trun}]. \tag{8} \] 在Eq. (8)裡我們將所有鄰居\(\mathcal{N}\)的貢獻都收在求和符號裡了，其中超參數 (hyper parameter) \(\lambda\)和\(\Psi_{\rm trun}\)分別用來表示懲罰的程度 (penalty) 與MSE的截斷值 (truncated value)。\(\Psi_{\rm trun}\)通常會影響修復後的灰度值和周圍之間差異的大小，截斷值愈大表示\(x_{i,j}\)之間愈smooth，愈小則表示愈sharp。通過Eq. (7)我們得到了對真實像素值\(x_i\)的估計為 \[ \hat{x}_i ={\rm argmin}~E(p_i,x_i,\mathcal{N}). \tag{9} \]　在IA7裡我們稱Eq. (9)中的\(E\)為能量函數，並把它和統計力學裡決定系統狀態的配分函數關聯起來，這裡我們重新解釋這個觀點，由於在Eq. (9)中\(p_i\)和\(\mathcal{N}\)都是給定的，而\(x_i\)共有256種可能，其中每一種值出現的機率可表示成 \[ p(x_i|p_i,\mathcal{N})=\frac{e^{-E(x_i)}}{\sum_{x_i=0}^{255}e^{-E(x_i)}}.\tag{10} \] 在等式右方我們就忽略不寫出能量函數裡那些變量\(p_i\)和\(\mathcal{N}\)了，讀者只要知道它們確實是存在的就好。
　如果還記得我們在IA8裡所提到的決定性 (determistic) 與隨機性 (probabilistic) 算法的概念的話，這裡透過Eq. (9)在每一次迭代中都強迫取使能量函數\(E\)最小的那個\(x_i\)的方式就稱作決定性算法，而如果我們採用Eq. (10)的方法計算出256個灰度值每個值出現的機率，並採樣出\(x_i\)的值 (Gibbs採樣，但通常會用multinomial分佈來實現，後述) 後完成一次迭代，即使採樣出來的\(x_i\)值並不一定滿足Eq. (9)的使能量最小也沒關係，這就稱作隨機性算法，熟悉機器學習的讀者也許會發現Eq. (10)和多分類模型中使用的softmax函數如出一轍。
　在IA8裡我們見識到了隨機性算法不容易被truncated在local maxima的泥淖裡，但在下一節的實作我們會發現隨機性算法在這裡並沒有顯著的優勢，亦即修復的結果並不會優於決定性算法，而且還要花費非常多的計算資源，在第四節裡我們會探討為何如此。好的，基本上已經完成理論的複習，並重整了一些在IA7中沒有提到的內容，接著我們要動手用python來實作彩圖去噪。

三、圖像去噪實作

在本節的實作裡，我只會講解如何建立計算Eqs. (9,10)的函數，然後概要的解釋如何透過這兩個函數來分別建立決定性和隨機性算法的圖像去噪機制。事實上，我寫了兩個class (收錄在imdenoise.py裡) 來實現它們，在操作的時候我會直接讀入這兩個類來操作，讀者只要知道如何操作這兩個class來去噪即可，背後如何寫這兩個類的腳本的知識可不必知道，而前面使用到的gaussian_noise和psnr函數實際上也都收在一個叫utils.py的腳本裡，讀者會發現這個腳本內還有生成其它噪聲像是shot noise (其實就是Poisson noise) 的函數。

3.1 構造能量函數與決定性算法中的像素取值

假設我們有一張帶有噪聲的灰階圖像如圖4，用陣列表示為img，假定該陣列的大小為img.shape=(R,C)，今天我們要修復位置在img[r,c]的目標像素，假定原來帶有噪聲的像素值img[r,c]\(=p\)，且鄰著img[r,c]的九宮格代表它目標像素的鄰域，由於只考慮四連通，所以只有橘色的部分是鄰居\(\mathcal{N}\)。
　假設它真正的不含噪聲的像素值是\(x\)，而\(x\)共有256種可能，我可以令x=np.arange(256)，所以Eq. (9)的能量函數就可以表示成 
def energy(p, nei, lam, cutoff, x = np.arange(256)):
    """
    Calculating the energies of 256 pixel values for
    the center pixel in the box, Eq. (7)
    """
    nei = np.asarray(nei)
    Phi = (x-p)**2
    Psi = lam*(np.clip((x-nei[0,1])**2,0,cutoff) +   \
                np.clip((x-nei[1,0])**2,0,cutoff) +  \
                np.clip((x-nei[1,2])**2,0,cutoff) +  \
                np.clip((x-nei[2,1])**2,0,cutoff))
    return (Phi+Psi)
其中函數內的變量Phi和Psi分別對應到Eqs. (6,8)。我們看一下energy函數的輸入，其中p不意外就是img[r,c]處帶有噪聲的原始像素值\(p\)，而nei則是圖4九宮格的陣列，但因為只用到4位鄰居，這就是為什麼nei只取了[0,1]、[1,0]、[1,2]和[2,1]四個位置，分別對應到九宮格橘色的方框。lam對到Eq. (8)中的\(\lambda\)而np.clip用來截掉超過cutoff的值，即\(\Psi_{\rm trunc}\)。原則上能量函數會回傳一個含有256個元素的能量陣列ene，其第0個位置表示像素值\(x=0\)所對應到的能量，第2個位置表示\(x=1\)所對應到的能量⋯⋯到第255個位置表示\(x=255\)所對應的能量，所以決定性算法Eq. (9)告訴我們最有可能的像素值\(x\)便是擁有能量最小的那個。假定能量陣列ene中在位置129處有能量最小，這告訴我們在[r,c]處除去噪聲後的像素值應該取129，這可以用img[r,c]=np.argmin(ene)來實現。

圖4. 圖像陣列操作示意圖，我們只取橘色的4個位置做鄰居\(\mathcal{N}\)

　這樣我們就完成決定性算法的介紹了，非常簡單吧，我們唯一要做的就是逐一把img中的每個位置都掃一遍，並用上面的方式來決定每個位置的像素值。不過在圖4中會發現如果[r,c]剛好像在[0,1]這種邊界的位置，我們取的九宮格會超過邊界，這時python就會鬼叫報錯，一個可行的方法是在img的邊界處更補上一個0，這樣陣列的大小就變成(R+2,C+2)，其中最外圍的像素值都是0，這個步驟叫做padding，如圖5所示，這與現在介紹卷積神經網路的padding是一模一樣的概念，另外padding可以用np.pad來處理。

圖5. 圖像陣列的padding示意圖

3.2 隨機性算法，像素值的採樣

原則上計算256個像素值每個出現的機率是來自Eq. (10)，而其中所需的輸入能量陣列nei與3.1節的結果相同，只是我們不直接取最小的那個能量，而是多做一步用nei所計算出的機率並使用Gibbs採樣來隨機選出可能的像素值。實現這個採樣的方法由下面的代碼給出 
def pix_sampling(ene):
    """
    Sampling possible pixel value by given energies for all 256 values, Eq. (10)
    ...
    """
    ener = np.asarray(ene-np.min(ene))
    numerator = np.exp(-ener)
    denominator = np.sum(numerator)
    prob = numerator/denominator # Eq. (10)
    return np.argmax(np.random.multinomial(1,prob))
其中numerator和denominator分別是計算Eq. (10)的分子和分母，相除後得prob是一個有256個值的機率陣列，第0個位置的值代表像素值0出現的機率，以此類推，與3.1節的能量陣列的解讀方式相同。最後我們將prob丟進numpy的多項式分佈multinomial中生成隨機數，原則上這個函數最後也會回傳一個共有256個元素的陣列，其中只有一個位置會出先1，其它位置都是0，這個出現1的位置就是我們的像素值，我們用np.argmax取出是1的位置。好了，就這樣我們也完成了隨機性算法的去噪方式，一樣我們需要掃過圖片的每個元素來完成一次迭代，腫膜樣，簡單吧！

3.3 Let's do it!

當然我不會詳細的介紹要怎麼架構修復的方式，但我把寫好的腳本放在github上，連同這篇專題的notebook載下來後利用 
import sys
sys.path.append('libs')
from denoise import denoise_det, denoise_prob
from utils import gaussian_noise,psnr
讀入我自己寫的去噪類，後綴det和prob分別表示用決定性和隨機性算法來去噪，當然還有我們生成高斯噪聲和計算PSNR的函數也不要忘記讀進來。從這裡開始我假設讀者都沒有也不需要像第1節一樣自己寫一個噪聲生成和計算PSNR的函數了，請直接用我寫好的腳本來處理，讀者只要去享受數學的威力即可，不需要現在架構腳本的泥淖裡 (除非你有興趣啦😂)！當然也不要忘記一併讀入cv2模組。

圖6. 左為園乾淨的圖像，右為受高斯噪聲污染後的圖像

　現在我們現在從隨附的資料夾images內輸入一張貓貓的照片kit.jpeg，並用下面的方式替它加入高斯噪聲 
kit = cv2.imread('images/kit.jpeg')
kit_noi = gaussian_noise(kit, sd = 15)
乾淨和被噪聲污染的貓貓照如圖6所示。接著我們要嘗試修復被污染的照片kit_noi，我們可以比較兩種算法的結果，如下 
kit_det = denoise_det(kit_noi)
kit_prob = deboise_prob(kit_noi)

# Iterate 3 times
its = 3
for i in range(its):
    kit_det.denoise()
    kit_prob.denoise()
在上面的代碼中當把受污染的圖片當作輸入放進denoise_det/prob類後就結束了，讀者原則上不需再輸入其它的參數，當然像Eq. (9)內的\(\lambda\)和\(\Phi_{\rm trunc}\)可以透過lam和cutoff來變更，它們預設分別是1和1000。當準備完成後，使用後綴denoise()來完成一次去噪的迭代過程，在這個的過程裡，每個位置的像素值都會依照Eq. (10)的方式來決定可能是不含噪聲的值\(x\)。通常一次迭代的結果不見得會太好，在上面我們迭代了3次，並將圖畫在圖7中，可以使用後綴status()來叫出這一次迭代的結果，例如使用plt.imshow(kit_prob.status()[:,:,::-1])來顯示去噪後的圖樣。

圖7. 去噪後的結果，左為使用決定性算法，右為使用隨機性算法

　沒有意外的，在圖7的兩張成果圖後噪聲消去的差不多了，PSNR也從原本的24 dB上下大舉躍升到31 dB左右，代表我們修復的結果很不錯，劣化降低非常多。大體上這便是實作去噪的主要內如了，在下一節裡我們來探討一下這些結果的意義。

四、結果探討

在上面的計算中，我們取了迭代3次，所以你484在想如果我把迭代多次一點，會不會像IA8的binary image那樣變得更好。Hmm，這是個好問題，但原則上如果優化不會再增加的話，迭代在多次也是浪費時間和計算資源而已，至於優化是否有變化當然可以用PSNR指標來看待，我將迭代10來次每次PSNR的值記錄如圖8。

圖8. 不同算法的PSNR vs. 迭代次數的變化

　從圖8很明顯的看出來，不管是哪種算法，通常在第4-5次迭代後PSNR便呈現持平的狀態，所以原則上去噪不需要迭代非常多次，通常5次或者多一點差不多就可以停止了，不過這當然是要依據照片被噪聲污染的程度，一個好的方法是追蹤每次迭代的PSNR，如果近3次迭代變化率小於某個閾值 (threshold value)，那麼我們就停止迭代，因為即使像素值還在變化，但這樣的變化並不會提升PSNR的指標，所以不需要再多花費資源在計算上。
　另外圖8也顯示出一個很顯著的結論，那就是隨機性算法並沒有顯著的優勢，比如PSNR勝過決定性算法之類的。如果替迭代一次計時的話，這張\(600\times 200\)的kit.jpeg利用決定性算法花費了5.3秒而隨機性算法需耗費9.4秒，可想而知隨機性算法非常的耗時，若圖片再增更大則需花費更多的資源，所以去噪並沒有特別的必要使用隨機性算法。
　而這是為什麼呢？原因出在能量陣列，如果讀者有興趣取出某個位置所計算的能量陣列ene來看的話，會發現這些元素少從值數百到數萬不等 (在我考慮四連通的算法且取\(\Phi_{\rm trunc}=\infty\)時能量的上限會是\(5\times255^2\)，你知道為什麼嗎?)，當我們利用Eq. (10)計算機率時別忘了它是放在指數裡，所以即使兩個像素值\(x=128\)和\(x=200\)所對應的能量值分別為\(E(x=128)=355\)和\(E(x=200)=367\)看起來很靠近，但機率通常是能量取指數，即\(p(x=128)\propto e^{E(x=128)}=e^{-367}\)，相較於\(p(x=200)\propto e^{-355}\)是非常的suppressed的 (差了\(\frac{p(x=128)}{p(x=200)}=e^{-367+355}=e^{-12}=e^{-|\Delta E|}\))，所以即使用採樣的方式，其結果與決定性算法直接選像素質\(x=128\)這個擁有最小能量的方式基本上無二致，這就是為什麼隨機性和決定性算法兩者去噪的結果是不分軒輊的，但有一點隨機性算法是完全被打趴的劣勢，那就是計算效能，如同前述的計時結果，這時候除非必要，我們應採取決定性算法即可。
　眼睛比較好的讀者仔細看我們修正的圖片，雖然PSNR顯著的提高表示劣化降低了，但仍會發現我們的圖看起來好像有些細節失去了 (和Ground truth比起來)，一個形容就像是畫的比較擬真的水彩畫，這是由於我們的能量函數中\(\Psi\)的中心像素值\(x_i\)和它鄰邊的差異\(x_j\)取的是二次方，統計學裡這叫\(L_2\)-norm loss function，這種解釋方式是用歐氏距離，真正的\(x_i\)應該是和它的4個鄰居的距離最靠近 (也就是擁有最少的方均根誤差)，或者說這5個像素值應分屬同一個群 (cluster)，且該群之間像素值的差異是遵守高斯分佈。講這麼多的統計學語言，\(L_2\)-norm就預設了這幾個相鄰像素值被Gaussian smearing了，所以這些像素值看起來都很接近，這也解釋了為什麼結果看起來像水彩畫一樣。但真正的自然圖片存在著非常多的細節，包含了邊角之類的，這些都代表了它們之間的差異不遵守高斯分佈，單純選用\(L_2\)-norm loss function反而會損失這些細節，但討論更細緻的loss function的選擇不是這篇專題所要做的，有興趣的讀者可以參考影像處理或是機器學習的專書。
　另外還有一點我忘了提，那就是在pix_sampling函數裡我一開始減去一個np.min(ene)這是因為機率是求在指數的能量，而當能量為數百的時候這個機率即使是所有像素值中最大的，它仍然超過機器本身的精度 (32-bit float)，所以在指數內扣掉等同是在Eq. (10)上分子分母同除一個相同的因子，這不影響結果，但會讓電腦可以計算，否則機率陣列會全部回傳0，根本無法採樣。

關於腳本的使用方式

在github上的主資料夾內會有一個IA10_image_denoise.ipynb的notebook檔以及兩個資料夾images和libs，主要的結構如下
  • 📄IA10_image_denoise.ipynb
  • 📂images
    • 📄kit.jpeg
    • 📄lena.jpeg
  • 📂libs
    • 📄denoise.py
    • 📄noisegen.py
    • 📄utils.py
範例以及函數的基本操作會透過notebook介紹，但這裡我大概簡述一下我提供的函數的功能及使用方式。

denoise.py

denoise.py內可操作的內容主要有兩個class與一函數，使用方式如下
  1. denoise_det(im, lam = 1, cutoff = 1000)
    使用決定性算法的去噪class，輸入im為圖像陣列，可以是彩色或是灰階，建議使用cv2.imread()的輸入。lam和cutoff分別對應到Eq. (8)的\(\lambda\)和\(\Psi_{\rm trunc}\)。
    操作方法：
    • denoise()：無輸入，完成一次整張圖片的去噪
    • status()：無輸入，回傳完成一次迭代後去噪的成果，可以藉由執行denoise()多次來提升去噪的結果
  2. denoise_prob(im, lam = 1, cutoff = 1000)
    同denoise_det()，但使用隨機性演算法
  3. denoise_multif(imls)
    使用Eq. (3)透過多張同場景的噪聲圖來去噪，輸入imls為含有多張噪聲圖的list，每個元素是cv2所讀入的影像陣列
在這個腳本裡還有其它函數主要是去噪class的subroutines，即是文章所提的生成256個像素值對應的能量陣列的函數energy()，以及如何從能量陣列取樣出最有可能的像素值pix_sampling()，它們通常不能單獨使用。

noisegen.py

noisegen.py內含有4種生成噪聲的方式，使用方式如下
  1. gaussian_noise(im, mean = 0, sd = 10)
    輸入圖像陣列im，mean和sd分別表示噪聲的均值和標準差，生成一張帶有高斯噪聲的im
  2. poisson_noise(im, a = 5)
    生成一張帶有Poisson噪聲的im
  3. saltpepper_noise(im, noi_f = 0.2, pep_f = 0.5)
    生成一張帶有椒鹽噪聲的im，noi_f控制im中有多少百分比被椒鹽噪聲污染，pep_f控制椒鹽噪聲中胡椒 (即全黑的椒點) 佔的比例
  4. uniform_noise(im, amp = 10)
    生成一張帶有均勻噪聲的im，噪聲的均值0，振幅為amp

utils.py

utils.py僅含一個計算PSNR的函數，使用方式如下
  1. psnr(im, tr_im)
    輸入im是帶有噪聲的圖片，tr_im是不帶噪聲的im，回傳PSNR值，單位是dB
好了，所有在github上關於這個專題所提供的腳本以及其函數的使用方式大概就是這樣，範例操作請看notebook檔吧，Enjoy!



Source on github: https://bit.ly/2GoJQvv
(請下載這個位置裡所有的檔案並放到同一個資料夾下)

留言

張貼留言

這個網誌中的熱門文章

物理雜記3 - Ising模型與模擬退火

圖片
Ising模型主要是用來描述物質的鐵磁性 (ferromagnetism)，是由德國物理學家Ernst Ising在1920左右提出了處理一維的解析解。該模型包含了每個晶格上原子的自旋 (spin) 參數，兩兩相鄰的原子彼此間可以是偏好處在共同的態上 (都是自旋向上)  或是處在相反的態上 (一上一下)，而外部的磁場會對該晶格上的自旋狀態施加能量，導致其偏向某一個方向。在這篇文章裡，我們不重複Ising當年複雜的數學解析解，而且該解也只能適用在1D的狀態，對於2D以上的系統目前並沒有解析解，所以我們將要另闢蹊徑，採用MCMC模擬的方式來決定一個系統終態每個晶格上自旋的排列情形。 一、Ising模型簡介 考慮一個由圖1所構成的晶格系統，晶格上的點都代表一顆原子，而黑白則表原子的自旋狀態，我們用黑表示向下而白表示向上，所以這個binary系統可以用更簡潔的\(\pm 1\)來表示每個原子可能的自旋 (即這個系統的樣本空間\(\Omega=\{+1,-1\}\))，另外為了方便，我們採用物理學的習慣將自旋表示為希臘字母\(\sigma_{n}\)，下標\(n\)可以想像成在 n 處原子的自旋態。原則上討論的都是一個由費米子 (fermions) 構成的系統，所以自旋應該用\(\pm \frac{1}{2}\)來表述，但在計算機上我們用\(\pm 1\)比較乾脆且簡潔。另外補充一下冷知識自然界已知的物質幾乎由費米子構成，僅有媒介作用力的粒子是玻色子 (bosons)，它們的差別在於費米子有半整數 (即不能約化成整數的分數) 自旋而玻色子則是整數自旋。 圖1. 晶格系統，黑白點分別代表該處原子的自旋狀態，黑表向下，白表向上。紅色的箭頭指出中心原子的4位鄰居\(m\) 　Ising模型假定了每個晶格上的原子的自旋都會受其鄰居 (neighbors \(\mathcal{N}\)) 的狀態影響，例如圖1中間白色的原子雖然他在這一刻是處在\(+1\)的狀態，但根據其週邊四個鄰居 (分別是3白1黑) 的況，它仍有一定的機會會在下一刻翻轉自旋的狀態。 　對一個給定位置 n 的原子，它和它的周圍鄰居構成的局部區域 (local area) 形成的一個系統的能量由Ising模型給出 \[ E_n = -\sum_{m\in \mathcal{N}} J
閱讀完整內容

PS7 - 小專題: 為何高斯分佈擁有最大熵?

圖片
高爾頓板，或稱做quincunx，是由Galton爵士 (1822-1911) 所發明用來驗證中央極限定理 (CLT) 的裝置。它也是第一位討論均值回歸 (regression to the mean) 現象的科學家，但同時亦是優生學 (eugenics) 這樣錯誤觀念的推廣者 在 上篇雜記 的結尾我們提到了對於一個趨於穩態的系統擁有最大的熵，而這也是統計力學的直接應用，用物理的角度來講便是當系統平衡時 Gibbs 自由能不再變化 (\(\Delta G =0\))，而在此時會對應到系統的熵有全域極大值 (global maximum)。故大自然很多現象都已經處於穩態的狀況下，我們對這個現象的母體做抽樣，所得到的觀測樣本應也有最大熵，則描述這些觀測樣本最好的分佈便是符合該條件下擁有最大熵的機率分佈。而指數族的分佈符合最大熵的特點，我們在這個小專題內便嘗試證明高斯與二項式分佈在連續和離散的隨機變量下的條件擁有最大熵的性質。 一、機率熵 在 上篇雜記 裡裡我們探討了夏儂熵，或者稱作資訊量\(I\)的意含，所以一個系統內 (我們這裡就不再用系綜這種彆扭的歷史名詞了) 每個成員出現的機率是\(p_i\)，則對該系統的熵可以用夏儂熵公式來表示 \[ I(p) = -\sum_i p_i \ln p_i. \tag{7.1} \] 所以對某一組 連續 機率分佈函數\(p(x)\)，其中\(x\)是隨機變量，則將(7.1)的求和符號改成積分後有 \[ I(p) = -\int p(x)\ln p(x) dx, \tag{7.2} \] 稱為機率分佈\(p(x)\)的 機率熵 ，而積分的上下限則看隨機變量的範圍，例如高斯分佈就是\(\pm \infty\)。如果今天我們有兩組機率分佈\(p(x)\)和\(q(x)\)，我們便可以定義它們之間的 相對熵 ，或稱作 KL散度 (Kullback-Leibler divergence)，定義為 \[ D_{KL}(q,p)\equiv\int q(x)\ln \left(\frac{q(x)}{p(x)}\right) dx = I(q,p)-I(q)\geq 0, \tag{7.3} \] 是 正定的 (positively defined)，其中 \[ I(q,p)=-\int q(x)\ln p(x) d
閱讀完整內容

IA9a - K-means分群與EM算法: 理論

圖片
市場內分門別類的產品 原則上 K-means分群 是 EM算法 (expectation-maximization algorithm) 的一個特例，相較於EM從似然函數 (likelihood) 的方式出發並導出不同資料點對不同群 (或稱聚類、簇，cluster) 的貢獻度，K-means單純用歐氏距離來分類，雖然這加快了K-means算法的計算效率，但顯然的這麼naïve的方法會產生明顯的缺陷。但仍讓我們會先從簡單的K-means開始，並探討K-means算法的結果在某些情況下和人類肉眼辨識的結論會產生明顯的歧義，接著再推進到通用的EM算法，除了 python 的demo外，我也會證明在數學上如何從EM得到K-means的極限 (僅限高斯混合的情況)，而在這個part內，讓我們先從數學理論的角度來面。 一、K-means分群 今天如果我有一大坨資料共120個點如圖1所示，從資料的散佈趨勢來看，很明顯的分成兩群，那我們是不是可以有一個演算法來自動化將資料點如同我們臆測的一樣分成兩袋？一個最簡單的發法便是使用K-means分群，而這個大寫K表示群的數量，以圖1而言，我們猜測K=2，所以是個2-means分群問題。 圖1. 資料散佈圖 (普魯士藍) 與初始化的紅綠兩點 　所以為了達成分群，我們需要先初始化兩個點，這兩個點代表兩群的中心。事實上，初始化的位置不必要一定剛好在兩群的中心，可以是任意挑選的地方，但為了計算方便，我們會將起始化的兩點開始在靠近資料點的平均值附近，如同圖1中紅綠兩點。我們將這兩個初始化的點叫作\(p_k\)其中\(k=1,2\)表示第\(k\)群。另外將資料點的集合用\(\mathbf{D}=\{d_1,d_2\cdots d_n \cdots d_{120}\}\)表示，所以我們計算每個點\(d_n=(x_{d_n},y_{d_n})\)和\(p_k=(x_{p_k},y_{p_k})\)的歐氏距離 \[ \ell_{n,k}=\sqrt{(x_{d_n}-x_{p_k})^2+(y_{d_n}-y_{p_k})^2}. \tag{1} \] 若該點和\(p_1\)的距離小於和\(p_2\)的距離 (\(\ell_{n,1}<\ell_{n,2}\))，則我們說\(d_n\)比較靠近\(p_1\)且應該
閱讀完整內容