zehua

菜单

Administrator
发布于 2026-09-09 / 0 阅读
0
0

信号处理与成像系统 TP1（续）：PALM 二维定位算法

#实验报告#M2 ATSI

引言

在TP1中，我们详细推导了PALM（Photo-Activated Localization Microscopy）技术的一维定位算法，包括信号模型的建立、最大似然（ML）估计器的构造以及Cramér-Rao下界（CRLB）的计算。这些方法成功实现了亚像素精度的荧光分子定位。然而，实际的PALM成像是在二维空间中进行的，一个荧光分子的位置由两个坐标参数 (\tau, \eta) 描述，分别对应图像的行和列方向。本额外TP的目标是将一维ML估计方法扩展到二维情况，构造完整的二维ML估计器，并通过实验验证其性能。

二维观测模型

信号模型的建立

考虑一个荧光分子位于位置 (\tau, \eta)，其中 \tau 是行方向（x 方向）的位置，\eta 是列方向（y 方向）的位置。当成像系统采集图像时，第 (i, j) 个像素（i 是行索引，j 是列索引）观测到的信号为：

s_{ij} = a r_{ij}(\tau, \eta) + b_{ij}
  • a 是信号幅度，反映了荧光分子的发光强度
  • r_{ij}(\tau, \eta) 是第 (i, j) 个像素上的二维离散化点扩散函数（PSF）
  • b_{ij} 是加性高斯白噪声，满足 b_{ij} \sim \mathcal{N}(0, \sigma^2)，且各像素之间的噪声相互独立。

二维PSF的可分离性

由于PSF是各向同性（isotropic）高斯，这意味着二维PSF可以表示为两个一维PSF的乘积：

r_{ij}(\tau, \eta) = r_i(\tau) \cdot r_j(\eta)

其中 r_i(\tau) 是行方向的一维离散化PSF，r_j(\eta) 是列方向的一维离散化PSF。每个一维PSF都可以通过误差函数解析计算。具体地，对于行方向：

r_i(\tau) = \frac{1}{2} \left[ \text{erf}\left(\frac{(i+1)Dx - \tau}{\sqrt{2}\sigma}\right) - \text{erf}\left(\frac{iDx - \tau}{\sqrt{2}\sigma}\right) \right]

其中 Dx 是像素尺寸（离散化步长），\sigma 是PSF的高斯标准差，与半峰全宽（FWHM）w 的关系为 \sigma = wDx/(2\sqrt{2\ln 2})。列方向的PSF r_j(\eta) 使用相同的公式，只需将 \tau 替换为 \eta，i 替换为 j。

这个可分离性质在计算上非常有利。虽然我们需要构造一个 M \times N 的二维PSF矩阵（M 和 N 分别是ROI的行数和列数），但实际上只需要计算 M 个 r_i(\tau) 值和 N 个 r_j(\eta) 值，然后通过外积 r_i(\tau) \times r_j(\eta) 得到所有 r_{ij}(\tau, \eta)。这将计算复杂度从 O(MN) 降低到 O(M+N)（在误差函数计算方面）。

二维最大似然估计器

对数似然函数

在加性高斯白噪声假设下，观测数据 \{s_{ij}\} 的对数似然函数为：

\ell(\tau, \eta, a) = -\frac{1}{2\sigma^2} \sum_{i=1}^{M} \sum_{j=1}^{N} [s_{ij} - a r_{ij}(\tau, \eta)]^2 + \text{常数}

这个表达式涉及三个未知参数：两个位置参数 (\tau, \eta) 和一个幅度参数 a。

浓缩似然函数

直接优化三个参数会使问题复杂化。我们可以利用浓缩似然函数（concentrated likelihood）的技巧来简化问题，其基本思想是，对于给定的位置 (\tau, \eta)，通过解析求解得到最优幅度 a 。

对似然函数关于 a 求导并令其为零：

\frac{\partial \ell}{\partial a} = -\frac{1}{\sigma^2} \sum_{i=1}^{M} \sum_{j=1}^{N} r_{ij}(\tau, \eta) [s_{ij} - a r_{ij}(\tau, \eta)] = 0

解这个线性方程得到：

\hat{a}(\tau, \eta) = \frac{\sum_{i=1}^{M} \sum_{j=1}^{N} s_{ij} r_{ij}(\tau, \eta)}{\sum_{i=1}^{M} \sum_{j=1}^{N} r_{ij}^2(\tau, \eta)}

这个最优幅度是观测数据与PSF的内积除以PSF的范数平方，其集合意义是观测向量在PSF方向上的投影长度。

将最优幅度 \hat{a}(\tau, \eta) 代回原始似然函数，得到只依赖于位置参数的浓缩代价函数：

C(\tau, \eta) = \sum_{i=1}^{M} \sum_{j=1}^{N} [s_{ij} - \hat{a}(\tau, \eta) r_{ij}(\tau, \eta)]^2

二维ML估计器定义为使浓缩代价函数最小化的位置：

(\hat{\tau}_{ML}, \hat{\eta}_{ML}) = \arg\min_{(\tau, \eta)} C(\tau, \eta)

这是一个二维无约束非线性优化问题。联合优化两个参数 (\tau, \eta) 的好处是能够充分利用二维观测数据中的所有信息，包括两个方向之间存在的相关性。

从理论角度看，如果PSF不是完美的圆对称，或者存在某些系统性的像差使得 x 和 y 方向耦合，那么Fisher信息矩阵会包含非对角元素，联合优化能够利用这些耦合信息，从而达到理论上的最优性能（即达到CRLB）。

算法代价函数计算

在算法中，对于每个候选位置 (\tau, \eta)，代价函数 C(\tau, \eta) 的计算利用了二维高斯PSF的可分离性质，即分别计算行方向和列方向的一维离散化PSF。

二维PSF矩阵通过外积运算构造：

r_{ij}(\tau, \eta) = r_i(\tau) \times r_j(\eta)

随后，最优幅度估计值 \hat{a}(\tau, \eta) 通过观测数据与PSF的归一化内积计算得到。最终的代价函数值为残差平方和：

C(\tau, \eta) = \sum_{i=1}^{M} \sum_{j=1}^{N} [s_{ij} - \hat{a}(\tau, \eta) r_{ij}(\tau, \eta)]^2

实验结果

实验使用的PALM数据包括

  • 一个测试图像（ImageTest，尺寸70×100像素，包含17个已知位置的荧光分子）
  • PALM图像序列（ImagesPALM，尺寸70×100×593，共593帧）
  • 真实坐标数据（i_molecules和j_molecules）用于验证算法精度。

测试图像上的定位精度

首先在测试图像上验证算法，候选点检测阶段找到17个候选点，与真实数量一致，表明没有遗漏。对每个候选点应用二维ML估计器后，得到17个亚像素位置估计。

为了评估定位精度，我们将每个真实位置与最近的估计位置进行匹配，计算所有匹配对的定位误差（欧氏距离），然后求均方根误差（RMSE）：

\text{RMSE} = \sqrt{\frac{1}{N_{\text{matched}}} \sum_{k=1}^{N_{\text{matched}}} [(\hat{\tau}_k - \tau_k^{\text{true}})^2 + (\hat{\eta}_k - \eta_k^{\text{true}})^2]}

实验结果为RMSE = 0.1572 pixels

image-20251016153049027

从图像可以观察到，算法成功检测并定位了所有17个荧光分子，且精准的找到了荧光分子的真实坐标，且我们可见荧光分子的分布较为分散，没有PSF的重叠。每个分子周围的灰度分布呈现典型的高斯PSF特征，中心亮度最高，向外逐渐衰减，这对应了我们之前提到的荧光分子寻找方式（亮度最大值）

超分辨率图像重建

在完成所有593帧PALM序列的荧光分子定位后，获得了10,081个亚像素精度的分子位置坐标，我们需要将其叠加重建为超分辨率图像。

关于超分辨率重建的放大倍数，我们设定为8倍，也就是说重建图像的尺寸从原始的 70 \times 100 像素扩展到 560 \times 800 像素，原始图像中的每个像素被细分为 8 \times 8 = 64 个超分辨率像素，提高了空间采样密度。然后我们需要将定位得到的荧光分子亚像素坐标 (\hat{\tau}, \hat{\eta}) 需要映射到高分辨率网格上。

最后结果如下

image-20251016153418353

信号处理与成像系统 TP1：PAL...
信号处理与成像系统 TP2：环形二...

评论