Phase-space Deconvolution for Light FIeld Microscopy算法解读

Phase-space Deconvolution for Light FIeld Microscopy算法解读

Code

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)$表示轴向z处的横截面上的三维坐标点的横向坐标。
• $x=(x_1,x_2)$表示传感器平面上微透镜的中心点坐标。
• $u=(u_1,u_2)$表示传感器平面上的任意像素相对微透镜中心的横向坐标，即横向偏移量。

注意：

• 图中，MLA上不同颜色的传感器像素具有不同的空间频率。
• 对比（b）和（c），相比于有周期性变化的空间PSF，右图中相空间PSF更加平滑连续

2.2.空间域PSF转化为相空间PSF

首先，数学表示传统的LF和相空间PSF
空间域PSF：$h_l(x+u,p,z)$。其中$x+u$决定了传感器上像素点的坐标，$(p,z)$决定了三维空间内体素的坐标。
相空间PSF：$h_p(x,u,p,z)$。其中$(x,u)$决定了传感器上像素点的坐标，$(p,z)$决定了三维空间内体素的坐标。

核心原理：由于每个微透镜在LFM原理图中都是局部傅里叶变换，我们可以直接对LFM捕获的像素进行重新排列，从而获得4D相空间的测量结果，而不需要按以往方法先标定图像。

从空间中轴向深度$z$的横截面上$(p_1,p_2)$处的体素发出的光经过LFM系统到达传感器平面上$x^{\prime }=(x_1^{\prime },x_2^{\prime })$处的过程可以用一个带目标约束的球面复合场表示$U_z\left({\mathbf{x}}^{\prime }-\mathbf{p}\right)$。因此，相空间PSF：$h_p(\mathbf{x},\mathbf{p},z,\mathbf{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​​∥∥​Kz​Fωx′⋅x​​{Uz​(x′−p)⋅t(x′−x)}⋅s(u)∥∥​22​dωx′−x​,
其中，$t,s$表示2D矩形空间采样窗。$F$表示傅立叶变换。

要在算法层面实现空域PSF到相空间PSF的转化，我们还需要对其成像过程进行离散化，每个像素的采样（即相空间PSF）可以表达为：
$H_{i,j,k}={\int }_{{\beta }_i}\left({\int }_{{\alpha }_j}{\int }_{{\gamma }_k}{\left|h_p(\mathbf{x},\mathbf{p},z,\mathbf{u})\right|}^{2}d\mathbf{u}d\mathbf{x}\right)d\mathbf{p}dz$,

其中，${\alpha }_j$表示在相空间$(j_1,j_2)$位置的空间像素的面积，这个坐标系是在微透镜中心$N_j=(N_{j1},N_{j2})$处建立的。$bet{a}_i$表示物体的体素体积，坐标位置为$i=(i_1,i_2,i_3)$。${\gamma }_k$表示传感器平面像素的面积，坐标位置是$k=(k_1,k_2)$，其对应于这个位置子孔径的特定空间频率范围。

算法实现层面上分为3步：

• 第一，固定深度z。以不同微透镜相对MLA中心的方向给不同视角的空域PSF施加一个偏移量，大小范围为$(-(Nnum-1)/2,(Nnum-1)/2)$。
• 第二，按相对微透镜中心的方向（一个相空间的状态）提取patch（大小$patch\_size=Nnum$）。以上图为例，视角数为$13\times 13$，$patch\_size=13$，patch的个数为15。
• 第三，不同视角下具有相同$(u_1,u_2)$——即具有相同的相对于微透镜中心的方向，此时空间频率相一致——的patch按倒序排列在一起的到相空间PSF。

实现代码见Code。

2.3 空间域$LF$转换为相空间$WDF$

主要分为三步：

1. 将空间域$LF$按照传感器平面坐标$x$和$u$进行重新排布，获得低分辨率的$WD{F}_L$。
   注意：这时我们发现由于微透镜阵列的物理限制（传感器只能捕捉微透镜区域的光线），无法对$x$完全采样，这表明相空间的采样较低。这正是伪影问题的来源。
2. 考虑到相空间的平滑约束，我们可以对$WD{F}_L(x,u)$进行三次插值，提高分辨率。
3. 应用子孔径滤波器（一个基于坐标$u$的低通滤波器）来消除插值的边缘膨胀。
   因此，用于下一步Deconvolution的高分辨率相空间$WDF$可计算得到：
   $$WDF(\mathbf{x},\mathbf{u})=cubic(WD{F}_L(\mathbf{x},\mathbf{u}))\otimes {\text{ filter }}_{\mathrm{u}}$$
   其中，$cubic(\ast )$代表2D三次插值，可以提升相空间的平滑约束。$filte{r}_u$表示一个关于子孔径频率$u$的低通滤波器。因此，成像模型可以转换为在相空间域的表示形式：
   $$WDF(\mathbf{x},\mathbf{u})={\iint }_{z,p}g(\mathbf{p},z){\left|h_p(\mathbf{x},\mathbf{p},z,\mathbf{u})\right|}^{2}dpdz$$

2.4 相空间Deconvolution实现3D重建

根据2.2.节公式2和2.3.节公式2，相空间$WDF(x,u)$和物体$g(p,z)$可以描述为：
$$\begin{array}{c}{Y}_{j,k}={\int }_{{\alpha }_j}{\int }_{{\gamma }_k}WDF(\mathbf{x},\mathbf{u})dudx,\\ X_1={\int }_{{\beta }_1}g(\mathbf{p},z)dpdz\end{array}$$
这时，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\times 13$）。散粒噪声是重建体中最主要的噪声，一般符合Poisson分布，因此成像模型可以进一步表示为：
$$\begin{array}{rl}& Y_{j,k}=Pois(\sum _i{X}_i\cdot H_{i,j,k})\end{array}$$
先验知识：由于相位空间中不同PSF的能量差异，不同的空间频率分量具有不同的散粒噪声分布。具体而言，高空间频率分量的PSF通常比低空间频率分量的PSF强度低，即散粒噪声方面的SNR较低。
为了利用这些先验知识，我们分别将重建问题作为不同空间频率分量的优化来解决。（这是受用于相干成像的Ptychographic方法的启发）在实践中，我们使用Richardson-Lucy算法按权重数$w_k$顺序更新具有不同空间频率分量的体积，以平衡不同空间频率分量的散粒噪声分布。在大量迭代$iter$中，我们通过以下迭代公式更新体积：
$$\begin{array}{rl}{X}_1^{(k,\text{ iter })}& \leftarrow \left(1-c{w}_k\right)\cdot X_1^{(k-1,iter)}\\ & +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,tter)}\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{array}$$
其中，$\odot$表示点积算子，$./$表示点除算子。$c$是一个常数，用于调整收敛速度（一般$c=80$）。$X_i^{(k,iter)}$表示第$iter$次迭代来自第k个空间频率分量的体积更新。$iter$表示最大迭代次数。
对应空间频率$u_k$的权重$w_k$可以基于PSF的能量分布计算出：
$$w_k=\frac{{\sum }_i{\Vert H_{i,k}\Vert }_1}{{\sum }_{k\in \Omega }{\sum }_i{\Vert H_{i,k}\Vert }_1}.$$
便利所有的空间频率$u_k$后，我们的算法会进入下一次迭代，即$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.