影像重構

由 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)：

[ ]

( ) ( )

P y =T δ x (3.4)

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 x

object x

psf x

訊雜比(Signal-to-noise ratio)大幅下降，將造成後續影像處理上困難，因此需 用另一方法得到 PSF。

在此論文中使用的方法為使用兩次(3.6)式，首先使用一個已知其形貌

15

的樣品，在此稱為校正樣品(Calibration sample)，經由光學顯微系統量測，

用數學式描述如(3.9)式。

Calibration Calibration

psf ⊗sample =image _(3.9)

將得到的影像(Calibration image)與校正樣品做反摺積運算可得 PSF，此方法 不需製作如數學定義般點狀樣品。

上式中，J₁為第一類貝索函數(Bessel functions of the first kind)，N 為數值孔 徑(NA, Numerical aperture)，λ₀為光源的波長，

r

為在像平面上的一點到光 軸的距離。將(3.10)式傅立葉轉換，可以得到光傳遞函數(OTF, Optical transfer function)：

0 1, 1/ 2 cut-off frequency)，其物理意義為，以頻譜空間討論繞射極限問題：

0

2

cut off

frequency N

λ

− = (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)：

17 續與週期性，兩空間週期皆為N ，或稱為窗口大小(Window size)，(3.13)、

(3.14)式可以簡化為只討論一個窗口，即稱為離散傅立葉轉換(DFT, Discrete Fourier transform)：

1

18

19

Time domain Frequency domain

( )a

add add

add

add add

add 同理，傅立葉逆轉換使用 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

1.寫出似然函數L(likelihood function)，其定義為當系統參數(Parameter) 為 p時量到Data的條件機率有多少：

21

3-4-2 最大期望演算法(EM)

EM 用迭代法(Iterative method)尋找系統參數，使似然函數最大值，如 前段假設系統有一個未知參數會影響到測量數據，詳細步驟如下：

1.將 p_k 當成系統參數，寫出似然函數L，如(3.28)式。

2.尋找一個新的參數p_k+1代入步驟 1 的結果大於將 pk_{代入步驟 1 的結}

果，表示系統在參數 pk₊1_{的情況下，比在參數為} p_{k的情況下更有可能} 量到已知Data：

( ) ( )

1 1

( _k ) | _k ( _k) | _k

L p ₊ =P Data p ₊ ≥ L p = P Data p (3.30) 3.除非 p_k 收斂，否則將 p_k+1取代 p_k回到步驟 1。

3-4-3EM-MLE

不難看出 MLE 困難之處，為如何表示系統的似然函數，與在系統的參 數較多，並且L對 p偏微分式較複雜的情況下，即使列出方程式，也未必 能求出未知參數。而 EM 則具有計算量較大，難寫出迭代方程式，與同樣 需寫出系統似然函數等問題。

為解決以上問題，Y. Vardi 等人將光學系統引入[16] ，再將 EM 與 MLE 兩種方法結合起來，首先假設光學系統對測量到的訊號，會有一雜訊分佈，

其分佈符合帕松分佈(Poisson distribution)，在 4-5 雜訊分析章節會再次驗證 此假設的合理性，帕松分佈函數如下：

22 將上述參數代入帕松分佈函數(3.31)，可得到在樣品為 object [x]的條件下，

系統的實際值：

image x

psf x x object x x

noise

psf x x object x

P e

image x

′

po noise

noise

23

image x

psf x x object x x

x

psf x x object x

L e

image x

′

通常 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 object x

image x psf x x object x

psf x x object x

′

∂ = − −

∂  − ′ ′ 

∑ ∑

^(3.39)

24

object 0 α

∂ ≤

∂ (3.43)

因此

α

與 object [x]關係為負相關(Negative correlation)，示意圖如圖 3-3

1

α

object Answer

object_Ans

圖 3-3

α

與 object 關係圖 考慮某次迭代：

k Ans

object >object _(3.44)

其對應

α

的量值將介於 0 到 1 之間，將object_k 乘以

α

可得到更接近

objectAns的object_k+1，如圖 3-4 左；同理可知，如果object_k 小於object_Ans， 其對應

α

的量值將大於 1，將object_k 乘以

α

也可得到更接近object_Ans的

25

objectk

objectk

α ×

1

α

object Answer

objectk

objectk

α ×

object x object x

k object x

object x object x image x

psf x x object x

26 L object

image x psf x x object x image x

′ (Stirling's approximation)：

ln( !) n ≈ n ln n n −

_(3.49) L object

image x psf x x object x image x image x image x

′

object x object x image x

psf x x object x

27 L object

psf x x object x image x

image x image x

′

image x image x

psf x x object x psf x x

image x

′

28

3-5 影像重構程式

Start

End Ture

False k←k+1

i = N-1?

object x object x image x

psf x x object x objectc

Solve object through EM-MLE Solve PSF through

Inverse filtering

1

IMAGE k

k k

OBJECT k object x IMAGE k image x

← ℑ

← ℑ

δ δ< 0

圖 3-5 影像重構程式流程圖

29

影像重構程式以 Labview^TM程式撰寫而成，分成以逆濾波器(Inverse filtering)反解 PSF，與以 EM-MLE 反解未知樣品兩個部分，程式流程圖如 圖 3-5，首先在程式開始後，先輸入校正樣品object x_c[ ]與校正樣品顯微鏡 影像image x_c[ ]，窗口大小皆為N，經過傅立葉轉換相除完成反摺積，再將 頻域空間將空間截止頻率k₀以上的頻率濾除，也就是令實虛部皆為零，將 此結果逆轉換後得到 PSF。

第二部分將第一部分求得的 PSF 代入 EM-MEL，也就是(3.46)式，由兩 個迴圈所構成，第一個為 For loop，確保實空間窗口內每一個x_i都會被迭代 一次，第二個為 While loop，計算(3.45)中每次疊造成object_k 的變化，如果 變化小於目標δ₀，代表繼續迭代下去，產生變化量已經不在討論的範圍內，

則程式停止；如變化小於或等於目標δ₀則繼續迭代。

30

在文檔中 利用反摺積求得點擴散函數進行光學顯微鏡影像還原術 (頁 23-40)