SIC ITUR AD ASTRA

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

在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...

一張被亂塗鴉的臭貓貓照片 將刮花或是有噪點的圖像修復是一個有意思個課題，已經有許多的算法來幫助將受損的圖像復原，例如使用流體動力學的Navier-Stokes方程將像素色彩視為是一種流，這是 python-OpenCV 內建的修復方式 ( cv2.inpaint )。而在在這篇文章中我們要另闢蹊徑，採用機率的觀點，從鄰邊像素的值來建立一個簡單圖像修復的方法，除了灰階的照片外，更將之應用於彩色照片上。對於那些不感興趣怎麼建立復原的方式，但感興趣直接使用修復的腳本，我建立的一個 python 的 class 叫 imres  ( IM age RES toration的縮寫) 存放在 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\)是黃色，那麼他就...

高爾頓板，或稱做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...

PASIEKA / SCIENCE PHOTO LIBRARY Data analysis is simply a dialogue with the data —S. F. Gull 一、機率到底告訴我們了什麼? 假定我們有兩枚六面骰A和B每面為標記為1~6的整數，我們如何確定兩枚骰子都是公正的? 即每個面骰出的機率皆為1/6。古典機率的概念告訴我們若擲兩枚骰子多次，並記錄每面出現的次數，如果一枚骰子是公正的，那麼理論上每面出現的次數相同或接近相同。 圖1.  骰A (粉藍) 與骰B (粉紅) 各擲50次每面出現的次數 　現在我擲50次並記錄兩枚骰子每個面出現的次數，並用將其畫成直方圖如圖1所示，從記錄的結果來看，我可以斷定哪個骰子是公正的嗎? 看起來好像不行，那如果我改成擲100次、1000次就可以分別公不公正了嗎? 換個方向來看，也許我們不在乎骰子是否公正，我們在乎的是用我們觀測的資料集\(D\)，它是否可以告訴我們骰子每面出現的機率分佈是什麼? 都是同一組資料，一個是基於古典推論 (classical inference)，另一個則是貝氏推論 (Bayesian inference)。對於古典推論，即頻率論而言，我們有以下原則 (tenets)： 機率對應到的是事件出現的頻率，它們是真實世界的 客觀呈現 (objective properties) 參數是固定的 (未測量前是未知的)，比如六面骰每面擲出的機率可以是這裡提到的參數。這裡的參數對應的是柏拉圖世界的量，如果是一枚六面公正骰，那麼每面擲出的機率一定是1/6，不可能有波動 (fluctuation) 當我們重複步驟愈多次 (如擲骰的次數)，那麼每面出現的頻率應該是一個可以被定義的量 (即收斂，如果該量會產生波動則不收斂) 但對於貝氏推論而言，它則是基於如下的態度 (stance)： 機率是對 主觀信念 (subjective beliefs, prior ) 的描述，而不受限於資料出現的頻率長遠來看是否是一個定值 我們對於模型參數 (例如每個面出現的機率) 的推論給出的是它位於“某個值”的機率分佈，這個分佈代表的是我們對這個參數有多清楚 (不清楚的話該參數就會是一個很大的範圍) 關於古典推論的原則是沒有問題的，非黑即白，但對於貝氏...