發表文章

IA7&PS9 - 貝氏推論下的圖像修復技術

圖片
一張被亂塗鴉的臭貓貓照片 將刮花或是有噪點的圖像修復是一個有意思個課題,已經有許多的算法來幫助將受損的圖像復原,例如使用流體動力學的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\)是黃色,那麼他就...

IA6 - 小專題: 進擊的黴菌

圖片
初始3個菌點經過80單位時間後形成的菌落 利用電腦模擬真實世界的情況,不僅僅是電腦科學家,還包含了物理學家、生物學家、化學家和多種領域專家有志一同的想法。自然界千變萬化,圖騰多樣,但或許它們都可以從幾個簡單的規則當中得到答案, 這便是所謂的湧現現象 (emergence phenomenon)。 我們來看個簡單的例子,黴菌菌落的增值過程。 一、菌落的消長 用一個\(N\times N\)的區域代表一整片的培養區域如圖1,這個區域裡有比例\(q\)的部份是培養基質可以生長出黴菌 (淺綠),剩下的區域則不能長出黴菌 (白色),假設長出和不能長出的區域是隨機分佈的,那麼我隨意在這片區域灑上幾個菌點 (墨綠色),是問它們發展成菌落的過程為何? 圖1. 培養區域,淺綠是含有生長基質的區域,白色不能生長,墨綠色表示初始的5個隨機菌點  這個問題和Prof. Downey在他的著作Think Complexity的第七章在描述滲透的情況基本上一模一樣,除了關鍵的黴菌生長需要引入額外的判定條件外,所以我會依照Prof. Downey的方式來做,最後在引入這個判定生長的條件。欲執行這個雜記的notebook,我需要下載Prof. Downey所寫的一個 python 腳本 Cell2D.py ,可以從雜記的notebook裡直接點連結。這個腳本主要是需要其中的 Cell2DViewer 類來將黴菌菌落的消長輸出成動畫,如果只需要用 imshow 看特定時刻的結果,引入 Cell2D.py 不是必須的,但如果想執行本篇雜記的notebook第一個擴散的例子,則需要有其中的 Cell2D 類。 1.1 預備知識 如果說以圖1其中一個淺綠的基質為中心 (Eq. (1)的\(?\)處),那麼它週邊\(3\times 3\)的區域內4個東南西北 (不含對角線方向:東南、東北……等等) 的方向可能含有黴菌或不含有黴菌,如下 \[ \begin{array}{ccc} \times & {\rm 黴} & \times\\ {\rm 無} & ? & {\rm 無}\\ \times & {\rm 無} & \times \end{array}\tag{1} \] 所以我們想知道如果它周圍有一個以上的菌點,...

IA5 - 元胞自動機

圖片
生命遊戲 (Game of Life, GoL),取自維基百科 元胞自動機 (cellular automaton,複數automata,縮寫CA) 由John von Neumann在上個世紀50年代所提出,希望能用來模擬生物自我複製的情況,然後到了70年代由John Conway設計了生命遊戲 (Game of Life, GoL) 並刊上了科學美國人後才逐漸吸引其他研究者的目光。而在80年代由天才科學家,也是知名數學套件Mathematica®的開發者Stephen Wolfram利用統計力學的方式做了一系列的研究後,我們才對這類型的計算理論與結果有了系統性的理解。 一、元胞自動機的演化規則 CA通常由包含了\(N\)個元素的一維陣列所組成,每個元素不是\(1\)便是\(0\)。用比較圖像化的形式,將每個元素都想像成一個細胞,每個細胞不是黑\(1\)就是白\(0\)的兩種狀態,而整個陣列便代表在時間\(t\)時由\(N\)個細胞所構成的假想生物體的狀態,當然這個假想生物體就是由一系列\({\tt 101000100}\cdots\)來描述。  現在這個生物體的狀態由時間\(t\)演化到\(t+1\),那麼在\(t+1\)時這個生物體的每個細胞是處在\(1\)或\(0\)則由上一步的它自己與它左右兩個鄰居的狀態來決定,如下: \[ \begin{array}{cccccc} t & \cdots & {\tt 1} & {\tt 0} & {\tt 1} & \cdots\\ & & \searrow & \downarrow & \swarrow\\ t+1 & & & ? \end{array} \] 可以看見在\(t\)時三個相鄰的細胞狀態分別是\({\tt 101}\),則處在\(0\)的那個細胞它在\(t+1\)時則由它自己原來的狀態\(0\),以及鄰近那兩個處在\(1\)的細胞來決定。  所以說任三個相鄰細胞的狀態,可以決定中間那個細胞在下一個時刻時的狀態,由於每個細胞可以處在\(0\)和\(1\)兩種情況,三個細胞變有\(2^3=8\)種可能的組合: \[ \begin{array}{ccccccccc} t & ...

IA4 - 編碼、解碼與糾錯: 以二進制訊號為例

圖片
位於美國新墨西哥州的甚大天線陣列 (Very Large Array, VLA) 資訊 (information) 的傳遞在人類文明演進的歷史中扮演著極重要的角色,不僅僅是知識、技能,更是讓文化以及思想的瑰寶得以傳承,縱然生命會消逝,但前人所建立起的經驗卻能夠代代相傳給往後的人們。這些知識、技能和思想都可以歸納成資訊的一種,而在資訊遞交給後代的過程就稱做傳輸 (transmission)。傳輸不見得是保真的,亦即後代的人們所讀到前人的著作可能有多種版本,例如莎士比亞 (W. Shakespear) 的著作全集就有多種版本在流通,又或是1631年的邪惡聖經 (Wicked Bible) 事件,印刷廠將十誡中的 Thou shalt not commit adultery (不可姦淫) 誤植成 Thou shalt commit adultery (應當姦淫)。這些都是資訊在傳遞的過程當中發生失真,大部分的失真皆無大礙,並不影響資料原來意義的解讀,但少部分的失真卻是重大的,可能導致文意生變,甚至完全與當時書寫者的核心概念大相徑庭。所以既然失誤不可避免,有沒有什麼方法能讓我們盡可能的重建出原始資料的樣貌呢? 一、訊號的傳遞 讓我們來看看資料是怎麼傳輸的,以通電話為例 (圖1),當你朋友撥手機給你時,他的聲音會被手機錄下並編碼,接著通過基地台、電信機房,再到離你最近的基地台將訊號廣播出去,你的手機收到訊號後將之解碼重新轉成類比音訊,再透過喇叭將你朋友的聲音播放給你聽,同時地,若換你說話,則這個過程將逆向發生! 圖1. 聲音訊號如何透過手機傳到朋友的手機或話筒上,補充一下當代手機的編解碼的方式是利用類似貝氏機率的方式去重建可能的原始訊號,不是這篇專題內提到的\(R_r\)編碼  但注意訊號在傳遞的過程可能失真,就像你聽到手機播放你朋友的聲音時,總會覺得不像是親耳聆聽他的聲音一樣清晰,像隔著一層簾幕一樣,或甚至音調也有一些變化,有時背景還會傳來微若的雜訊滋滋聲。這便是圖1中間紫色虛線框裡描述的,當訊號透過層層基地台、電信機房時,原始的訊號發生了失真,可能被干擾了,或位元流失、誤植⋯等等,如果你的手機收到這樣破損的訊號,它甚至可能沒辦法重新解讀成類比音訊撥放給你聽。這時候,在一開始你朋友的手機做了一件事,那就是不僅僅是將你朋友的聲音送出去...

IA3 - 再看\(n\)維球上的點,淺談均勻性與不對稱性

圖片
哈伯深空圖像 (NASA/ESA/S. Beckwith(STScI) and The HUDF Team) 根據宇宙學原理 (cosmological principle),一般相信宇宙在大尺度底下 (>Mpc) 是均勻 (homogeneous) 且各向同性 (isotropic) 的,亦即宇宙不存在中心,就像在球的表面上行走一樣,我們不會找到哪個地方在球面上有特別突出的重要性一樣,另外在宇宙各處的物質密度也應該是均等的。這樣均勻、公正且不偏頗的特性,同樣的也在統計取樣中,特別是沒有任何先驗資訊底下扮演了重要的角色。這裡我們回顧上個專題I&A2中生成任意維度單位球面上點的例子,我們以為生成的點是均勻地散布在球面上,但實際上這是一個錯覺,為什麼我們生成\((x,y)\)的座標明明就是均勻的分佈,但最後會導致不均勻,或不對稱性的發生呢? 一、均勻性 讓我們回想一下上一個專題談到的均勻性,我們用二維單位圓上的點的例子來舉例比較容易想像。見下圖1,為了產生在單位圓上的點,我們可以隨機生成\((x,y)\)分別在\((-1,1)\)之間的點並正歸化它們,以形成單位圓上的點,而向量\((x,y)\)則用於決定在極座標上的角度\(\phi\)。可是從圖1會很明顯的發現紫色的點它們都是對應到相同的\(\phi\)角,而綠色的則是角度為\(\pi/2\),但由於\(\phi\)的位置的點可以超出單位圓,所以產生在\(\phi\)這個角度的單位圓上的點,實際上是比角\(\pi/2\)上的綠點來得多。亦即I&A2的算法2其實生成的並不是均勻分布的點,比如說在\(\phi=(\frac{\pi}{4},\frac{3\pi}{4},\frac{5\pi}{4},\frac{7\pi}{4})\)處的點 (二維的例子剛好在這裡的點是最多的) 會多於\(\phi=(0,\frac{\pi}{2},\pi,\frac{3\pi}{2})\)的點 (這裡的點是最少的)。 圖1. 藉由平均分佈生成單位原上的點,以及各個角度\(\phi\)所對應點出現的機率  我們可以重新run一次I&A2的算法2,取\({\tt dim}=2\),然後利用 \[ \phi = \tan^{-1}\left(\frac{x}{y}\right...

IA2 - \(n\)維球上的點

圖片
馬可夫過程 (Markov process) 生成三維球面上的點,此時點的位置不是均勻分佈,而是猶如醉漢走路 (drunkard's walk) 的情況。馬可夫過程在針對自然界運作或是社會結構演變的模擬中扮演著核心的角色! 一、二維圓盤的case 這次的專題非常簡單,我們先從二維的情況開始,二維球面上的點即為是分散在圓周上的點。我們知道如何產生散布在\((x,y)\)上的點 (圖1),但如何將之收縮到單位圓的圓周上呢? 圖1.  均勻散布在邊長為1的正方形內的藍點以及單位圓 (紅) 的範圍  在圖1我們產生了1000個點,以及對應單位圓圈起來的面積。這次我們並不是要估計\(\pi\),而是要找出哪些點分佈在紅色的圓周上,亦即產生出散布在單位圓上的點。你可能想說可以用以下的方式來實現 (以二維的case為例): 產生\(N\)的點的座標分別在\((-1,1)\)之間 \(x={\tt rand}(-1,1)\) \(y={\tt rand}(-1,1)\) \({\tt if}~(x^2+y^2=1):~{\tt return}~(x,y)\) \({\tt else:~return}~0\) 上述的算法在第2步就是去篩出隨機生成的\(N\)個點中,哪些是符合半徑等於\(1\)的狀況 (即這些點在紅色的圓上),如果不等於\(1\)則將之刪去。我不做計算,但我可以打包票,上面的算法接受率 (acceptance rate) 是超級低,生成數以萬計的點可能都不見得有一個是符合半徑等於\(1\)的條件。  但換個角度來看生成的點,點座標\((x,y)\)其實也代表著向量,既然單位圓有半徑為\(1\),我們只要將每個點的向量都做歸一化的動作,我們就可以得到單位向量\((x^\prime,y^\prime)\) (即向量距原點的長度為\(1\)),所以可以保證歸一化後的單位向量一定會在圖1的紅色圓上,而未歸一化前的點座標則用來決定在極座標上的角度\(\phi\),即 \[ \phi=\tan^{-1}\left(\frac{y}{x}\right). \] 示意圖如圖2。這樣不僅保證所有的點都在單位圓上,又由於點是隨機生成的,在圓周上的位置\(\phi\)也是在\(0\)和\(2\pi\)之間均勻分佈。 圖2. 歸...

IA1 - 布豐投針,圓周率\(\pi\)的估計

圖片
Mathematica®生成400根針的隨機投放示意圖 布豐投針問題 (Buffon's needle problem) 是這樣敘述的: 有一房間的地板由數片平行並列的木板所構成,板與板之間溝紋的距離皆相同,如上圖棕色隔線的間距。試問隨機扔一針,該針觸及溝紋的機率為何? 這項問題主要由18世紀的布豐伯爵所提出的,並由義大利數學家Mario Lazzarini在1901年實地進行拋針實驗,共拋出了3,408根針(🥴)並計算出\(\pi\)的近似值約為\(355/113\),參見 維基百科 。在我們繼續追本溯源這個歷史問題以前,讓我們先看看當代最常見估計圓周率的方法。 一、灑點估計法 這個方法概念非常的簡單,我們將產生\(N\)個座標\((x,y)\)隨機散佈在\(0\)和\(1\)之間的點,所這些點將被包在以原點為中心,邊長為2的正方形的第一象限中 (first quarter)。接著我們計算出每個點距原點的長度\(d=\sqrt{x^2+y^2}\),若\(d<1\)則紀錄,反之則丟棄,最後我們紀錄了\(n\)個點,這些點都被包在以原點為中心,半徑為\(1\)的\(1/4\)圓內,見圖1。 圖1. 在第一象限中藍色正方形的面積為\(1\)而紅色區域的面積為圓的\(1/4\)  現在我們相信這些處於單位圓內點的數量和全部點的數量的比值應等於那\(1/4\)圓面積和正方形面積之比,即 \[ \frac{n}{N}=\frac{A_圓}{A_方}.\tag{1} \] 而\(A_正 = 1^2\)和\(A_圓=\pi 1^2/4\),所以我們便得到了 \[ \pi = 4\frac{n}{N},\tag{2} \] 正是\(\pi\)的估計。我們預測若\(N\)愈大,則其估計值愈接近真正的圓周率。另外\(\pi\)可以解析表示成 \[ \pi = \int_{-1}^{1}\frac{dx}{\sqrt{1-x^2}}. \tag{3} \] 上式是由現代分析學之父卡爾·魏爾斯特拉斯 (Karl T. W. Weierstraß, 1815-1897) 在1841年所定義的。  我們接著可以透過 rand 來實現這個簡單的蒙地卡羅模擬,腳本如下: import numpy as np import matp...