IA7&PS9 - 貝氏推論下的圖像修復技術
 |
| 一張被亂塗鴉的臭貓貓照片 |
將刮花或是有噪點的圖像修復是一個有意思個課題，已經有許多的算法來幫助將受損的圖像復原，例如使用流體動力學的Navier-Stokes方程將像素色彩視為是一種流，這是python-OpenCV內建的修復方式 (cv2.inpaint)。而在在這篇文章中我們要另闢蹊徑，採用機率的觀點，從鄰邊像素的值來建立一個簡單圖像修復的方法，除了灰階的照片外，更將之應用於彩色照片上。對於那些不感興趣怎麼建立復原的方式，但感興趣直接使用修復的腳本，我建立的一個python的class叫imres (IMage REStoration的縮寫) 存放在imtools.py中，可以直接將受損的照片輸入並修復，操作方法可以直接跳第四節的部份，執行實作的腳本需要python-OpenCV套件cv2來讀取並做初步的圖像操作。本文也是IA和PS兩系列的交叉集，其中需要有統計力學配分函數 (partition function) 的概念，但我想就算不懂，還是不影響閱讀的 😂
一、都是像素惹的禍
在本文開始的cover image顯示了一張加菲貓在舔冰淇淋的圖片，但當中有部份被青色的塗鴉刮花了，我們的目的便是想修復這張受損的圖片，補齊缺失的像素值。我們都知道數位圖片是由像素構成的，每個像素給定一個顏色，例如cover image就由1000⨉773個像素構成 (通常是寬W⨉高H)。這個成像可以視為是當初入射的光子打在感光元件CCD上並被紀錄輸出成一張數位相片，假定每個像素是由一個帶有顏色的光子 (頻率) 所構成1，我們可以問一個問題：當給定某個位置的像素顏色時，其周圍的像素最有可能是什麼顏色？
 |
| 圖1. 觀測到像素\(z_i\)，對應到真正乾淨但未知圖片的\(x_i\)像素，而\(x_i\)與它鄰邊像素\(x_j\)是有關連的，這關連的連通關係構成了一張馬可夫毯 (Markov blanket) |
這個問題我想從看圖著手會比較容易點，在圖1中心的\(z_i\)是觀測到的像素，它可能是有雜訊或被刮花，總之是corrupted，不過我們相信他來自一個真正但未知的像素\(x_i\)，而\(x_i\)與其鄰邊 (周遭) 的顏色必然是有一定相關連的，比如說\(x_i\)是黃色，那麼他就有很高的機率是被黃色或是黃色相仿的顏色所構成的鄰居\(x_j\)所包圍。在圖1中我們用實現表示相關連，例如\(x_i\)與東西南北四個方向的像素質有關連，這變構成了一個四連通的關係，若把虛線也考慮進去，則\(x_i\)變與周圍八個方向都有關連，是一個八連通關係，而這個連通的關係用圖表示像個毯子一樣，所以我們叫它馬可夫毯 (Markov blanket)，至於為什麼叫馬可夫毯呢？這只能留待後面一點再解釋了。
故一開始的問題可以表示成：當觀測到\(z_i\)時，它所對應到某個\(x_i\)的機率可以表示成 \[ P(x_i|z_i),\tag{1} \] 上式表一條件機率，而最有可能為真的\(x_i\)的便是使Eq. (1)有最大值，或等價於求2 \[ {\rm argmax~} P(x_i|z_i). \tag{2} \] 在繼續討論以前，我們先將條件機率Eq. (1)用貝氏公式改寫成 \[ P(x_i|z_i) \propto P(z_i|x_i)P(x_i). \tag{3} \] 這樣改寫的好處是要直接從\(z_i\)逆相工程回推\(x_i\)可能有困難，因為他可能對應到很多種情況，而且每種情況的機率可能都相仿變得無法決定。
但改寫後變成求當給定某個\(x_i\)後，它演化成\(z_i\)的機率\(P(z_i|x_i)\)。如果是8-bit色彩深度則\(x_i\)有\(2^8=256\)種可能，每一種被選到的機率是\(P(x_i)\)，這便是\(x_i\)的先驗機率 (prior probability)。至於先驗機率\(P(x_i)\)怎麼決定，原則上引用John Skilling在1984年的結論3：在沒有任何先驗資訊底下 (即完全地無知ignorant)，只要擁有最大熵 (maximum entropy, MaxEnt) 的分佈，就是最合理的prior假設，在這個問題則表示我們要先產生出一張prior image，而這張圖像每個像素的顏色都是在0~255之間隨機取值，這樣的prior image滿足MaxEnt的條件，看起來大概就像是圖2一樣，這樣充滿雜訊的圖片就是我們潛在可能是真實的圖片中每個像素\(x_i\)構成的prior image，我們會用它們做之後的圖像復原。
 |
| 圖2. 滿足MaxEnt的prior image (500⨉500)，左右各表灰階和RGB彩色，色彩深度每個通道階是8-bit |
另外\(P(z_i|x_i)\)這個描述\(x_i\)演化成\(z_i\)的機率稱做似然函數 (likelihood function)，而Eq. (1)的\(P(x_i|z_i)\)則是\(P(z_i|x_i)P(x_i)\)的後驗機率 (posterior probability)，故能使後驗機率最大的\(x_i\)則為數學上最有可能是corrupted的\(z_i\)所對應到的真實\(x_i\)，這個問題現在是一個求解最大後驗 (maximum a posterior)，或是最大化似然函數 (maximum likelihood) 的過程。
1 實際上每個像素是由很多個光子擊中CCD的最小成像單元，這個單元會構成輸出相片的一個像素，不過為了方便起見我們視一顆光子會生成一個像素
2 argmax和argmin表示找到一個\(x_i\)的敘述 (argument) 使得該式有最大 (max) 或最小 (min)
3 J. Skilling and R. K. Bryan, MNRAS 211, 111 (1984)
2 argmax和argmin表示找到一個\(x_i\)的敘述 (argument) 使得該式有最大 (max) 或最小 (min)
3 J. Skilling and R. K. Bryan, MNRAS 211, 111 (1984)
二、圖像到底是什麼？淺談Markov隨機場
在談論修復圖片以前，我們先來看看圖片是什麼？對人類來說，圖片是一幅還有色彩與特徵的平面形狀，但這是像人類一樣以全局觀看的方式得到圖像的意義，而對於電腦，它本身並沒有所謂全局觀看的局面，就算他將整張圖片都掃過一遍，它也看不出色彩或是特徵。電腦真正看到的，是一張充滿格點的巨大網路，每個隔點都是由一個數字構成，或者更精確的來說是陣列。比如說我們將圖片用cv2讀入：
所以對電腦來說，它看到的圖片是一個陣列，每個陣列上的元素都在0~255中取值，而圖片刮花的地方因為毀損而不存在任何灰度值了，或者可以說那邊的顏色是沒有意義的，可以想成該處的值是NaN。而我們生成的第1張prior image代表我們刮花地方原本灰度值的初步猜測，而這些一開始的猜測會依據周遭沒有被刮花的像素的灰度值做調適，然後再修正其自身的灰度值，得到第一次迭代 (iteration) 後修復的照片，見圖4。如果一次迭代後對刮花部份的顏色修正不滿意，我們就將第一次迭代後的圖像當作新的prior image送進第二次迭代中，不斷重複這個過程直到我們對刮花部份修正的結果滿意而停止。當然也有可能當迭代夠多次後，往後的迭代在刮花部份的值都不再變化，我們就應該停止收手了，因為這張圖片上每個刮花處的像素灰度質已經達成穩定，或者用物理的術語來說該系統已經達成穩態 (stationary state)，不會再變化了，再多的計算不過是浪費時間而已。
言歸正傳，再未達成穩態前每次迭代後刮花部份的灰度值都是基於上一步在同個像素的灰度值與其相鄰觀測到的像素灰度值放入某個正規化函數後於0~255之間隨機取值 (見圖4)，而與上上步或更早之前迭代沒有直接的關連了 (decoupled)，所以這一連串的過程就是個馬可夫鏈 (Markov chain)，而圖片是二維的陣列，我們稱二維陣列是一個場 (field)，既然灰度值是隨機的，則稱做隨機場 (random field)，但這個隨機又符合馬可夫鏈的性質，且迭代多次會得到多個隨機場 (random fields)，整個修復的過程我們就稱做馬可夫隨機場 (Markov random fields, MRF)，見圖5，所以這也回答了圖1中為什麼我們叫連通關係是一個馬可夫毯的緣故。另外由於每次迭代所修復得到的圖片 (一個MRF) 都不是原來可以直觀的結果，而是原本被隱蔽起來的，所以這種又稱做隱馬可夫隨機場 (Hidden Markov random fields, HMRF)。
import cv2
img = cv2.imread("/pic/at/somewhere/img.jpeg")
type(img)
img.shape
 |
| 圖3. RGB影像陣列，數字代表該色彩通道內各個的8-bit灰度值 |
所以對電腦來說，它看到的圖片是一個陣列，每個陣列上的元素都在0~255中取值，而圖片刮花的地方因為毀損而不存在任何灰度值了，或者可以說那邊的顏色是沒有意義的，可以想成該處的值是NaN。而我們生成的第1張prior image代表我們刮花地方原本灰度值的初步猜測，而這些一開始的猜測會依據周遭沒有被刮花的像素的灰度值做調適，然後再修正其自身的灰度值，得到第一次迭代 (iteration) 後修復的照片，見圖4。如果一次迭代後對刮花部份的顏色修正不滿意，我們就將第一次迭代後的圖像當作新的prior image送進第二次迭代中，不斷重複這個過程直到我們對刮花部份修正的結果滿意而停止。當然也有可能當迭代夠多次後，往後的迭代在刮花部份的值都不再變化，我們就應該停止收手了，因為這張圖片上每個刮花處的像素灰度質已經達成穩定，或者用物理的術語來說該系統已經達成穩態 (stationary state)，不會再變化了，再多的計算不過是浪費時間而已。
 |
| 圖4. 透過prior image猜測刮花部份的顏色和觀測到未刮花部份的顏色來對刮花的部份做修復 |
言歸正傳，再未達成穩態前每次迭代後刮花部份的灰度值都是基於上一步在同個像素的灰度值與其相鄰觀測到的像素灰度值放入某個正規化函數後於0~255之間隨機取值 (見圖4)，而與上上步或更早之前迭代沒有直接的關連了 (decoupled)，所以這一連串的過程就是個馬可夫鏈 (Markov chain)，而圖片是二維的陣列，我們稱二維陣列是一個場 (field)，既然灰度值是隨機的，則稱做隨機場 (random field)，但這個隨機又符合馬可夫鏈的性質，且迭代多次會得到多個隨機場 (random fields)，整個修復的過程我們就稱做馬可夫隨機場 (Markov random fields, MRF)，見圖5，所以這也回答了圖1中為什麼我們叫連通關係是一個馬可夫毯的緣故。另外由於每次迭代所修復得到的圖片 (一個MRF) 都不是原來可以直觀的結果，而是原本被隱蔽起來的，所以這種又稱做隱馬可夫隨機場 (Hidden Markov random fields, HMRF)。
 |
| 圖5. 圖像修復的過程，每一次迭代後的像素灰度值都僅與上一步該處的值有關，對單個像素而言整個過程構成一個馬可夫鏈，而整張二維的圖像陣列 (場) 的修復過程則是馬可夫隨機場 |
三、推論色彩
從圖5的示意圖來看，我們先隨機給定某個刮花處像素的灰度值 (淡紫色)，這就是我們的prior，而其值在0~255隨機取符合MaxEnt的條件，接著依照週邊觀測的情況 (observed，prior左邊半透明淺綠色)，prior對其自身做修正得到了橘色，這個橘色就是prior和observed透過最大化Eq. (3)右式得到的posterior，這就是一次迭代後修復的結果，但如果我們對這個結果仍不甚滿意，就接著以橘色當prior配上observed再進行第二次迭代得到了posterior淺橘色，不斷重複這個過程直到滿意或顏色不再變化為止！
 |
| 圖6. cover image中修復過程眼睛部份的放大圖，左邊為刮花部份的prior image，右邊為多次迭代後的結果 |
舉cover image的貓貓為例，左眼部份已經被刮花了，為了修復貓貓的眼睛，在被刮花的部份我們填上先驗的像素灰度值，看起來如圖6左圖，那些花花綠綠七彩顏色的便是對刮花的像素初步猜測滿足MaxEnt條件所填上的prior image，因為只有括花的部份需要修復，我們只對刮花處套用prior image，其他沒有被刮花的部份則維持不動，而右邊為一次迭代後依據鄰近像素修正後的結果，整張圖看起來像圖7，右邊修復的部份的灰度值都是對應到左圖刮花部份的prior和\(P(z_i|x_i)\)的posterior。
 |
| 圖7. 同圖6，prior image和一次迭代後的修復圖 |
3.1 求後驗\(P(x|z)\)
所以由於\(z_i\)是給定的，那麼求最有可能的\(x_i\)就像是統計力學裡給定巨觀條件後，微觀狀態最有可能的分佈，這原則上是一個Gibbs distribution:
\[
P(x|z) =\frac{e^{-E(x,z)}}{\mathcal{Z(z)}}\propto e^{-E(x,z)},\tag{5}
\]
其中\(E(x,z)\)表能量函數而\(\mathcal{Z(z)}\)是配分函數。所以Eq. (2)最大化\(P(x|z) \)變對應到最小化能量函數，或是
\[
{\rm argmin~} E(x,z). \tag{6}
\]
在圖像重建的方式中，能量函數最常取的形式是
\[
E(x,z) = \sum_{i \in\mathcal{V}}\Phi_i(x_i,z_i)+\sum_{(i,j)\in \mathcal{E}}\Psi_{ij}(x_i,x_j,w).
\]
在上式中奇怪的表示\(i \in\mathcal{V}\)代表對這張圖的每個像素\(i\)逐元素計算求和，而\((i,j)\in \mathcal{E}\)表示計算每個\(i\)和週邊像素的關係，即圖1的馬可夫毯。刮花修復的優點是我們不用管\(\Phi\)這個修正項，其恆為0 (可以見A. Blake et al.的著作MRF for vision and image processing)，所以真正重要的是
\[
\Psi_{ij}(x_i,x_j)=\min[(x_i-x_j)^2,\psi_{\rm cut}].\tag{7}
\]
可以見到重建刮花的部份與刮花部份的觀測值\(z_i\)是無關的，這也是合理，因為前面說過刮花的部份本來就不存在可見的像素，因為已經毀損了，其值可用NaN代替，直到我們將它用prior image初步亂猜的像素灰度值取代為止。將Eq. (6)套入Eq. (7)後我們得到修復後的像素\(x_i\)必須滿足條件
\[
x_i = {\rm argmin}\sum_{(i,j)\in \mathcal{E}} \min[(x_i-x_j)^2,\psi_{\rm cut}], \tag{8}
\]
其中\(j\)是與\(i\)連通的鄰近像素灰度值，如果是四連通則任意\(i\)會對應到4個\(j\)，八連通則8個，多連通則更複雜，稍候我們要建造的是20連通的馬可夫毯。在Eq. (8)中的\(\psi_{\rm cut}\)表示截斷值，它決定了\(x_i\)和鄰近\(x_j\)灰度值的差值平方是否可以任意取或者必須小於某個截斷值。
到這裡我們就初步完成理論的描述了，接下來我們就要用這些方式去建造一個能夠做圖像復原的python腳本啦！
首先需要有受損的圖片，可以在github的連結內將libs和images這兩個資料夾連同本篇對應的ipynb下載下來，首先我們需要引入imres和cv2模組，cv2幫助我們較好的操作與讀取圖像。將下載下來的3個東西放在同個資料夾內並在ipynb內執行
所以接著我們讀入影像檔，以images資料夾內的img1_xxx為例，示範如下
原則上輸入mask是必要的，否則無法定位刮花的位置，兩張圖的大小一定要相等，才能正確的對應位置，否則會報錯。這與cv2.inpaint修復相片的輸入一致。另外mask用cv2.imread讀入時不一定要是bmp，原則上他可以是任何cv2接受的圖檔，白色的部位指出刮花的位置。輸入後轉成binary image，刮花的地方binary value等於1。
接著我們要將受損的照片和遮罩引入imres類做初始化的動作，像這樣：
通常2~3次迭代後圖片就可以復原的不錯，當然如果差強人意，則可以迭代更多次，通常我們會一開始就選擇迭個5次之類的，可以用如下方式實現：
可以點圖放大後發現復原的蠻不錯的，另外在右眼的部份因為毀損的比較嚴重，原則上重建無法很靠近原圖，這是情有可原的，當然還有臭貓2幾個字還隱約可見😂。圖10和圖11是另外兩張迭代5次的復原圖，可以比較看看。
另外可以發現在爛水果左下番茄的部份，有蠻明顯接縫的痕跡，這是因為我們將RGB三通到都採取相同的截斷值\(\psi_{\rm cut}\)，這也許不是最好的方法，而是不同通道應該選擇不同的值，這樣接縫的痕跡也許會比較顯著的消去，這個部份就是要修改imres讓它不要這麼自動的一次處理3個通道，把每個通道拆開來分別取不同的\(\psi_{\rm cut}\)，不過我實在太懶了沒什麼興趣這樣做，雖然一開始我是分開不同通道處理來寫imres類沒錯，但為了懶人我只好將它自動化壓在一起了！
如果還有時間寫下一個有趣的計算機視覺，那我我覺得應該來寫寫超解析相片 (super-resolution image)，這就像看CSI影集當中將低解析度的相片重製成高解析度相片的技術，這應該也是非常有趣的，也可以用貝氏推論，或是現在很火的deep learning的方法來做！
Source on github: IA7_img_restoration.ipynb
到這裡我們就初步完成理論的描述了，接下來我們就要用這些方式去建造一個能夠做圖像復原的python腳本啦！
四、實作修復圖像
怎麼建立修復的python類imres這邊我就不解釋了，原理是基於第三節的方法，有興趣可以直接看imtools.py這個腳本，現在我們只要知道怎麼操作它就行了！首先需要有受損的圖片，可以在github的連結內將libs和images這兩個資料夾連同本篇對應的ipynb下載下來，首先我們需要引入imres和cv2模組，cv2幫助我們較好的操作與讀取圖像。將下載下來的3個東西放在同個資料夾內並在ipynb內執行
import cv2
import sys
sys.path.append("./libs")
from imtools import imres
所以接著我們讀入影像檔，以images資料夾內的img1_xxx為例，示範如下
img1 = cv2.imread("images/img1_clean.jpeg")
img1_s = cv2.imread("images/img1_scratch.jpeg")
img1_m = cv2.imread("images/img1_mask.bmp", 0) # 0 changes it to greyscale image
 |
| 圖8. 左邊是刮花的貓貓照，右邊是定位刮花的位置 |
原則上輸入mask是必要的，否則無法定位刮花的位置，兩張圖的大小一定要相等，才能正確的對應位置，否則會報錯。這與cv2.inpaint修復相片的輸入一致。另外mask用cv2.imread讀入時不一定要是bmp，原則上他可以是任何cv2接受的圖檔，白色的部位指出刮花的位置。輸入後轉成binary image，刮花的地方binary value等於1。
接著我們要將受損的照片和遮罩引入imres類做初始化的動作，像這樣：
img1res = imres(img1, mask = img1_m, cutoff = 1500)
plt.imshow(img1res.status()[:,:,::-1])
通常2~3次迭代後圖片就可以復原的不錯，當然如果差強人意，則可以迭代更多次，通常我們會一開始就選擇迭個5次之類的，可以用如下方式實現：
# Iterate 5 times
for i in range(5):
img1res.restore()
plt.subplot(121)
plt.imshow(img1res[:,:,::-1])
plt.title("Restored image (5 iterations)")
plt.axis("off")
plt.subplot(122)
plt.imshow(img1res.status()[:,:,::-1])
plt.title("Original image")
plt.axis("off")
 |
| 圖9. 迭代5次後的復原圖與原圖比較 |
可以點圖放大後發現復原的蠻不錯的，另外在右眼的部份因為毀損的比較嚴重，原則上重建無法很靠近原圖，這是情有可原的，當然還有臭貓2幾個字還隱約可見😂。圖10和圖11是另外兩張迭代5次的復原圖，可以比較看看。
 |
| 圖10. 另一張臭貓貓復原圖 |
 |
| 圖11. 爛水果復原圖 |
五、小結
在這篇文章裡，我簡介了如何透過機率和配分函數的概念來修復受損的圖像，原則上只要受損的位置不算太大，這個方法都可以不錯的復原，但如果受損的位置太大，或原始的照片本身就有非常多的細節的話，那麼復原的結果就會差強人意，如果不是作為鑑識證據 (forensics evidence) 的話，那麼可以考慮使用紋理合成 (texture synthesis) 的方法，但我不會稱這個方法叫作修復，這比較像是複製貼上，從旁邊未受損的像素當中選一塊區域最符合邊界接縫處然後複製，再拼貼到受損的部位，這比較像是讓電腦天馬行空任意剪裁貼圖塞進去，不過這個方法基本上很看不出痕跡，算是有意思的一種計算機藝術創作。另外可以發現在爛水果左下番茄的部份，有蠻明顯接縫的痕跡，這是因為我們將RGB三通到都採取相同的截斷值\(\psi_{\rm cut}\)，這也許不是最好的方法，而是不同通道應該選擇不同的值，這樣接縫的痕跡也許會比較顯著的消去，這個部份就是要修改imres讓它不要這麼自動的一次處理3個通道，把每個通道拆開來分別取不同的\(\psi_{\rm cut}\)，不過我實在太懶了沒什麼興趣這樣做，雖然一開始我是分開不同通道處理來寫imres類沒錯，但為了懶人我只好將它自動化壓在一起了！
如果還有時間寫下一個有趣的計算機視覺，那我我覺得應該來寫寫超解析相片 (super-resolution image)，這就像看CSI影集當中將低解析度的相片重製成高解析度相片的技術，這應該也是非常有趣的，也可以用貝氏推論，或是現在很火的deep learning的方法來做！
Source on github: IA7_img_restoration.ipynb
留言
張貼留言