文章目录
Phase-space Deconvolution for Light FIeld Microscopy算法解读
1.摘要
本文记录了笔者研读Phase-space Deconvolution for Light FIeld Microscopy时,对算法原理和复现的个人理解。这篇文章创新性地利用光场相空间的平滑性的先验知识,通过在相空间域中对成像过程进行建模,将空间域中非均匀的点扩散函数PSF转换为更小尺寸的相空间域的均匀点扩散函数。这一创新有效地解决了传统的基于波动光学传播模型进行deconvolution重建时,由于系统误差(如光学像差和背景荧光)导致的周期性伪影,进而导致对比度降低的问题。
2.算法原理
2.1.LFM系统架构
如上图(a)所示,是光场显微镜(Light Field microscopy,LFM)的系统结构,这是一个LF1.0的系统,微透镜阵列MLA位于中继平面上。为了方便描述,我们先定义坐标系:
- p = ( p 1 , p 2 ) p=(p_1,p_2) p=(p1,p2)表示轴向z处的横截面上的三维坐标点的横向坐标。
- x = ( x 1 , x 2 ) x=(x_1,x_2) x=(x1,x2)表示传感器平面上微透镜的中心点坐标。
- u = ( u 1 , u 2 ) u=(u_1,u_2) u=(u1,u2)表示传感器平面上的任意像素相对微透镜中心的横向坐标,即横向偏移量。
注意:
- 图中,MLA上不同颜色的传感器像素具有不同的空间频率。
- 对比(b)和(c),相比于有周期性变化的空间PSF,右图中相空间PSF更加平滑连续
2.2.空间域PSF转化为相空间PSF
首先,数学表示传统的LF和相空间PSF
空间域PSF:
h
l
(
x
+
u
,
p
,
z
)
h_l(x+u,p,z)
hl(x+u,p,z)。其中
x
+
u
x+u
x+u决定了传感器上像素点的坐标,
(
p
,
z
)
(p,z)
(p,z)决定了三维空间内体素的坐标。
相空间PSF:
h
p
(
x
,
u
,
p
,
z
)
h_p(x,u,p,z)
hp(x,u,p,z)。其中
(
x
,
u
)
(x,u)
(x,u)决定了传感器上像素点的坐标,
(
p
,
z
)
(p,z)
(p,z)决定了三维空间内体素的坐标。
核心原理:由于每个微透镜在LFM原理图中都是局部傅里叶变换,我们可以直接对LFM捕获的像素进行重新排列,从而获得4D相空间的测量结果,而不需要按以往方法先标定图像。
从空间中轴向深度
z
z
z的横截面上
(
p
1
,
p
2
)
(p_1,p_2)
(p1,p2)处的体素发出的光经过LFM系统到达传感器平面上
x
′
=
(
x
1
′
,
x
2
′
)
x'=(x'_1,x'_2)
x′=(x1′,x2′)处的过程可以用一个带目标约束的球面复合场表示
U
z
(
x
′
−
p
)
U_{z}\left(\boldsymbol{x}^{\prime}-\boldsymbol{p}\right)
Uz(x′−p)。因此,相空间PSF:
h
p
(
x
,
p
,
z
,
u
)
h_{p}(\boldsymbol{x}, \boldsymbol{p}, z, \boldsymbol{u})
hp(x,p,z,u)可以数学表示为:
h
p
(
x
,
p
,
z
,
u
)
=
∫
ω
x
′
−
x
∥
K
z
F
ω
x
′
⋅
x
{
U
z
(
x
′
−
p
)
⋅
t
(
x
′
−
x
)
}
⋅
s
(
u
)
∥
2
2
d
ω
x
′
−
x
,
h_{p}(\boldsymbol{x}, \boldsymbol{p}, z, \boldsymbol{u})=\int_{\omega_{\mathbf{x}^{\prime}-\mathrm{x}}}\left\|K_{z} F_{\omega_{\boldsymbol{x}^{\prime} \cdot \boldsymbol{x}}}\left\{U_{z}\left(\boldsymbol{x}^{\prime}-\boldsymbol{p}\right) \cdot t\left(\mathrm{x}^{\prime}-\boldsymbol{x}\right)\right\} \cdot s(\boldsymbol{u})\right\|_{2}^{2} d \boldsymbol{\omega}_{\mathrm{x}^{\prime}-\mathrm{x}},
hp(x,p,z,u)=∫ωx′−x∥∥KzFωx′⋅x{Uz(x′−p)⋅t(x′−x)}⋅s(u)∥∥22dωx′−x,
其中,
t
,
s
t,s
t,s表示2D矩形空间采样窗。
F
F
F表示傅立叶变换。
要在算法层面实现空域PSF到相空间PSF的转化,我们还需要对其成像过程进行离散化,每个像素的采样(即相空间PSF)可以表达为:
H
i
,
j
,
k
=
∫
β
i
(
∫
α
j
∫
γ
k
∣
h
p
(
x
,
p
,
z
,
u
)
∣
2
d
u
d
x
)
d
p
d
z
H_{i, j, k}=\int_{\beta_{i}}\left(\int_{\alpha_{j}} \int_{\gamma_{k}}\left|h_{p}(\boldsymbol{x}, \boldsymbol{p}, z, \boldsymbol{u})\right|^{2} d \boldsymbol{u} d \boldsymbol{x}\right) d \boldsymbol{p} d z
Hi,j,k=∫βi(∫αj∫γk∣hp(x,p,z,u)∣2dudx)dpdz,
其中, α j \alpha_j αj表示在相空间 ( j 1 , j 2 ) (j_1,j_2) (j1,j2)位置的空间像素的面积,这个坐标系是在微透镜中心 N j = ( N j 1 , N j 2 ) N_j=(N_{j1},N_{j2}) Nj=(Nj1,Nj2)处建立的。 b e t a i beta_i betai表示物体的体素体积,坐标位置为 i = ( i 1 , i 2 , i 3 ) i=(i_1,i_2,i_3) i=(i1,i2,i3)。 γ k \gamma_k γk表示传感器平面像素的面积,坐标位置是 k = ( k 1 , k 2 ) k=(k_1,k_2) k=(k1,k2),其对应于这个位置子孔径的特定空间频率范围。
算法实现层面上分为3步:
- 第一,固定深度z。以不同微透镜相对MLA中心的方向给不同视角的空域PSF施加一个偏移量,大小范围为 ( − ( N n u m − 1 ) / 2 , ( N n u m − 1 ) / 2 ) (-(Nnum-1)/2,(Nnum-1)/2) (−(Nnum−1)/2,(Nnum−1)/2)。
- 第二,按相对微透镜中心的方向(一个相空间的状态)提取patch(大小 p a t c h _ s i z e = N n u m patch\_size=Nnum patch_size=Nnum)。以上图为例,视角数为 13 × 13 13\times13 13×13, p a t c h _ s i z e = 13 patch\_size=13 patch_size=13,patch的个数为15。
- 第三,不同视角下具有相同 ( u 1 , u 2 ) (u_1,u_2) (u1,u2)——即具有相同的相对于微透镜中心的方向,此时空间频率相一致——的patch按倒序排列在一起的到相空间PSF。
实现代码见Code。
2.3 空间域 L F LF LF转换为相空间 W D F WDF WDF
主要分为三步:
- 将空间域
L
F
LF
LF按照传感器平面坐标
x
x
x和
u
u
u进行重新排布,获得低分辨率的
W
D
F
L
WDF_L
WDFL。
注意:这时我们发现由于微透镜阵列的物理限制(传感器只能捕捉微透镜区域的光线),无法对 x x x完全采样,这表明相空间的采样较低。这正是伪影问题的来源。 - 考虑到相空间的平滑约束,我们可以对 W D F L ( x , u ) WDF_L(x,u) WDFL(x,u)进行三次插值,提高分辨率。
- 应用子孔径滤波器(一个基于坐标
u
u
u的低通滤波器)来消除插值的边缘膨胀。
因此,用于下一步Deconvolution的高分辨率相空间 W D F WDF WDF可计算得到:
W D F ( x , u ) = c u b i c ( W D F L ( x , u ) ) ⊗ filter u \left.W D F(\boldsymbol{x}, \boldsymbol{u})=cubic(WDF_L(\boldsymbol{x}, \boldsymbol{u})\right) \otimes \text { filter }_{\mathrm{u}} WDF(x,u)=cubic(WDFL(x,u))⊗ filter u
其中, c u b i c ( ∗ ) cubic(*) cubic(∗)代表2D三次插值,可以提升相空间的平滑约束。 f i l t e r u filter_u filteru表示一个关于子孔径频率 u u u的低通滤波器。因此,成像模型可以转换为在相空间域的表示形式:
W D F ( x , u ) = ∬ z , p g ( p , z ) ∣ h p ( x , p , z , u ) ∣ 2 d p d z W D F(\boldsymbol{x}, \boldsymbol{u})=\iint_{z,p} g(\boldsymbol{p}, z)\left|h_{p}(\boldsymbol{x}, \boldsymbol{p}, z, \boldsymbol{u})\right|^{2} d p d z WDF(x,u)=∬z,pg(p,z)∣hp(x,p,z,u)∣2dpdz
2.4 相空间Deconvolution实现3D重建
根据2.2.节公式2和2.3.节公式2,相空间
W
D
F
(
x
,
u
)
WDF(x,u)
WDF(x,u)和物体
g
(
p
,
z
)
g(p,z)
g(p,z)可以描述为:
Y
j
,
k
=
∫
α
j
∫
γ
k
W
D
F
(
x
,
u
)
d
u
d
x
,
X
1
=
∫
β
1
g
(
p
,
z
)
d
p
d
z
\begin{gathered} Y_{j, k}=\int_{\alpha_{j}} \int_{\gamma_{k}} W D F(\boldsymbol{x}, \boldsymbol{u}) d u d x, \\ X_{1}=\int_{\beta_{1}} g(\boldsymbol{p}, z) d p d z \end{gathered}
Yj,k=∫αj∫γkWDF(x,u)dudx,X1=∫β1g(p,z)dpdz
这时,2.3.节公式2的成像模型可以离散化为:
Y
j
,
k
=
∑
i
X
i
⋅
H
i
,
j
,
k
,
k
=
{
Ω
∣
(
k
1
,
k
2
)
∈
Ω
}
,
\begin{aligned} &Y_{j, k}=\sum_{i} X_{i} \cdot H_{i, j, k}, \\ &k=\left\{\Omega \mid\left(k_{1}, k_{2}\right) \in \Omega\right\}, \end{aligned}
Yj,k=i∑Xi⋅Hi,j,k,k={Ω∣(k1,k2)∈Ω},
其中,
Ω
\Omega
Ω是所有空间频率分量的集合,数量是每个微透镜后的像素总数(比如
13
×
13
13\times13
13×13)。散粒噪声是重建体中最主要的噪声,一般符合Poisson分布,因此成像模型可以进一步表示为:
Y
j
,
k
=
P
o
i
s
(
∑
i
X
i
⋅
H
i
,
j
,
k
)
\begin{aligned} &Y_{j, k}=Pois(\sum_{i} X_{i} \cdot H_{i, j, k}) \end{aligned}
Yj,k=Pois(i∑Xi⋅Hi,j,k)
先验知识:由于相位空间中不同PSF的能量差异,不同的空间频率分量具有不同的散粒噪声分布。具体而言,高空间频率分量的PSF通常比低空间频率分量的PSF强度低,即散粒噪声方面的SNR较低。
为了利用这些先验知识,我们分别将重建问题作为不同空间频率分量的优化来解决。(这是受用于相干成像的Ptychographic方法的启发)在实践中,我们使用Richardson-Lucy算法按权重数
w
k
w_k
wk顺序更新具有不同空间频率分量的体积,以平衡不同空间频率分量的散粒噪声分布。在大量迭代
i
t
e
r
iter
iter中,我们通过以下迭代公式更新体积:
X
1
(
k
,
iter
)
←
(
1
−
c
w
k
)
⋅
X
1
(
k
−
1
,
i
t
e
r
)
+
c
w
k
⋅
X
i
(
k
−
1
,
t
ter
)
⊙
(
∑
j
(
Y
~
j
,
k
⋅
/
(
∑
i
X
i
(
k
−
1
,
t
t
e
r
)
⋅
H
i
,
j
,
k
)
)
⋅
H
i
,
N
j
−
j
,
k
)
⋅
/
(
∑
j
1
⋅
H
i
,
N
j
−
j
,
j
,
k
)
\begin{aligned} X_{1}^{(k, \text { iter })} & \leftarrow\left(1-c w_{k}\right) \cdot X_{1}^{(k-1, i t e r)} \\ &+c w_{k} \cdot X_{i}^{(k-1, t \text { ter })} \odot\left(\sum_{j}\left(\tilde{Y}_{j, k} \cdot /\left(\sum_{i} X_{i}^{(k-1, t t e r)} \cdot H_{i, j, k}\right)\right) \cdot H_{i, N_{j}-j, k}\right) \cdot /\left(\sum_{j} 1 \cdot H_{i, N_{j}-j, j, k}\right) \end{aligned}
X1(k, iter )←(1−cwk)⋅X1(k−1,iter)+cwk⋅Xi(k−1,t ter )⊙(j∑(Y~j,k⋅/(i∑Xi(k−1,tter)⋅Hi,j,k))⋅Hi,Nj−j,k)⋅/(j∑1⋅Hi,Nj−j,j,k)
其中,
⊙
\odot
⊙表示点积算子,
.
/
./
./表示点除算子。
c
c
c是一个常数,用于调整收敛速度(一般
c
=
80
c=80
c=80)。
X
i
(
k
,
i
t
e
r
)
X_i^{(k,iter)}
Xi(k,iter)表示第
i
t
e
r
iter
iter次迭代来自第k个空间频率分量的体积更新。
i
t
e
r
iter
iter表示最大迭代次数。
对应空间频率
u
k
u_k
uk的权重
w
k
w_k
wk可以基于PSF的能量分布计算出:
w
k
=
∑
i
∥
H
i
,
k
∥
1
∑
k
∈
Ω
∑
i
∥
H
i
,
k
∥
1
.
w_{k}=\frac{\sum_{i}\left\|H_{i, k}\right\|_{1}}{\sum_{k \in \Omega} \sum_{i}\left\|H_{i, k}\right\|_{1}} .
wk=∑k∈Ω∑i∥Hi,k∥1∑i∥Hi,k∥1.
便利所有的空间频率
u
k
u_k
uk后,我们的算法会进入下一次迭代,即
i
t
e
r
+
1
iter+1
iter+1,直到算法收敛。整个算法流程如下图所示。
3.笔者按
这个方向同道甚寡,希望有同志一起学习交流鸭~
参考文献
[1] Lu Z, Wu J, Qiao H, et al. Phase-space deconvolution for light field microscopy[J]. Optics express, 2019, 27(13): 18131-18145.