國立交通大學
電子物理系所
碩士論文
利用反摺積運算求得點擴散函數進行
光學顯微鏡影像還原術
Optical microscope image restoration via measuring the point
spread function by deconvolution method
研究生：羅仁宏
指導教授：吳光雄 教授
利用反摺積求得點擴散函數進行光學顯微鏡影像還原術
Optical microscope image restoration via measuring the point spread function by deconvolution method
研究生：羅仁宏 Student：Ren Hong Luo 指導教授：吳光雄 Advisor：Kaung Hsiung Wu
國 立 交 通 大 學 電 子 物 理 系
碩 士 論 文
A Thesis
Submitted to Department of Electrophysics College of Science National Chiao Tung University
in partial Fulfillment of the Requirements for the Degree of Master
in
Electrophysics
June 2012
Hsinchu, Taiwan, Republic of China
i
利用反摺積運算求得點擴散函數進行光學顯微鏡影像還原術
研究生：羅仁宏 指導教授：吳光雄 國立交通大學 電子物理系所 中文摘要 本研究探討在不改變光學顯微鏡系統下，求得更接近樣品訊號，且突 破繞射極限的反摺積顯微術。原理為將光學顯微鏡量測之影像視為點擴散 函數(PSF)與樣品訊號之摺積，在已知樣品訊號的情況下，可透過反摺積運 算求得 PSF。本論文使用兩種方式定義光學顯微鏡系統的 PSF，其一如前述 以反摺積運算解得 PSF，相異於多數文獻使用樣品尺度小於繞射極限的顯微 鏡影像作為 PSF。經過將定義之 PSF 代入詳細推導並修正之最大期望演算 法與最大似然估計(EM-MLE)後可求得樣品實際訊號。比較此訊號可知，以 我們的方法解得之 PSF 代入 EM-MLE 後，在形貌與演算速度上都較多數文 獻使用的方法優異。最後以模擬確認其差異確實來自於取得 PSF 方法的不 同。另外本研究亦指出系統雜訊為帕松雜訊，並進一步討論雜訊對於反摺 積顯微術的影響。ii
Optical microscope image restoration via measuring the point spread
function by deconvolution method
Student：Ren Hong Luo Advisor：Kaung Hsiung Wu
National Chiao Tung University Institute of Electrophysics
Abstract
In our research, we use deconvolution microscopy system to obtain sample signal which is closer to the original sample morphology and non-diffraction limited in the maintenance of microscope optical system. The principle of
deconvolution microscopy is consider microscopy image as a convolution result of point spread function (PSF) and sample signal. In the case of known sample signals, PSF can be obtained by deconvolution algorithm. In this experiment, we used two method to define PSF of microscope optical system. First we bring PSF into expectation–maximization algorithm and maximum-likelihood
estimation(EM-MLE) after derivation the actual sample signal can be obtained. The second method is directly use the microscope image of the sample which is less than the diffraction limit as PSF. Result shows that not only morphology but also algorithmic speed used in the first method are both superior to the second method. Finally, we use computer simulation to confirm the result of the
difference really dependent on two method we use to obtain the sample signal. In addition, in our research also find out the noise in our system can be describe as Poisson noise and further discuss the influence to the deconvolution
iii
致謝
非常感謝吳光雄教授與劉佳明教授毫無保留的指導，兩位教授在硬體 設備、樣品、研究方向等，都非常有耐心地給予許許多多的幫助。也謝謝 陳洋元博士與許世英教授幫忙製作樣品，因為這些人的幫助，使得這篇碩 論得以順利完成。另外謝謝父親、母親與兄長從小到大的照顧與教誨，使 我有幸攻讀研究所與面對問題。 接著謝謝小趴給我實驗建議、儀器支援與陪我戰，謝謝JGB 淡定哥讓 我們超歡樂，謝謝豆姐與淑雯姐嘻嘻哈哈，謝謝聖仁哥開導我，使我不白 目，謝謝凱耀哥嗆某某，謝謝GG 跟小武關心，謝謝呆做實驗的陪伴，謝謝 小白碩一拉我讀書讓我得以歐趴，謝謝外國人家瑄的英文支援，謝謝曼麗 幫我剖析人格，與最後關頭當服務生幫買，謝謝凱迪陪喝，謝謝凱力常約 投飲料與幫掃圖，謝謝蔡婷與薇薇常常在走廊遇到可以扯兩句，也謝謝學 弟妹有時候聽我講垃圾話，最後謝謝強者我同學 BOBO 哥，常常請我吃飯 還陪我喇迪賽，並且開導我如何看事，謝謝女朋友恬妤陪我走過這一段風 風雨雨，也謝謝國陞耍白爛，另外謝謝仕賢學長常常逼我在考試前做實驗， 謝謝龍貓在空中飛。iv
內文目錄
中文摘要...i 英文摘要...ii 致謝...iii 圖目錄 ... vi 第一章 導論 ... 1 第二章 光學顯微鏡系統與成像原理 ... 5 2-1 光學顯微鏡成像原理 ... 5 2-2 光學顯微鏡系統 ... 6 2-2-1 光源模組 ... 7 2-2-2 樣品載台模組 ... 7 2-2-3 顯微鏡訊號偵測模組 ... 8 2-2-4 儀器控制與訊號處理模組 ... 9 2-3 顯微鏡影像 ... 11 第三章 影像重構 ... 13 3-1 線性非時變系統 ... 13 3-2 空間截止頻率 ... 15v 3-3 傅立葉轉換 ... 16 3-3-1 離散傅立葉轉換 ... 16 3-3-2 快速傅立葉轉換 ... 18 3-4 最大期望演算法與最大似然估計 ... 20 3-4-1 最大似然估計(MLE) ... 20 3-4-2 最大期望演算法(EM) ... 21 3-4-3 EM-MLE ... 21 3-4-4 微分式詳細推導 ... 26 3-5 影像重構程式 ... 28 第四章 分析結果與討論 ... 30 4-1 樣品 ... 30 4-2 實驗與分析結果 ... 31 4-3 分析結果討論 ... 39 4-4 雜訊影響 ... 40 4-5 雜訊分析 ... 43 4-6 迭代收斂速度比較 ... 46 4-7 不同 PSF 經 EM-MLE 結果比較 ... 48 第五章 結論與未來展望 ... 50 參考文獻 ... 51
vi
圖目錄
圖 1-1 Blind deconvolution 反解前(a)及反解後(b)比較圖[7] ... 3
圖 2-1 共焦顯微鏡原理示意圖 ... 5 圖 2-2 多功能顯微鏡光路圖 ... 6 圖 2-3 多功能顯微鏡功能圖 ... 7 圖 2-4 SPAD 工作示意圖 ... 9 圖 2-5 掃描模式示意圖 ... 10 圖 2-6 觸發訊號圖 ... 11 圖 2-7 SEM 與廣場顯微鏡比較圖 ... 11 圖 2-8 SEM 與掃描顯微鏡比較圖 ... 12 圖 2-9 螢光共焦顯微鏡影像圖 ... 12 圖 3-1 算符結果圖[15] ... 18 圖 3-2 FFT 蝶狀圖 ... 19 圖 3-3
α
與 object 關係圖 ... 24 圖 3-4 迭代過程關係圖 ... 25 圖 3-5 影像重構程式流程圖 ... 28vii 圖 4-1 金相顯微鏡樣品影像圖 ... 30 圖 4-2 掃描顯微鏡影像剖面圖 ... 31 圖 4-3 樣品(150 nm)SEM 圖 ... 32 圖 4-4 手動輸入樣品形貌圖 ... 32 圖 4-5 反解 PSF 結果圖 ... 33 圖 4-6 EM-MLE 與原始影像比較圖 ... 34 圖 4-7 樣品(1040 nm)SEM 圖 ... 35 圖 4-8 結構(1040 nm)經 EM-MLE 與原始影像比較圖 ... 35 圖 4-9 樣品(180 nm)SEM 圖 ... 36 圖 4-10 結構(180 nm)經 EM-MLE 與原始影像比較圖 ... 36 圖 4-11 樣品(1000 nm)SEM 圖 ... 37 圖 4-12 結構(1000 nm)經 EM-MLE 與原始影像比較圖 ... 37 圖 4-13 樣品(雙結構)SEM 圖 ... 38 圖 4-14 結構(雙結構)經 EM-MLE 與原始影像比較圖 ... 38 圖 4-15 1040 nm 反解結果和 PSF 摺積運算與原始影像比較圖 ... 39 圖 4-16 雙結構和 PSF 摺積運算與原始影像比較圖 ... 39 圖 4-17 1040 nm 結果頻譜振幅圖 ... 40 圖 4-18 PSF 與 Object 頻譜振幅比較圖 ... 42
viii 圖 4-19 反解 PSF 受雜訊影響結果圖 ... 42 圖 4-20 測量雜訊位置圖 ... 43 圖 4-21 (A)點取樣次數與訊號關係圖 ... 44 圖 4-22 (A)點雜訊分佈圖 ... 44 圖 4-23 平均與變異數比較圖 ... 46 圖 4-24 迭代次數 k 與 [ ]
δ
k 的關係圖 ... 47 圖 4-25 方法一與方法三反解 1040 nm 結構比較圖 ... 48 圖 4-26 模擬方法一與方法三反解 1040 nm 結構比較圖 ... 491
第一章 導論
人類對於肉眼無法看見的細微結構，總是充滿著強烈好奇心，因此顯 微鏡不論在醫學、物理、材料、生物等各領域都扮演相當重要的角色，從 歷史最悠久的光學顯微鏡到近代根據 Louis de Broglie 物質波理論，以高能 電子束成像的電子顯微鏡(Electron microscope)；或是使用量子力學中的穿 隧效應，以穿隧電流大小來判斷探針與樣品距離的掃描式穿隧探針顯微鏡 (STM, Scanning tunneling microscope)，都證明了顯微技術因應科學研究而不 斷地被創新且蓬勃發展。 雖然光學顯微鏡解析度不如電子顯微鏡的奈米等級，更無法達到 STM 原子級的解析度，但是電子顯微鏡幾乎都需要將樣品至於真空中，而 STM 只能觀察樣品的表面結構，使得光學顯微鏡因非破壞量測、穿透深度可達 數釐米等優點，目前還是被廣泛使用在各領域中，並朝著得到更佳空間解 析度、時間解析度、影像對比度、樣品穿透深度邁進，我們的研究則著重 於空間解析度部份做改進。 光學顯微鏡解析度的極限在 1884 年由德國科學家 Ernest Abbe 提出， 根據 Rayleigh criterion，橫向解析度在圓形光瞳、單色光等條件，且不考慮 像差的情況下[1] ，為 0.61 乘光源波長再除以數值孔徑(NA, Numerical aperture)。所以提高解析度直接的做法，就是改用短波長光源，或是增加透 鏡的數值孔徑。但目前這兩個方法都已經達到了極限，當光源波長達到紫2 外線波段時，就會對生物樣品染色體造成破壞；而數值孔徑中鏡口角與介 質折射率兩個參數，也都已經因為製程技術等原因達到極限，因此過去科 學家一直努力找尋其他方法增加空間解析度。 有一種增加空間解析度的方法，是使用 J. W. Goodman 在 1968 年提出 的光學顯微鏡成像原理[1] ，其原理為樣品訊號與點擴散函數(PSF, Point spread function)做摺積(Convolution)運算即為光學顯微鏡成像，進而發展出 反摺積演算法(Deconvolution Algorithms)，可在不需要改變光學顯微鏡系統 的情況下演算出比較好的影像[2] ，最早的反摺積演算法為 D. A. Agard [3] 在 1984 年提出最小相鄰反摺積(Nearest-Neighbor Deconvolution)，此演算法 假設光學顯微鏡焦平面影像的訊號，一部分其實來自焦平面上方或下方結 構，也就是縱向空間解析度不足造成影像模糊，因此將影像扣除上方或下 方結構與 PSF 摺積結果，可以得到更清晰且縱軸解析度更高的結果。
在 D.A. Agard 團隊獲得成功之後，A. Erhardt [4] 等人改進了最小相鄰 反摺積，將顯微鏡三維的 PSF 以傅立葉光學理論推導出，再使用其 PSF 反 解出顯微鏡三維清晰且高解析度影像，此方法稱為 Inverse filtering。
此時在數學領域中，A. P. Dempster[5] 等人將最大期望演算法(EM, Expectation maximization algorithm)與最大似然估計(MLE, Maximum
likelihood estimation)的方法結合起來，在只知道部分實驗數據的情況下， EM-MLE 可以估計出更合理的未知實驗數據。之後 L. A. Shepp 與 Y. Vardi[6]
3
在假設系統雜訊為帕松雜訊的情況下成功將 EM-MLE 用於斷層掃描影像解 析上，開啟 EM-MLE 還原影像的序幕。
目前最常被使用在光學顯微鏡的反摺積演算法為盲蔽反摺積演算法 (Blind deconvolution)，最早由 T. J. Holmes[7] 等人成功實現，其步驟為：
1.給定一個初始的 PSF。 2.用 PSF 與顯微鏡影像，以 EM-MLE 演算樣品原始樣貌。 3.用第 2 步驟求得的樣品原始樣貌與顯微鏡影像，以 EM-MLE 演算比步 驟 1 更合理的 PSF。 4.以步驟 3 求得 PSF 取代初始 PSF，回到步驟 2。 盲蔽反摺積演算法可大幅度提升光學顯微鏡影像清晰度，如圖 1-1。
圖 1-1 Blind deconvolution 反解前(a)及反解後(b)比較圖[7]
幾乎所有的反摺積演算法，第一步皆由透過實驗或是理論求得 PSF， 實驗求得方法為觀察一個尺度小於繞射極限的微小樣品，將取得影像直接 當成 PSF[8] ，但微小樣品難製備，且尺寸太小會造成影像訊號的訊雜比太
4 小，如果樣品尺寸太大又會造成影像與 PSF 差異過大，另外理論求得 PSF 方法是將孔徑值(NA)、光源波長等光學參數代入傅立葉光學計算[4] ，但可 惜無法考慮透鏡之球面像差等參數，將造成後續影像處理上困難。 本篇論文希望研究出一個新的方法，可以在不需要製作微小樣品的情 況下，以實驗複合理論的方法求得更合理的 PSF，進而使用此 PSF，用於 EM-MLE 回解樣品結構，除重現文獻中突破繞射極限的結果，更期望求出 更合理的樣品實際形貌。
5
第二章 光學顯微鏡系統與成像原理
本章將依序詳述論文中使用的光學顯微鏡原理、系統與顯微鏡影像。2-1 光學顯微鏡成像原理
光學顯微鏡為一相當廣泛的統稱，只要以光訊號放大樣品的影像，幾 乎都可稱為光學顯微鏡[9] ，我們在這裡介紹 Confocal 與 Conventional scanning microscopy。Confocal scanning microscopy[10] 示意圖如圖 2-1，其特點在光偵測器 前加一個針孔(pinhole)，此針孔可以將樣品焦平面下方與上方反射光濾掉， 只讓焦平面的光進入偵測器，大幅提升樣品的 Z 軸向解析度。 Sample Light Source Detector z x 圖 2-1 共焦顯微鏡原理示意圖
如將偵測器前的針孔移除，即為 Conventional scanning microscopy，其與共 焦顯微鏡最大不同在於縱軸解析度[11] 。兩種顯微鏡在同一時間點，都只 能測量樣品上小範圍的光強度訊號，因此需要移動樣品掃描，做出樣品座 標對光強度關係圖才可以得到顯微鏡影像。且一般都設計成使用同一個物
6
鏡來激發樣品與接受樣品訊號光，以降低光學系統的複雜度。
2-2 光學顯微鏡系統
本實驗所使用的顯微鏡，為實驗室自行開發的多功能顯微鏡系統(圖 2-2)，其四個功能分別為廣場顯微鏡(Wide field microscope)、掃描顯微鏡 (Conventional scanning microscope)、共焦顯微鏡(Confocal microscope)、共 焦螢光顯微鏡(Confocal fluorescence microscope)。另外多功能顯微鏡依其工 作原理可分成四個部份(圖 2-3)，分別為光源模組、樣品載台模組、顯微鏡 訊號偵測模組、儀器控制與訊號處理模組，本章將詳述各模組的工作原理。 ND stage PD L1 Laser OBJ1 Tstage Sample OBJ 2 PIstage DIM M M M M M Side view
Laser: CW diode Laser ND: ND filter
OBJ: Objective lens 1: 10x NA0.25 2: 100x NA0.80 3: 10x NA0.25 PH: Pinhole L1: Lens f=50mm M: Mirror I: Iris BSC:Beamsplitter cube B:Beamsplitter DIM:Dichroic Mirror PIstage: PI xyz stage Tstage: Tanlian x-y stage PH M BSC x y M BSC ND OBJ3 OBJ3 y z SPAD SPAD CCD ND I B Top view 圖 2-2 多功能顯微鏡光路圖
7 y z ND stage PD L1 Laser OBJ1 Tstage Sample OBJ 2 PIstage DIM M M M M M Side view PH M BSC x y M BSC ND OBJ3 OBJ3 y z SPAD SPAD CCD ND I B Top view Computer Counter card 光源系統 樣品載台系統 顯微鏡訊號偵測系統 儀器控制與 訊號處理系統 Multimeter GPIB DIO BNC RS-232USB RS-232 PCI BNC 圖 2-3 多功能顯微鏡功能圖
2-2-1 光源模組
光源模組功能為照明與激發樣品，基於雷射光(Laser)具有單色光、發散 角小、光斑分佈與理論符合等特性，選用 Cobolt 公司發展的光強度高穩定 雷射 Cobolt BluesTM，其波長 473nm，並在雷射出光後架設一組空間濾波器(Spatial filter)，以修正雷射光斑(Laser spot)的形狀與調整平行光，最後通過 由 OLYMPUS 研發長工作距離 100 倍物鏡 LMPlanFLN，孔徑值 0.8，局部 照明或激發樣品。
2-2-2 樣品載台模組
樣品載台模組由一個粗調二維平台，與一個細調二維平台堆疊而成， 細調平台為德國的 Physik Instrumente 公司所研發三軸壓電位移平台(PI
8
stage) P-563.3CD，其內建光學尺回饋平台目前位置，可達位置平方平均數 (Root mean square)正負 2 奈米的解析度。而粗調平台為台灣坦聯公司所發展 的兩軸位移平台(Tanlian stage)MR-IX71，使用步進馬達，可進行行程 20mm 長距離移動。
2-2-3 顯微鏡訊號偵測模組
訊號偵測模組根據顯微鏡功能的不同，偵測器也會有所不同。廣場顯 微鏡使用的是 Edmund 公司研發的感光偶合元件(CCD)，型號 NT63-848， 掃描顯微鏡使用 HINDS Instruments 發展的光導(Photoconductive)光電轉換 器 DET-200-002，再使用 Keithley-2000 數位電表(Multimeter)，將類比電壓 訊號轉成數位電壓訊號。而共焦顯微鏡與共焦螢光顯微鏡使用的偵測器較 為特別，為 MPD 公司發展的 PDF Series 單光子計數器(Single photon avalanche diodes)，其原理是利用 PN 二極體加以逆向偏壓[12] ，當逆向壓 大小到達崩潰電壓會使二極體導通，但如果施加逆向偏壓的時間在奈秒的 等級，二極體不會馬上有崩潰電流，而是進入一種稱之蓋格模式(Geiger mode)的暫態，如圖 2-4 開路，此模式極為不穩定，如有數光子進入二極體， 光電效應即可打破暫態平衡造成崩潰電流的產生，如圖 2-4 閉路，光強度越 強，產生崩潰電流的機率越高，再配合被動抑制電路(Passive quenching circuit)，將以導通的二極體再次回到開路。簡而言之，此偵測器會在奈秒區 間偵測有無光訊號，偵測到光訊號機率與光訊號本身強度成正比。
9 I V 開路 閉路 VBD Linear amplifier
mode Geiger mode
APD SPAD
圖 2-4 SPAD 工作示意圖
另外 PDF Series 搭載孔徑 62.5 μm 漸變式折射多模光纖(Multimode graded index)，其光纖頭可以當成共焦顯微鏡中的針孔[13] ，因此將物鏡收 集到的樣品反射光以 10 倍物鏡，孔徑值 0.25，耦合至光纖中即可完成共焦 顯微鏡，如在光纖耦合前加裝一個帶通濾光片(Band pass filter)，除去非樣 品螢光波長的光，即為共焦螢光顯微鏡。
2-2-4 儀器控制與訊號處理模組
儀器控制模組為自動測量模組，程式以 LabviewTM程式碼撰寫而成，根 據顯微鏡取得圖像的方法，可分成兩種模式，第一種稱為步進模式，步驟 為擷取來自樣品的反射光強度，移動樣品到下一目標，再擷取新訊號，重 複上述步驟即可以掃描出二維矩陣，其優點為擷取訊號時間較長，可以得 到較高的訊雜比。 第二種為掃描模式，此模式會控制樣品，向掃描軸以等速度運動移動， 再以固定週期，瞬間擷取從樣品反射光訊號(圖 2-5)，即可得到一維光學10
Time Position
samplesamplesamplesample
Gate in signal
圖 2-5 掃描模式示意圖
影像，更換掃描起始點後進行重複掃描，即可得到二維圖。但因儀器控制
時間尺度小於電腦可直接控制的 10-3秒等級，因此需利用 PI stage 控制器產
生觸發訊號(Trigger output signal)，如圖 2-6，輸入由美商國家儀器公司 (National Instruments)發展的 PCI-6602 計數卡(Counter card)，產生的 TTL 觸 發訊號(Gate in signal)，其訊號高準位時間長短可調整。將此訊號輸入 MPD，可使 MPD 在 Gate in signal 高準位時啟動。如在 MPD 啟動期間有光 訊號，MPD 就會不停輸出脈衝訊號(SPAD signal)，只要以計數卡計算在 Gate in signal 高準位期間內 SPAD signal 的脈衝個數就可以知道光強度。
11 Gate in signal Trigger output signal SPAD signal Light ON OFF Time 圖 2-6 觸發訊號圖 掃描模式好處是可以大幅提升掃描的速度，掃描模式取6000 點所花時 間約比步進模式快約一百倍，但訊雜比相較於步進模式不佳，因此需根據 實驗需求的不同，選擇不同的模式。
2-3 顯微鏡影像
本節將上述顯微鏡的影像依序列出，以了解不同顯微鏡定性上有何異 同。廣場顯微鏡結果如圖 2-7 右，與左邊掃描式電子顯微鏡(SEM)圖比較可 知解析度小於 500 nm(感謝交通大學材料所廖晨宏學長提供樣品)。 圖 2-7 SEM 與廣場顯微鏡比較圖12 掃描顯微鏡在步進模式結果如圖 2-8 右，右下角黑線為 4 μｍ比例尺， 左為 SEM 結果(感謝交通大學電物系許世英老師提供樣品)。 圖 2-8 SEM 與掃描顯微鏡比較圖 共焦顯微鏡與螢光共焦顯微鏡系統架構相當相似，在此只列出後者的 結果。螢光共焦顯微鏡結果如圖 2-9，使用掃描模式，左下角白線為 500 nm， 樣品為 Sigma Aldrich 公司所生產的乳膠粒，大小約為 100 nm。 圖 2-9 螢光共焦顯微鏡影像圖
13
第三章 影像重構
由 J. Goodman 的著作[1] 與 Agard, D. A.等人的研究[14] ，我們可知光 學顯微鏡，可近似為一個線性非時變系統(Linear shift invariant system)，本 章將詳述在此近似下如何將顯微鏡得到的影像，反解出更接近樣品實際形 貌的影像，最後說明反解過程所使用到的傅立葉轉換。
3-1 線性非時變系統
如一個系統，輸入 O (x)，輸出 I (y)且滿足：[
]
( ) ( ) I y =T O x (3.1) 其中T 為一個函數，如滿足：[
1( ) 2( )]
[
1( )]
[
2( )]
T aO x +bO x =aT O x +bT O x (3.2) 其中 a、b 各為一常數，則此系統具有線性(Linear)，另外若滿足：[
]
( ) ( ) I y+α =T O x+α (3.3) 則此系統具有非時變(Time-invariant system)，如系統同時滿足線性(3.2)式， 與非時變(3.3)式兩個特性，就是線性非時變系統(簡稱 LSI 系統)。LSI 系統可以簡單表示T，方法為將脈衝函數(通常選用 Kronecker delta
function)輸入系統，將會得到有特別意義的輸出，此輸出稱為脈衝響應 (Impulse response)：
[
]
( ) ( )
14
脈衝響應與系統輸入，輸出可以表示成摺積運算(Convolution)：
I = ⊗ P O (3.5)
如將光學顯微鏡系統導入上述的模型中，可了解光學顯微鏡的成像 (image)，即可看成點擴散函數(PSF, Point spread function)與測量的樣品訊號 (object)做摺積運算：
[ ]
[ ]
[ ]
psf x
⊗
object x
=
image x
(3.6) 其中psf x
[ ]
為點狀樣品通過光學系統繞射的結果，如得知光學顯微鏡的 PSF，依照線性 LSI 系統的特性與(3.7)式的 Convolution 定理：{
A B}
{ } { }
A B ℑ ⊗ =ℑ × ℑ (3.7) 即可使用(3.8)式推出原始物體(object)傅立葉轉換結果，經由傅立葉逆轉換 即可得原始物體(object)形貌。{
[ ]
}
{
{
[ ]}
}
[ ]
image xobject x
psf x
ℑ ℑ = ℑ (3.8) 欲得知光學顯微鏡 PSF，最簡單且直接方法是用一個如數學定義般點 狀樣品，放入光學顯微鏡中掃描，將掃描到的成像當成 PSF。但在實驗上 製作如數學定義一般點狀樣品有技術上困難，且其顯微鏡訊號會很微弱， 訊雜比(Signal-to-noise ratio)大幅下降，將造成後續影像處理上困難，因此需 用另一方法得到 PSF。 在此論文中使用的方法為使用兩次(3.6)式，首先使用一個已知其形貌15 的樣品，在此稱為校正樣品(Calibration sample)，經由光學顯微系統量測， 用數學式描述如(3.9)式。 Calibration Calibration
psf
⊗
sample
=
image
(3.9) 將得到的影像(Calibration image)與校正樣品做反摺積運算可得 PSF，此方法 不需製作如數學定義般點狀樣品。3-2 空間截止頻率
一個無像差、單色光且圓形光瞳的 PSF 可表示成下式：(
)
2 1 0 02
2
/
( )
2
/
J
Nr
psf r
Nr
p
λ
p
λ
=
(3.10)上式中，J1為第一類貝索函數(Bessel functions of the first kind)，N 為數值孔
徑(NA, Numerical aperture)，
λ
0為光源的波長，r
為在像平面上的一點到光 軸的距離。將(3.10)式傅立葉轉換，可以得到光傳遞函數(OTF, Optical transfer function)： 0 1, 1/ 2 ( ) ( ) ( ), ( ) 4 0, 1/ 2 x s OTF s f s x N xλ
< = = > ∏
∏
(3.11) 可知 OTF 在頻率s
一定值以上都等於零，此頻率稱為空間截止頻率(Spatial cut-off frequency)，其物理意義為，以頻譜空間討論繞射極限問題： 0 2 cut off N frequency λ − = (3.12)16
3-3 傅立葉轉換
我們使用快速傅立葉轉換(Fast Fourier Transform)處理反褶積運算
(deconvolution)，其在處理離散數據，具有簡單化反褶積運算與較離散傅立 葉轉換(Discrete Fourier Transform)計算快等優點，本節將由基本傅立葉轉換 出發，延伸至可討論離散數據的離散傅立葉轉換。
3-3-1 離散傅立葉轉換
傅立葉轉換[15] 可將在實空間的數據 f (x)，轉換到頻率空間 F (s)，其在數 據連續且無周期性的條件下(Function on R)，傅立葉正、逆轉換表示如式 (3.13)、(3.14)： 2 ( ) ( ) isx s f x ∞ F s e p ds =−∞ =∫
(3.13) 2 ( ) ( ) isx x F s ∞ f x e− p dx =−∞ =∫
(3.14) 但數值分析無法處理積分區間無窮大，所以需引入週期性，即數據實 空間訊號具有 p 的週期性使 f (x + p)=f (x)，滿足上述條件即可推導出對應頻 譜空間將不連續，其頻譜空間為： [ ], 0, 1, 2, 3, , F k k = ± ± ± 3 ±∞ (3.15) 中括號表示 k 為整數，且(3.13)、(3.14)可以變為： 2 / [ ] ( ) ikx p k F k e f x p ∞ =−∞ =∑
(3.16)17 2 / 0 1 [ ] p ( ) ikx p x p F k f x e− p dx = =
∫
(3.17) 再考慮數據在實空間因取樣(sample)造成的不連續性，實空間函數可以 表示成： [ ], 0, 1, 2, , f n n= ± ± 3 ±∞ (3.18) 可推導出在頻率空間將有週期性 F (s + p)=F (s)，則(3.13)、(3.14)式可以變 為： 2 / 0 [ ] p ( ) isn p s f n F s e p ds = =∫
(3.19) 2 / 1 [ ]( )
isn p n f n e pF s
∞ − p =−∞=
∑
(3.20) 因此如數據同時具有週期與不連續性，在頻率空間也會同時具有不連 續與週期性，兩空間週期皆為N ，或稱為窗口大小(Window size)，(3.13)、 (3.14)式可以簡化為只討論一個窗口，即稱為離散傅立葉轉換(DFT, Discrete Fourier transform)： 1 2 / 0 [ ][ ]
N ikn N k F k ef n
p − ==
∑
(3.21) 1 2 / 0 1 [ ][ ]
N ikn N n f n e NF k
p − − ==
∑
(3.22)18
3-3-2 快速傅立葉轉換
快速傅立葉轉換(FFT)利用了離散傅立葉的兩個性質，假設一組離散週 期數據(a, b, c, d)傅立葉轉換後為(A, B, C, D)，DFT 具有 Translation rule：
2 3 2 /
, , , )
, ,
,
)
(
d a b c
(
A w B w C w D
,
w e
− pi Nℑ
=
=
(3.23)上式中 N=4 代表窗口大小，另外 DFT 還具有 Zero packing rule：
( , 0, , 0, , 0, , 0)
a
b
c
d
( , , , , , , ,
A B C D A B C D
) / 2
ℑ
=
(3.24) 將上述兩規則以算符ϒ Ε−1, −1與Ζ ℜ2, 2表示： 1 − ϒ(a, b, c, d)
ℑ(d, a, b, c)
ℑ(A, B, C, D)
1 − Ε(A, wB, w
2C, w
3D)
2 Ζ(a, b, c, d)
ℑ(a, 0, b, 0, c, 0, d, 0)
ℑ(A, B, C, D)
2 ℜ(A, B, C, D, A, B, C, D)
圖 3-1 算符結果圖[15] 再利用單一數傅立葉轉換不變：( )
α
( )
αℑ
=
(3.25) 可將(a, b, c, d)的傅立葉正轉換，依循圖 3-2 FFT 蝶狀圖運算，此類運算流 程稱為快速傅立葉正轉換。19 1 − Ε 2 ℜ 1 − Ε 1 − Ε 2 ℜ 2 ℜ
Time domain Frequency domain
( )a 2 Ζ ( )c ( , 0)c 1 − ϒ 2 Ζ ( , 0)a (0, )c ( , )a c ( )b 2 Ζ ( )d ( , 0)d 1 − ϒ 2 Ζ ( , 0)b (0, )d ( , )b d 2 Ζ ( , 0, , 0)b d 1 − ϒ 2 Ζ ( , 0, , 0)a c (0, , 0, )b d add add add ( , , , )a b c d ( )a ( )c 1 ( , ) 2 c c 1 ( , ) 2 a a 1 ( , ) 2 c ωc 1 ( , ) 2 a+c a+ωc ( )b ( )d 1 ( , ) 2 d d 1 ( , ) 2 b b 1 ( , ) 2 d ωd 1 ( , ) 2 b+d b+ωd 1 ( , , , ) 4 b+d b+ωd b+d b+ωd 1 (*,*,*,*) 4 add add add ( , , ,A B C D) 2 ℜ 2 ℜ 2 ℜ 1 (*,*,*,*) 4 圖 3-2 FFT 蝶狀圖 同理，傅立葉逆轉換使用 Modulation rule： 2 3 2 / ) )
( ,
a w b w c w d
,
,
( , , ,
B C D A w
,
e
− pi Nℑ
=
=
(3.26) 與 Summation rule： ((a c,b d+ ) / 2) ( , )A C ℑ + = (3.27) 也可寫出快速傅立葉逆轉換，上述正逆轉換即稱為快速傅立葉轉換。20 由上述討論可知 DFT 與 FFT 求出的答案一致，且在 N 為以 2 為底的指 數情況下，執行 DFT 所需複數加法量為N2，而 FFT 的複數加法量為 2 log N N ，這也是 FFT 較 DFT 計算速度快的原因。
3-4 最大期望演算法與最大似然估計
要了解如何將 EM-MLE 用在光學顯微鏡影像還原上，需分別了解 MLE 與 EM 兩種方法，再討論將顯微鏡成像如何使用 EM-MLE。3-4-1 最大似然估計(MLE)
MLE 為估計未知數據的方法，假設所討論的系統，有一個參數 p未知， 這參數影響測量出的已知數據Data，MLE 會求出當未知參數為多少時， 最有可能量到目前已知的數據，詳細的步驟如下： 1.寫出似然函數L(likelihood function)，其定義為當系統參數(Parameter) 為 p時量到Data的條件機率有多少：(
|)
L = P Data p (3.28) 2.使用微分法求似然函數最大值，似然函數對參數一次與二次偏微分，當 找到一個 p0使L對 p偏微分等於零，且二次偏微分小於或等於零，即 可知當系統參數為 p0時量到已知數據機率為最大，以數學式表示： 0 2 2 ( ) 0 0 p Max p L L p = p ∂ = ∧ ∂ ≤ ∂ ∂ (3.29)21
3-4-2 最大期望演算法(EM)
EM 用迭代法(Iterative method)尋找系統參數，使似然函數最大值，如 前段假設系統有一個未知參數會影響到測量數據，詳細步驟如下： 1.將 pk 當成系統參數，寫出似然函數L，如(3.28)式。 2.尋找一個新的參數pk+1代入步驟 1 的結果大於將 pk代入步驟 1 的結 果，表示系統在參數 pk+1的情況下，比在參數為 pk的情況下更有可能 量到已知Data：(
)
(
)
1 1 ( k ) | k ( k) | k L p + =P Data p + ≥ L p = P Data p (3.30) 3.除非 pk 收斂，否則將 pk+1取代 pk回到步驟 1。3-4-3EM-MLE
不難看出 MLE 困難之處，為如何表示系統的似然函數，與在系統的參 數較多，並且L對 p偏微分式較複雜的情況下，即使列出方程式，也未必 能求出未知參數。而 EM 則具有計算量較大，難寫出迭代方程式，與同樣 需寫出系統似然函數等問題。為解決以上問題，Y. Vardi 等人將光學系統引入[16] ，再將 EM 與 MLE 兩種方法結合起來，首先假設光學系統對測量到的訊號，會有一雜訊分佈， 其分佈符合帕松分佈(Poisson distribution)，在 4-5 雜訊分析章節會再次驗證 此假設的合理性，帕松分佈函數如下：
22 ( ) ! k P k e k λ λ − = (3.31) 此式子在描述，當系統的真實值(或期望值)為 λ 時，測量結果為k 的機率P 為多少。另外再假設樣品為 object [x]，且系統的 PSF 為 psf [x]，量到的訊 號為imagenoise( )x 。注意已考慮測量到的影像包含雜訊，因此： [ ] [ ] noise[ ] object x ⊗ psf x ≠image x (3.32) 將上述參數代入帕松分佈函數(3.31)，可得到在樣品為 object [x]的條件下， 系統的實際值： [ ] [ ] [ ] [ ] x object x psf x psf x x object x ′ ′ ′ ⊗ =
∑
− (3.33) 因此在位置x0量到imagenoise[x0]的條件機率為： 0 0 [ ] 0 [ ] [ ] [ ] [ ] [ ]! noise x image x psf x x object x x noise psf x x object x P e image x ′ ′ ′ − − ′ − ′ ′ ∑ =∑
(3.34) 可以將(3.34)式視為此系統的測量影像在x0點的似然函數：(
)
0 int ( [ 0]) [ 0] | [ 0] po x noiseL object x = P image x object x (3.35)
考慮樣品object x[ 0,x1]同時符合x0與x1量到的數據，其似然函數為個別機 率的相乘：
(
)
(
)
2 int 0 1 0 0 1 1(
[
,
])
[
] |
[
]
[ ] |
[ ]
po noise noiseL
object x x
P image
x
object x
P image
x
object x
=
×
(3.36)
23 object [x]，符合整個
x
空間中量到的數據，其似然函數為： [ ] [ ] [ ] [ ] [ ] [ ]! x image x psf x x object x x x psf x x object x L e image x ′ ′ ′ − − ′ − ′ ′ ∑ =∏
∑
(3.37) 為簡化問題，使用對數是一個遞增函數的性質，改討論似然函數取自然對 數的最大值，取代直接討論似然函數的最大值，將ln L對object x( 0)偏微 分與二次偏微分，詳細的推導將在 3-4-4 討論，在這裡先列出結果為： 通常 image [x]物理量與光強度等價，psf [x]、object [x]分別與光強度、反射 率等價，所以 image [x]、psf [x]、object [x]都不會出現小於零的情況，因此 可知(3.39)恆小於或等於零，ln L(object [x])函數為凸函數(concave function)，因此似然函數最大值發生在(3.38)等於零時，且通常 psf [x]為歸一 化函數(Normalized function)，使： 0 [ ] 1 x psf x− x =∑
(3.40) 因此似然函數最大值發生在 object [x]滿足： 0 0 0 0 ln ( [ ]) [ ] [ ] [ ] [ ] [ ] [ ] x x x L object x psf x x object x image x psf x x psf x x object x ′ ∂ = − − + ∂ − ′ ′ −∑
∑ ∑
(3.38) 2 2 0 0 2 2 0 ln ( [ ]) [ ] [ ] [ ] [ ] [ ] x x L object x image x psf x x object x psf x x object x ′ ∂ = − − ∂ ′ ′ − ∑
∑
(3.39)24 0 [ ] [ ] 1 [ ] [ ] x x image x psf x x psf x x object x ′ − = ′ ′ −
∑ ∑
(3.41) 此時不知道 object [x]等於多少，會滿足(3.41)，所以嘗試使用 EM 寫出迭代 關係式，首先由(3.41)式中令： 0 [ ] [ ] [ ] [ ] x x image x psf x x psf x x object xα
′ − = ′ ′ −∑ ∑
(3.42) 又因為(3.39)式得知： 0 object α ∂ ≤ ∂ (3.43)因此
α
與 object [x]關係為負相關(Negative correlation)，示意圖如圖 3-31
α
object Answer objectAns 圖 3-3α
與 object 關係圖 考慮某次迭代： k Ans object >object (3.44) 其對應α
的量值將介於 0 到 1 之間，將objectk 乘以α
可得到更接近 Ansobject 的objectk+1，如圖 3-4 左；同理可知，如果objectk 小於objectAns，
25 1 k object + ，如圖 3-4 右。 1
α
object Answer objectk objectkα
×
1α
object Answer objectk objectkα
×
圖 3-4 迭代過程關係圖 注意迭代的過程中，在 psf 或 image 與實際情況誤差太多的情況下，可能發生objectk 不斷的在objectAns左右跳來跳去，且距離objectAns越來越遠的情
況，因此迭代為注意每次迭代造成objectk 的變化是否漸減，需討論(3.45) 式，這裡分母選用objectk+1較接近實際
object
： 1 1 [ ] [ ] [ ] [ ] k k x k x object x object x k object x δ + + − =∑
∑
(3.45) 只要δ[ ]k 恆大於δ[k +1]，即可確定迭代過程有在逼近objectAns。 綜合以上 3-4-3 的討論可以知道 EM-MLE 最後可以簡化成一個簡單迭 代方程式： 雖然推導結果與大部分文獻[17] [18] 使用的公式不同： 0 1 0 0 [ ] [ ] [ ] [ ] [ ] [ ] k k x k x psf x xobject x object x image x
psf x x object x + ′ − = ′ ′ −
∑ ∑
(3.46)26 不同之處在於分子 psf 的座標，在我們的系統中，使用自行推導出公式較合 理，詳細內容將於 4-6 節討論。
3-4-4 微分式詳細推導
本段將推導由(3.37)得到(3.38)與(3.39)的過程，對(3.37)取自然對數可 得： [ ] [ ] ln ( ) [ ] [ ] [ ln ] ln( [ ]!){
}
x x x psf x x object x L object image x psf x x object x image x ′ ′ ′ ′ − − + = ′ ′ − −∑
∑
∑
(3.48) 因為通常 image [x]為光計數的結果，數量級約103，因此可使用斯特靈近似 (Stirling's approximation)：ln( !)
n
≈
n
ln
n n
−
(3.49) (3.48)式可變為： [ ] [ ] ln ( ) [ ] [ ] [ ln ] [ ]ln( [ ]) [ ]{
}
x x x psf x x object x L object image x psf x x object ximage x image x image x
′ ′ ′ ′ − − + = ′ ′ − − +
∑
∑
∑
(3.50) 再整理成： 0 1 0 0 [ ] [ ] [ ] [ ] [ ] [ ] k k x k x psf x xobject x object x image x
psf x x object x + ′ − = ′ ′ −
∑ ∑
(3.47)27 [ ] [ ] ln ( ) [ ] [ ] [ ] ln [ ] [ ]
{
}
x x x psf x x object x L object psf x x object x image x image x image x ′ ′ ′ ′ − − + = ′ ′ − +∑
∑
∑
(3.51) 將(3.51)式對object x[ 0]偏微分，討論更動x0處object x[ 0]的量值對似然函數 的影響： 0 0 0 0 ln ( [ ]) [ ] [ ] [ ] [ ] [ ] [ ] [ ] [ ]{
}
x x L object x sf x x p object x image x image x psf x x object x psf x x image x ′ ∂ − = − + ∂ × ′ ′ − − ∑
∑
(3.52) 整理可得(3.38)式。 將(3.38)式對object x[ ]1 偏微分： 令x0 = x1可以得到(3.39)式。 2 0 1 2 0 1 [ ] [ ] [ ] ln [ ] [ ] [ ] [ ] x x image x psf x x psf x x Lobject x object x psf x x object x
′ − − ∂ =− ∂ − ′ ′
∑
∑
(3.53)28
3-5 影像重構程式
Start End Ture False k←k+1 i = N-1? Ture i←i+1 i←0 False 1 [ ] [ ] [ ] [ ] [ ] [ ] N i i i k k N x k x psf x xobject x object x image x
psf x x object x + ′ − ← ′ ′ −
∑
∑
Input objectcSolve object through EM-MLE Solve PSF through Inverse filtering 1 0 0 [ ] if [ ] [ ] 0 if c c IMAGE k k k psf x OBJECT k k k − ℑ ≤ ← > { } { } [ ] [ ] [ ] [ ] c c c c OBJECT k object x IMAGE k image x ← ℑ ← ℑ 0 δ δ< 圖 3-5 影像重構程式流程圖
29 影像重構程式以 LabviewTM程式撰寫而成，分成以逆濾波器(Inverse filtering)反解 PSF，與以 EM-MLE 反解未知樣品兩個部分，程式流程圖如 圖 3-5，首先在程式開始後，先輸入校正樣品object xc[ ]與校正樣品顯微鏡 影像image xc[ ]，窗口大小皆為N，經過傅立葉轉換相除完成反摺積，再將 頻域空間將空間截止頻率k0以上的頻率濾除，也就是令實虛部皆為零，將 此結果逆轉換後得到 PSF。 第二部分將第一部分求得的 PSF 代入 EM-MEL，也就是(3.46)式，由兩 個迴圈所構成，第一個為 For loop，確保實空間窗口內每一個xi都會被迭代
一次，第二個為 While loop，計算(3.45)中每次疊造成objectk 的變化，如果
變化小於目標δ0，代表繼續迭代下去，產生變化量已經不在討論的範圍內，
30
第四章 分析結果與討論
4-1 樣品
為縮短測量的時間，降低光學系統上環境、機械與光源造成影像誤差， 本論文暫且將推導、實驗與樣品都簡化為一維。一維樣品希望如(3.6)式般， 在 PSF 不為零的空間內，光學顯微鏡影像只會與樣品在直角座標上一軸的 位置有關係，與另外兩軸無關。在此特別感謝中研院陳洋元博士提供符合 條件之樣品，該樣品在金相顯微鏡(物鏡孔徑值 0.9)下影像如圖 4-1，白色線 長皆為 15 μm，線寬如右圖右側標示：150 nm
180 nm
1040 nm
150 nm
180 nm
1000 nm
(5)
(1)
(2)
(3)
(4)
圖 4-1 金相顯微鏡樣品影像圖 而製作線的材料為鋁金屬，基板(圖 4-1 中灰色部分)為 ITO 玻璃基板，線與 基板都刻意選用可導電材料，使樣品可放入掃描式電子顯微鏡(SEM)獲得更 高解析度影像，以確認影像重構校正 PSF 部分中的樣品形貌。31
4-2 實驗與分析結果
樣品經掃描顯微鏡使用步進模式掃描結果如圖 4-2 所示，由左到右分別 是 1040 nm、150 nm、180 nm、1000 nm 與 150 nm 和 180 nm 雙結構，點與 點之間間距 10 nm，掃描方式由圖 4-1 上方往下方掃描，0
10
20
30
40
50
60
0.00
0.05
0.10
0.15
0.20
0.25
0.30
0.35
(5) (4) (2) (3)S
igna
l(
V
)
Position (µm)
(1) 圖 4-2 掃描顯微鏡影像剖面圖 首先使用第(2)條線 150 nm 結構解出系統的 PSF，流程如圖 3-5 中的 Solve PSF through Inverse filtering 部分，窗口大小N為 1000，結構形貌 SEM 結果如圖 4-3，手動輸入樣品形貌object xc[ ]如圖 4-4，圖中只畫出較高反射率附
近的區域，以方便觀察，另外高反射率代表樣品上鋁線，可將雷射光通過 物鏡聚集在 1040 nm 鋁線上，實驗測得反射率為 23.64%；而低反射率代表
32 樣品上基板區域，將雷射光斑聚焦於基板處，可以測得反射率為 2.056%。 圖 4-3 樣品(150 nm)SEM 圖 18.2 18.3 18.4 18.5 18.6 0 5 10 15 20 25 R ef le ctiv ity ( % ) Position (µm) 圖 4-4 手動輸入樣品形貌圖
33
將 150 nm 顯微鏡影像，與object xc[ ]代入逆濾波器運算，且將k0根據(3.12)
式設為 6
3.38 10× Hz，反解出 PSF 結果如圖 4-5 所示：
34 圖 4-5 下為放大上圖 PSF 左半邊結果，由圖中可推算出 PSF 半徑為 400 nm。 雖然(3.40)式要求，PSF 需為歸一化常數，但觀察最終迭代式(3.46)可以發 現，PSF 歸一與否不影響迭代結果。 接著進行圖 3-5 中 EM-MLE 解析樣品部分，窗口大小 N 皆為 1000，如 果每次迭代變化佔迭代結果小於 0.01%，即(3.45)式中 [ ]δ k 小於 0.0001，則 迭代結束。結果如圖 4-6 所示，黑線為反解後結果，點為原始光學顯微鏡影 像，因 150 nm 結構已經用於反解 PSF，因此沒有反解的影像。
10
20
30
40
50
0.00
0.05
0.10
0.15
0.20
0.25
0.30
0.35
0.40
S
igna
l(
V
)
Position (µm)
OM Image
EM-MLE
圖 4-6 EM-MLE 與原始影像比較圖 接著列出各段反解結果與反解前影像比較圖，各影像的收斂速度將在 4-5 節比較：35 1040 nm 結構 SEM 圖如圖 4-7，顯微鏡影像使用圖 4-2(1)號影像，反解 結果與反解前影像如圖 4-8，收斂於第 236 次迭代，反解影像半高寬(FWHM) 為 1060 nm，原始影像半高寬為 1110 nm。 圖 4-7 樣品(1040 nm)SEM 圖 3 4 5 6 7 8 9 10 11 12 13 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 S igna l ( V ) Position (µm) EM-MLE OM Image 圖 4-8 結構(1040 nm)經 EM-MLE 與原始影像比較圖
36 180 nm SEM 圖如圖 4-9，顯微鏡影像使用圖 4-2(3)號影像，反解結果 與反解前影像如圖 4-10，收斂於第 251 次迭代，反解影像半高寬為 190 nm， 原始影像半高寬為 420 nm。 圖 4-9 樣品(180 nm)SEM 圖 25 26 27 28 29 30 31 32 33 34 35 0.00 0.05 0.10 0.15 0.20 0.25 S igna l ( V ) Position (µm) EM-MLE OM Image 圖 4-10 結構(180 nm)經 EM-MLE 與原始影像比較圖
37 1000 nm SEM 圖如圖 4-11，顯微鏡影像使用圖 4-2(4)號影像，反解結果 與反解前影像如圖 4-12，收斂於第 192 次迭代，反解影像半高寬為 1050 nm， 原始影像半高寬為 1090 nm。 圖 4-11 樣品(1000 nm)SEM 圖 33 34 35 36 37 38 39 40 41 42 43 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 S igna l ( V ) Position (µm) EM-MLE OM Image 圖 4-12 結構(1000 nm)經 EM-MLE 與原始影像比較圖
38 雙結構 SEM 圖如圖 4-13，顯微鏡影像使用圖 4-2(5)號影像，反解結果 與反解前影像如圖 4-14，收斂於第 311 次迭代，反解影像半高寬為 180 nm 與 200 nm，原始影像半高寬為 410 nm 與 420 nm。 圖 4-13 樣品(雙結構)SEM 圖 45 46 47 48 49 50 51 52 53 54 55 0.00 0.05 0.10 0.15 0.20 0.25 S igna l ( V ) Position (µm) EM-MLE OM Image 圖 4-14 結構(雙結構)經 EM-MLE 與原始影像比較圖
39
4-3 分析結果討論
為確定解出影像合理性，我們將 PSF 與解出影像代入(3.6)式進行摺積 運算，理論上計算結果會等於顯微鏡影像，而實際運算結果，每一個結構 除小幅度雜訊震盪外，完全符合顯微鏡影像，這裡只列出代表的 1040 nm 與雙結構驗證結果，分別如圖 4-15、圖 4-16。 3 4 5 6 7 8 9 10 11 12 13 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 Convolution S igna l( V ) Position (µm) 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 OM Image 圖 4-15 1040 nm 反解結果和 PSF 摺積運算與原始影像比較圖 45 46 47 48 49 50 51 52 53 54 55 0.02 0.04 0.06 0.08 0.10 0.12 0.14 Convolution S igna l ( V ) Position (µm) 0.02 0.04 0.06 0.08 0.10 0.12 0.14 OM Image 圖 4-16 雙結構和 PSF 摺積運算與原始影像比較圖40 1040 nm 結構上方與結構邊緣，週期約 400 nm 震盪為影像失真，原因 為 3-4-3 所討論，因 PSF 會使高頻項遺失，如不把影像高頻項濾除會使反摺 積時在頻譜空間出現分母為零。實驗上 PSF 截止頻率有可能大於 6 3.38 10× Hz，但無法判別大於截止頻率的量值中雜訊佔了多少比例，如雜 訊所佔比例過多將造成 PSF 不合理，詳細於 4-4 節討論。 EM-MLE 可以還原小部分遺失高頻項[19] ，1040 nm 反解結果頻譜振 幅如圖 4-17，截止頻率為大於低通濾波截止頻率 6 3.38 10× Hz， 0 1x106 2x106 3x106 4x106 5x106 6x106 0 25 26 27 28 29 30 A m pl itude ( ar b. uni t) Frequency (1/m) 圖 4-17 1040 nm 結果頻譜振幅圖
4-4 雜訊影響
本節將討論大部分反摺積系統都會遇到的雜訊問題，假設系統符合摺41 積關係式(3.6)，但實驗會因雷射光源強度變化，光訊號轉電訊號會有熱雜 訊等因素，無法測量到實際
image x
[ ]
，而是測量到含有雜訊的影像，將雜 訊影像由一個簡單的式子描述： [ ][ ]
[ ]
Noise ximage
=
image x
+
Noise x
(4.1)別無選擇地將imageNoise[ ]x 當成
image x
[ ]
代入(3.6)式： [ ][ ]
[ ]
Noise xpsf x
⊗
object x
=
image
(4.2)假設
object x
[ ]
為已知，想解出psf x
[ ]
，將(4.2)式傅立葉轉換以進行反摺積運算：
{
psf}
{
object}
{
imageNoise}
ℑ × ℑ = ℑ (4.3)
假設頻譜上某一點：
{
image}
A iB, ,{
Noise}
AN iBN{
object}
C iDℑ = + ℑ = + ℑ = + (4.4) 則頻譜上此點ℑ
{
psf}
為：{
[ ]
}
{
image{
Noise}
}
objectpsf x
ℑ ℑ ℑ=
(4.5) 將(4.4)式代入(4.5)式：{
}
2 2 ( ) ( ) [ ] N N N N AC BD i BC AD A C B D i B C A D C D psf x ℑ + + − + + + − = + (4.6) 通常C和D很小時，根據(3.6)式得知 A和B也會很小，如此時雜訊也 很小，造成的影響不大；但如此時雜訊很大，根據(4.6)式ℑ{
psf}
解出結果42 較實際值大數至數萬倍，如圖 4-18，可發現圖中 PSF 頻譜中不合理突峰都 發生在 150 nm 結構頻譜振幅很小之頻率，此類放大雜訊的結果與文獻中的 敘述相吻合[20] ，而其頻譜振幅不合理突峰，將造成反解 PSF 帶有許多高 頻震盪項，如圖 4-19 所示。 0 1x107 2x107 3x107 4x107 5x107 0.0 0.2 0.4 0.6 0.8 1.0 1.2 PSF A m pl it ude ( ar b. uni t) Frequency (1/m) 0 1 2 3 4 5 Object 圖 4-18 PSF 與 Object 頻譜振幅比較圖 0 1 2 3 4 5 6 7 8 9 10 0.00 0.01 0.02 0.03
S
ig
na
l (
V
)
Position (µm)
圖 4-19 反解 PSF 受雜訊影響結果圖43
4-5 雜訊分析
由 4-4 雜訊影響相當劇烈，因此本節嘗試由實驗量測，分析出雜訊的種 類，實驗方法為停留在 1040 nm 結構上面的五個點，每個點各取一萬筆數 據，如圖 4-20 與圖 4-21，再將所測量數據中最大值與最小值之間分成 300 個區間，分別計算一萬筆數據中，有幾個數據量值符合各個區間內，將個 數與區間量值作圖即可得雜訊分佈圖，圖 4-20 中(A)點雜訊分佈圖如圖 4-22。 3 4 5 6 7 8 9 10 11 12 13 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 (E) (A) (D) (C) S igna l( V ) Position (µm) OM Image (B) 圖 4-20 測量雜訊位置圖44 0 2000 4000 6000 8000 10000 0.32 0.33 0.34 S igna l ( V )
Sampling times (count)
圖 4-21 (A)點取樣次數與訊號關係圖 0.320 0.322 0.324 0.326 0.328 0.330 0.332 0.334 0.336 0 20 40 60 80 100 120 T im es ( count ) Signal (V) 圖 4-22 (A)點雜訊分佈圖
45 觀察圖中可以懷疑雜訊為高斯分佈，或帕松分佈，根據劉佳明教授的著作 [21] ，光子進入光偵測器的過程幾乎都符合帕松過程，此過程機率分佈函 數的變異數(Variance)會與分佈函數期望值，也就是算數平均數(Arithmetic mean)相等，我們將圖 4-20 中五個點的數據都乘以同一常數後，將變異數依 循以下公式算出，
(
)
2 1 var(X) X N µ = − (4.7) 再與平均數作圖可得圖 4-23，比較可得圖 4-20 中五個點，平均量值均與變 異數相差不多，再如圖 4-22 分別畫出(B)、(C)、(D)、(E)雜訊分佈圖，都類 似圖 4-22 分佈，因此可假設整個顯微鏡系統雜訊為帕松雜訊。 需乘常數是因為，帕松分佈描述的過程是某段時間內會發生若干次事 件；對應本實驗中，為某段時間內光子有若干光子被訊號偵測模組偵測到， 但是因偵測器模組輸出為電壓訊號，造成訊號不等於被偵測到的光子數， 而是需要乘一個因數，此因數以列舉法求得為 65294.8。46 400 420 440 460 480 500 520 540 560 0 2000 4000 6000 8000 10000 12000 14000 16000 18000 20000 22000 24000 (A) (B) (C) (D) Mean Variance S igna l( V ) Position (µm) (E) 圖 4-23 平均與變異數比較圖
4-6 迭代收斂速度比較
本節將比較四種方法的收斂速度，並將四種結果分別列於表格一，四 種方法分別如下： 方法一：本實驗使用的方法，反解 PSF 再反解結構。 方法二：本實驗使用的方法，EM-MLE 改使用(3.47)式。 方法三：將 150 nm 結構的顯微鏡影像當成 PSF，EM-MLE 使用(3.46)式。 方法四：將 150 nm 結構的顯微鏡影像當成 PSF，EM-MLE 使用(3.47)式。 方法三與方法四目的在比較使用小於繞射極限樣品影像當作 PSF 與使用本 實驗方法求得有何異同，注意需將 150 nm 影像扣除 0.02 V 背景值，原因 為真實 PSF 離中心亮紋較遠處，應該會如圖 4-5 逐漸收斂到零。47 表格一 比較各結構在不同方法下收斂時間 收斂時間(次) 1040 nm 180 nm 1000 nm 150, 180 nm 方法一 236 251 192 311 方法二 244 272 205 327 方法三 351 380 284 455 方法四 360 402 289 463 經由表格一的比較，可以發現本實驗使用的方法收斂最快。雖然差異都在 數十次左右，但是在反解二維結構情況下，每次迭代計算量將龐大許多， 這時收斂次數多寡將大幅影響整體演算時間，四個方法反解出結構無明顯 差異。另外根據 3-4-3 節的討論，需確定(3.45)式的
δ
[
k
+
1]
要恆小於[ ]
k
δ
，這裡列出使用方法一迭代次數 k 與 [ ]δ
k 關係圖，由圖中可以確定每 次迭代造成變化的確越來越小，無 3-4-3 節不合理現象。50
100
150
200
250
300
0.000
0.002
0.004
0.006
0.008
0.010
1040 nm 180 nm 1000 nm 150, 180nmδ[
k]
(
ar
b.
uni
t)
k (count)
圖 4-24 迭代次數 k 與 [ ]δ
k 的關係圖48
4-7 不同 PSF 經 EM-MLE 結果比較
本節將比較前節中方法一與方法三解出結構有何異同，1040 nm 顯微鏡 影像反解出結構如圖 4-25： 304 305 306 307 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 Method 1 Method 3 S igna l( V ) Position (µm) 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 Check 圖 4-25 方法一與方法三反解 1040 nm 結構比較圖 經比較後發現方法一也較方法三結構中失真比較小，為確保失真不是由較 多次迭代所產生，將方法一在 EM-MLE 步驟迭代 351 次，並將結果畫在圖 4-25 中(圖例為 check)，觀察後可確定較大影像失真不是因方法三執行較多 次迭代造成。 為了確定上述差異，的確是由反解 PSF 與否造成，我們模擬實驗結果， 方法為依照(3.10)寫出 PSF，將 PSF 分別與 150 nm、1040 nm 結構摺積運算， 依照方法一與方法三。結果如圖 4-26 所示，方法一 EM-MLE 收斂在第 11649 次迭代，方法三收斂在第 177 次迭代，較實驗快收斂，而其他結構也都得 到相同結果。
303
304
305
306
307
0.00
0.05
0.10
0.15
0.20
0.25
0.30
0.35
0.40
S
igna
l(
V
)
Position (
µm)
Method 1 Method 3 圖 4-26 模擬方法一與方法三反解 1040 nm 結構比較圖 兩種方法解出各結構半高寬，無明顯差異，結果如表表格二所示： 表格二 比較各結構在不同方法下 FWHM FWHM 1040 nm 180 nm 1000 nm 150, 180 nm 方法一 1060 nm 190 nm 1050 nm 180, 200 nm 方法三 1070 nm 190 nm 1070 nm 180, 200 nm 顯微鏡影像 1100 nm 420 nm 1090 nm 410, 420 nm50
第五章 結論與未來展望
本論文使用反摺積運算，配合空間截止頻率低通濾波器，成功以一個 已知樣品校正反解光學系統 PSF。並將本實驗複合理論方法求得之 PSF，帶 入 EM-MLE 後可成功運算出解析度小於繞射極限的線寬，達到突破繞射極 限的目標，此突破繞射極限的結果定性上與文獻相同。 另外深入比較以本實驗與文獻中解出 PSF 的方法，因為無法由 PSF 圖 型比較出差異，因此需透過討論 EM-MLE 反解樣品結果，經討論發現我們 提出方法造成影像失真之情形較不嚴重。透過模擬可確認此差異的確來自 於求得 PSF 方法的不同。另外在 EM-MLE 的過程中，在停止迭代時機為， 每次迭代較前次迭代變化小於整體結構 0.01%的條件下，迭代收斂次數可減 少 97 到 144 次，還原影像的速度提升。 最後在推導過程中修正 EM-MLE 迭代方程式，本實驗中使用修正迭代 方程式，在同樣停止迭代條件下，可以使迭代過程次數減少 8 到 21 次，同 樣可使還原影像速度提升，且還原之影像無差異。 在未來，我們期望可以製作三維樣品，將我們論文的方法推導至三維 空間，可建構出樣品的三維空間影像圖，並以光學顯微鏡標準校正片，定 量探討顯微鏡解析度、對比度、景深等參數。另外也可嘗試將空間頻率的 高頻項反解出，希望可反解出不失真的光學顯微鏡影像圖。最後希望可比 較此方法與目前普遍的盲蔽反摺積演算法(Blind deconvolution)有何異同。51
參考文獻
[1] J. Goodman (2ed.), Introduction to Fourier Optics, Roberts and Company, United States, (2005).
[2] J. B. Pawley, Handbook of Biological Confocal Microscopy, Springer, Berlin, (2006).
[3] D. A. Agard, Annu. Rev. Biophys. 13, 191 (1984).
[4] A. Erhardt, G. Zinser, D. Komitowski, and J. Bille, Appl. Opt. 24, 194 (1985).
[5] A. P. Dempster, N. M. Laird, and D. B. Rubin, J. R. Stat. Soc. A: Stat. 39, 1 (1977).
[6] L. A. Shepp, and Y. Vardi, IEEE Trans. Med. Imaging 1, 113 (1982). [7] T. J. Holmes, and N. J. O'connor, J. Microsc. 200, 114 (2000).
[8] J. G. McNally, C. Preza, J. Conchello, and L. J. Thomas, J. Opt. Soc. Am. A 11, 1056 (1994).
[9] Y. Garini, B. J. Vermolen, and I. T. Young, Curr. Opin. Biotechnol. 16, 3 (2005).
[10] N. S. Claxton, T. J. Fellers, and M. W. Davidson, Der Hautarzt Journal 61, 421 (2010).
[11] D. K. Hamilton, and T. Wilson, J. Phys. E: Sci. Instrum. 19, 52 (1986). [12] S. Cova, M. Ghioni, A. Lacaita, C. Samori, and F. Zappa, Appl. Opt. 35,
52
[13] I. K. Ilev, and R. W. Waynant, Rev. Sci. Instrum. 71, 4161 (2000).
[14] D. A. Agard, Y. Hiraoka, P. Shaw, and J. W. Sedat, Methods Cell Biol. 30, 353 (1989).
[15] D. W. Kammler (2ed.), A First Course in Fourier Analysis, Cambridge University, United Kingdom, (2008).
[16] Y. Vardi, L. A. Shepp, and L. Kaufman, J. Am. Stat. Assoc. 80, 8 (1985). [17] C. H. Lee, H. Y. Chiang, and H. Y. Mong, Opt. Lett. 28, 1772 (2003). [18] G. M. P. V. Kempen, H. T. M. V. D. Voort, J. G. J. Bauman, and K. C.
Strasters, IEEE Eng. Med. Biol. Mag. 15, 76 (1996). [19] J. A. Conchello, J. Opt. Soc. Am. A 15, 2609 (1998).
[20] J. G. McNally, T. Karpova, J. Cooper, and J. Conchello, Methods 19, 373 (1999).
[21] J. M. Liu (2ed.), Photonic Devices, Cambridge University Press, United States, (2005).
![圖 1-1 Blind deconvolution 反解前(a)及反解後(b)比較圖[7]](https://thumb-ap.123doks.com/thumbv2/9libinfo/8240935.171323/13.892.140.790.458.877/圖11Blinddeconvolution反解前a及反解後b比較圖7.webp)
