SIC ITUR AD ASTRA

我的網誌移到 github 上囉! 以後的新文章都會發佈在 https://yenhsunlin.github.io/

一張帶有高斯噪聲的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_

為了IA9內容的完整性，我還是趕快把如何從高斯混合EM過渡到K-means的證明補齊吧，免得時間一久不但忘了也懶了😂 　讓我們從置信度開始，它是這樣表示的 \[ \gamma(z_{nk})=\frac{\pi_k \mathcal{N}(x_n|\mu_k,\sigma_k^2)}{\sum_{k^\prime}\pi_{k^\prime}\mathcal{N}(x_n|\mu_{k^\prime},\sigma_{k^\prime}^2)}\tag{1} \] 用來闡明該點\(x_n\)是屬於第\(k\)群的可能性，注意分母的\(k^\prime\)是個dummy variable，並且在上式中的參數捨去了粗體寫法，讀者必須知道它們都可以直接推廣到向量的形式。所以我們便可以計算出第\(k\)群的中心位置\(\mu_k\): \[ \mu_k=\frac{1}{N_k}\sum_{n} \gamma(z_{nk})x_n \tag{2} \] 和 \[ N_k=\sum_{n}\gamma(z_{nk}).\tag{3} \] 由於高斯分佈完整的形式為 \[ \mathcal{N}(x_n|\mu_k,\sigma_k^2)=\frac{1}{\sqrt{2\pi}\sigma_k}\exp \left[-\frac{(x_n-\mu_k)^2}{2\sigma_k^2}\right]\tag{4} \] 並將它代入Eq. (1)稍做整理後得到 \[ \gamma(z_{nk})=\frac{\beta_k \exp(-\beta_k^2 \ell_{nk}^2)}{\sum_{k^\prime} \beta_{k^\prime} \exp(-\beta_{k^\prime}^2 \ell_{nk^\prime}^2)}, \] 其中標準差的倒數\(\beta_k = 1/\sigma_k\)稱做 肛 …啊不好意思，是 剛度 (stiffness) 而\(\ell_{nk}\)即歐氏距離\(\sqrt{(x_n-\mu_k)^2}\)，又由於\(\pi_k\)是個常量 (\(z_k\)出現的先驗機率) 和我們在意的變數都沒有直接關聯，所以我將他等式中移除並不影響。 　現在我們就證明當剛度趨近無限大時，上述的高斯混合EM會變成K-mean的形式，故有

這一篇文章延續了上一篇IA9a的內容，但側重在 python 的實作，所以所有的節數、圖案以及公式編號都會依照IA9a的編號接序下來，所以重複的內容就直接引用，不在這裡重新贅述。 四、實作EM 4.1 高斯混合分佈 讓我們從1.1節的協方差矩陣開始，生成的圖形就如圖4 (腳本內加了seed random方便我們每次追蹤)，另外也看到了執行5次K-means算法後 (再多次也不會變了) 得到的分群結果如圖5，這顯然與我們的預期有極大的落差！ 　但我們將按照3.3結給出EM算法的偽代碼，假設資料群是高斯混合分佈的形式，重新執行一遍分群，我們將看到EM算法能很好的分出我們預期的兩群的結果。我們預計會有兩群，所以\(k=1,2\)，且第一步需要初始化參數組\(\boldsymbol{\theta}=\{\mu_k,\Sigma_k,\pi_k\}\)。初始化的過程是給定一個起始點，原則上隨機選擇的方式是合理的，不過為了計算上方便，可以讓初始的協方差矩陣都為單位矩陣，\(\pi_k\)則為每個群均分 (這裡就是每個都是0.5)，如下代碼 # Initializing parameters sigma1,sigma2 = np.eye(2),np.eye(2) pi1 = 0.5 pi2 = 1-pi1 N = len(cluster) 完成初始的參數選擇後，我們用偽代碼的第二步，也就是E-step計算置信度\(\gamma (z_{nk})\)。這裡我將Eq. (10)列出來方便解說 \[ \gamma (z_{nk})=\frac{\pi_k \mathcal{N}(\mathbf{x}|\mu_k,\Sigma_k)}{\sum_{k=1}^2 \pi_k \mathcal{N}(\mathbf{x}|\mu_k,\Sigma_k)}\tag{10}. \] 因為我們是假設兩群都是高斯分佈，所以對第一個群而言的分子就是 # Starting E-step: Calculating responsibilities # Numerator for k=1 k1 = pi1*multivariate_normal.pdf(cluster,p1,sigma1) 直接透過 python 腳本我們很容易理解了分子的意義便是計算每個資料點分別對應到\(\{

市場內分門別類的產品 原則上 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\)且應該

在IA4的文章裡我們提到了傳輸一張binary image後，有可能會受到干擾導致圖片被污染，裡面討論了復原的方法，像是一次傳輸多張照片，最後用投票的方式來決定該像素的顏色應該取\(+1\)還是\(-1\)，因為如果噪聲是隨機的，那麼不太可能同個位置超過一半的像素都是錯的。但這個方法有個缺點，就是必須一次傳輸多張相同的圖片才能達到良好的修復效果，可是這非常地消耗數據傳輸量以及講究頻寬的大小，這對應到成本的增加，尤其是當我們要傳輸高傳真的訊號時。那麼是不是有什麼方法，可以單從一張收到的影像中去回推原始未受到污染的圖片，我想要用兩篇文章來探討，我採用由簡入繁的方式，第一篇這裡我們先處理binary image的情況，在下篇我們再進一步處理彩色的照片。 一、從圖模型檢視雜訊圖像 在IA7裡我們談到了修復刮花圖像的技術，出於圖像毀損的狀態，所以我們當初設定是刮花處與真實圖像之間並無關聯，但對於雜訊圖像則情況不同，因為這些在圖像上的雜訊，可能是原來的像素加上一個擾動值，假設真實的像素值是\(x\)，那麼加上噪聲後新的像素\(y\)就可以表示成 \[ y = x + n\tag{1} \] 其中\(n\)表示噪聲。這種加上擾動值\(n\)最常假設的生成方式便是均值為0的高斯噪聲。 　從Eq. (1)中我們假設了一個簡單的關係式。由於一整張圖片是由非常多的像素所構成一個W×H的網格，或稱陣列，其中W和H分別表示圖像寬和高所對應的像素數量，我們可以用 圖模型 (graphical model, GM) 來對應真實無雜訊的圖像\(\mathbf{x}\)和觀測到存在噪聲的圖像\(\mathbf{y}\)之間的關係，見圖1，粗體字表示圖像是一個陣列而非純量。 圖1. 利用圖模型來推論真實圖像和雜訊圖像之間的關係 　我們認為真實圖像\(\mathbf{x}\)的每一個像素 (圖1中橘色實心的節點) 和周遭的像素都有關連，所以我們用紅色的實現將之連起，例如一張照片中一張人臉上任一處肉色的圖像一定和周圍皮膚的顏色有關聯，而不會和遙遠處天空的藍色有關聯。但我們觀測到的圖像\(\mathbf{y}\)已經被某些未知雜訊的生成機制所破壞，所以原則上這些像素點之間並不存在關聯，就如同圖1中藍色的空心節點，它們構成我們所見的圖案，彼此間沒有線將它們串連起來。而這個

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