激光扫描共聚焦显微镜是获得高分辨率荧光细胞图像和生物医学三维重建的最先进技术之一. 随着成像技术的进步, 针对微观细胞的研究也逐渐在深入. 这些研究需要涉及细胞的多维检测和分析、从二维细胞成像向三维细胞成像的演变、以及从传统空间成像到时间成像的扩展. 因此, 激光扫描共聚焦显微镜^[1]的研究对于生物医学及相关学科的发展具有重要意义. 三维可视化的关键技术在于提取一系列连续的二维切片, 通过激光扫描共聚焦显微镜重建厚生物组织的三维模型, 然后进行定量分析. 因此, 共聚焦图像复原方法对于活细胞的检测和三维结构的探索具有很高的价值.

由于光学成像系统中存在点扩散, 这会导致卷积效应和原始图像的退化, 因此难以对共聚焦弱荧光图像执行三维操作. 对于较厚的组织, 随着成像深度的增加, 共聚焦图像的退化会加剧. 因此, 有必要通过对图像退化模型执行去卷积复原来消除退化的干扰. 图像模糊可以表示为清晰图像和点扩散函数(Point Spread Function, PSF)的卷积加噪声, 图像复原则是通过去卷积方法估计清晰的原始图像的操作.

在图像去噪声、去模糊这一类问题中, 最常用到的是全变分(Total Variation, TV)正则项^[2], 这是因为TV模型具有较为良好的边缘保持特性^[3]. 但是, TV正则化方法适合倾向于分段常数解, 即它对细节和纹理的复原效果不佳, 相当容易形成阶梯效应.

截至目前, 关于模糊图像去卷积方法的研究大部分都假定噪声服从于高斯分布, 采取此类模型的优势是: 在所得的待优化问题中, 会有二次项的引入, 这将便于求解. 但在医学或全息^[4]成像、天文观测等特定应用领域中, 噪声往往是服从泊松分布的, 此时若再采取高斯噪声模型, 概率分布特征将不能得到有效的描述.

泊松图像反卷积过程是具有病态特征的, 特别在低信噪比的情况下, 经数次迭代后会产生高的估计噪声. 针对泊松噪声污染下的模糊的去卷积问题, 目前最常用的方法是Richardson–Lucy方法^[5], 它是期望最大化算法的一个特例. 但这种方法具有明显的阶梯效应. 董文德等采用最优化方法将自然图像进行训练, 得到专家场模型的所有滤波器, 使其具有更真实的概率分布特征, 以此实现对泊松噪声污染模糊图像的有效复原^[6]. Zhang等提出了一种基于小波框架的正则化模型, 以去除泊松噪声而得到去模糊图像^[7]. Liu等在基于高阶全变分的泊松图像去卷积中引入了空间自适应正则化参数^[8].

为提高所得共聚焦图像的复原精度和图像采集速度, 本文将采用一种基于Hessian矩阵范数的图像复原正则化方法.

1 基于Hessian矩阵范数正则化的图像复原模型

共聚焦图像的噪声主要来自两个来源: 信号本身和成像系统颤振. 泊松噪声污染下的图像退化模型可表示为:

${{y}} = P\left( {{{Af}}} \right)$

(1)

其中, ${{y}}$ 表示观测图像, ${{f}}$ 表示原始图像, ${{A}}$ 是线性模糊算子, $P\left( {{{Af}}} \right)$ 表示含参数的独立同分布泊松随机变量.

首先将传统的全变分梯度方法从一维向二维进行扩展, 定义了用于计算二阶导数的Hessian算子为:

${{Hf}} = \left[ {\begin{array}{*{20}{c}} {{{{f}}_{xx}}}&{{{{f}}_{xy}}}\\ {{{{f}}_{yx}}}&{{{{f}}_{yy}}} \end{array}} \right]$

(2)

基于Hessian矩阵范数正则化的图像复原模型表示为:

$\widehat {{f}} = \mathop {\arg \min }\limits_{{f}} \left\{ {\left[ {{{Af}} - {{y}}\ln \left( {{{Af}}} \right)} \right] + \mu {{\left\| {{{Hf}}} \right\|}_2}} \right\}$

(3)

所用的是F范数, 其具有凸性、尺度不变性、平移不变性和旋转不变性, 因此, 可作为推广的二阶TV.

2 图像复原模型的求解

为解耦复杂的目标函数, 采取交替方向算法^[9], 引入辅助变量 ${{{u}}_1}$ , ${{{u}}_2}$ 作为变量 ${{Af}}$ , ${{f}}$ 的近似, 实现将变量 ${{f}}$ 分离的目的, 当然, 同时也要补充一个二次项, 用于惩罚 ${{Af}}$ , ${{f}}$ 和 ${{{u}}_1}$ , ${{{u}}_2}$ 之间的差异. 所得解耦模型:

$\left( {\widehat {{f}},\widehat {{{{u}}_1}},\widehat {{{{u}}_2}}} \right) = \mathop {\arg \min }\limits_{{{f}},{{{u}}_1},{{{u}}_2}} \left\{ \begin{array}{l} \left[ {{{{u}}_1} - {{y}}ln\left( {{{{u}}_1}} \right)} \right] + \mu {\left\| {{{H}}{{{u}}_2}} \right\|_2}\\ + \dfrac{\alpha }{2}\left[ \begin{array}{l} \left\| {{{{u}}_1} - \left( {{{Af}} + {{{b}}_1}} \right)} \right\|_2^2\\ + \left\| {{{{u}}_2} - \left( {{{f}} + {{{b}}_2}} \right)} \right\|_2^2 \end{array} \right] \end{array} \right\}$

(4)

其中, $\alpha > 0$ 表示惩罚系数(拉格朗日参数), ${{{b}}_1}$ , ${{{b}}_2}$ 为对偶变量.

模型所涉及的以上变量之间的关系如图1所示.

可将上述模型的求解分成如下3个步骤:

(1)在固定 ${{u}}_1^{(k - 1)}$ , ${{u}}_2^{(k - 1)}$ , ${{b}}_1^{(k - 1)}$ , ${{b}}_2^{(k - 1)}$ 的前提下, 求解 ${{{f}}^{(k)}}$ .

${{{f}}^{(k)}} = \mathop {\arg \min }\limits_{{f}} \frac{\alpha }{2}\left[ \begin{array}{l} \left\| {{{u}}_1^{(k - 1)} - \left( {{{Af}} + {{b}}_1^{(k - 1)}} \right)} \right\|_2^2\\ + \left\| {{{u}}_2^{(k - 1)} - \left( {{{f}} + {{b}}_2^{(k - 1)}} \right)} \right\|_2^2 \end{array} \right]$

(5)

(2)对于步骤(1)所得的 ${{{f}}^{(k)}}$ , 可用来求解 ${{u}}_1^{(k)}$ , ${{u}}_2^{(k)}$ .

${{u}}_1^{(k)} = \mathop {\arg \min }\limits_{{{{u}}_1}} \left\{ \begin{array}{l} \left[ {{{{u}}_1} - {{y}}\ln\left( {{{{u}}_1}} \right)} \right]\\ + \dfrac{\alpha }{2}\left\| {{{{u}}_1} - \left( {{{Af}^{(k)}} + {{b}}_1^{(k - 1)}} \right)} \right\|_2^2 \end{array} \right\}$

(6)

${{u}}_2^{(k)} = \mathop {\arg \min }\limits_{{{{u}}_2}} \left\{ {\mu {{\left\| {{{Hu}_2}} \right\|}_2} + \frac{\alpha }{2}\left\| {{{{u}}_2} - \left( {{{{f}}^{(k)}} + {{b}}_2^{(k - 1)}} \right)} \right\|_2^2} \right\}$

(7)

(3)更新 ${{b}}_1^{(k)}$ , ${{b}}_2^{(k)}$ .

${{b}}_1^{(k)} = {{b}}_1^{(k - 1)} + {{Af}^{(k)}} - {{u}}_1^{(k)}$

(8)

$ {{b}}_2^{(k)} = {{b}}_2^{(k - 1)} + {{{f}}^{(k)}} - {{u}}_2^{(k)} $

(9)

2.1 ${{{f}}^{(k)}}$ 的求解

对式(5)使用快速傅里叶变换, 在频域中按照式(10)进行 ${{{f}}^{(k)}}$ 的求解.

${{{f}}^{(k)}} = {F^{ - 1}}\left\{ \begin{array}{l} \left[ \begin{array}{l} F{\left( {{A}} \right)^ * } \cdot \left( {F\left( {{{u}}_1^{\left( {k - 1} \right)}} \right) - F\left( {{{b}}_1^{\left( {k - 1} \right)}} \right)} \right)\\ + \left( {F\left( {{{u}}_2^{\left( {k - 1} \right)}} \right) - F\left( {{{b}}_2^{\left( {k - 1} \right)}} \right)} \right) \end{array} \right]\\ /\left[ {F{{\left( {{A}} \right)}^ * } \cdot F\left( {{A}} \right) + {{I}}} \right] \end{array} \right\}$

(10)

其中, ${\left( \cdot \right)^ * }$ 表示共轭运算, $F\left( \cdot \right)$ 表示傅里叶变换, ${F^{ - 1}}\left( \cdot \right)$ 表示傅里叶逆变换.

2.2 ${{u}}_1^{(k)}$ 的求解

对于式(6), 考虑代价函数:

$C\left( x \right) = x - {{y}}\ln \left( x \right) + \frac{\alpha }{2}{\left[ {x - \left( {{{Af}^{\left( k \right)}} + {{b}}_1^{\left( {k - 1} \right)}} \right)} \right]^2}$

(11)

然后依次计算它的一阶导数:

$\frac{{{\rm{d}}C\left( x \right)}}{{{\rm{d}}x}} = 1 - \frac{{{y}}}{x} + \alpha x - \alpha \left( {{{Af}^{\left( k \right)}} + {{b}}_1^{\left( {k - 1} \right)}} \right)$

(12)

二阶导数:

$\frac{{{{\rm{d}}^2}C\left( x \right)}}{{{\rm{d}}{x^2}}} = \frac{{{y}}}{{{x^2}}} + \alpha $

(13)

这是个凸优化问题, 由式(12)为零得:

$\alpha {x^2} + \left[ {1 - \alpha \left( {{{Af}^{\left( k \right)}} + {{b}}_1^{\left( {k - 1} \right)}} \right)} \right]x - {{y}} = 0$

(14)

求出其非负根. 又由于向量的最优解可在其每个分量都取最优解时得到, 通过组合若干个相互独立的最优化问题式(11), 可知式(6)的解为:

$ {{u}}_1^{(k)} \!= \!\frac{{\alpha \left( {{{Af}^{\left( k \right)}} \!+ \!{{b}}_1^{\left( {k - 1} \right)}} \right) \!-\! 1 \!+\! \sqrt {{{\left[ {1 \!-\! \alpha \left( {{{Af}^{\left( k \right)}} \!+\! {{b}}_1^{\left( {k - 1} \right)}} \right)} \right]}^2} \!+ \!4\alpha {{y}}} }}{{2\alpha }} $

(15)

2.3 ${{u}}_2^{(k)}$ 的求解

对于式(7)中的Hessian矩阵范数正则化项, 首先可将其等价表示为:

${\left\| {{{Hu}_2}} \right\|_2} \!=\! \mathop {\max }\limits_{{\omega } \in {B_{\infty ,2}}} {\left\langle {{{\omega }},{{Hu}_2}} \right\rangle _{{R^{N \times 2 \times 2}}}} \!=\! \mathop {\max }\limits_{{\omega } \in {B_{\infty ,2}}} {\left\langle {{{{H}}^ * }{{\omega }},{{{u}}_2}} \right\rangle _{{R^N}}}\!\!\!\!$

(16)

其中, ${B_{\infty ,2}} = \{ \left. {{{\omega }} = {{\left[ {{{\omega }}_1^{\rm{T}},{{\omega }}_2^{\rm{T}}, \cdots ,{{\omega }}_N^{\rm{T}}} \right]}^{\rm{T}}} \in {R^{N \times 2 \times 2}}} \right|\left\| {{{{\omega }}_i}} \right\| \le 1, $ $ \forall i \!=\! 1,2, \cdots ,N \}$ 表示单位范数球, H*为H的共轭算子.

式(7)即改写为:

${{u}}_2^{(k)} = \mathop {\arg \min }\limits_{u \in T} \left\{ \begin{array}{l} \tau \mathop {\max }\limits_{\omega \in {B_{\infty ,2}}} {\left\langle {{H^ * }{\omega },{u_2}} \right\rangle _{{R^N}}}\\ + \dfrac{1}{2}\left\| {{u_2} - \left( {{f^{(k)}} + {{b}}_2^{(k - 1)}} \right)} \right\|_2^2 \end{array} \right\}$

(17)

其中, T是用于约束解的范围的凸集, $\tau = \mu /\alpha $ .

使用梯度投影方法^[10]可得:

${{u}}_2^{(k)} = {P_T}\left[ {\left( {{{{f}}^{(k)}} + {{b}}_2^{(k - 1)}} \right) - \tau {{{H}}^ * }{{\omega }}} \right]$

(18)

其中, ${P_T}$ 是凸集T上的正交投影算子.

投影到单位范数球 ${B_{\infty ,2}}$ 上的算子可以定义为:

${P_{{B_{\infty ,2}}}}\left( {{{{\omega }}_i}} \right) = \left\{ {\begin{aligned} &{{{{\omega }}_i},\;\;\left\| {{{{\omega }}_i}} \right\| \le 1}\\ &{\frac{{{{{\omega }}_i}}}{{\left\| {{{{\omega }}_i}} \right\|}},\;\;\left\| {{{{\omega }}_i}} \right\| > 1} \end{aligned}} \right.$

(19)

关于 ${{\omega }} $ 的迭代过程如下:

${{v}}^{(1)} = \omega ^{(1)} = 0$

while i < Maxiter(最大迭代次数) do

$\;\;\;\;\omega ^{(i)} \leftarrow {P_{{B_{\infty ,2}}}}{{v}}^{(i)} + \dfrac{1}{{64\tau }}{P_T}(({{f}}^{(k)} + {{b}}_2^{(k - 1)}) - \tau {{{H}}^*}{{v}}^{(i)})$

${t_{i + 1}} \leftarrow \dfrac{{1 + \sqrt {1 + 4t_i^2} }}{2}$

$\;\;\;\;{{v}}^{(i + 1)} \leftarrow \omega ^{(i)} + \dfrac{{{t_i} - 1}}{{{t_{i + 1}}}}(\omega ^{(i)} - \omega ^{(i - 1)})$

end do

${{u}}_2^{(k)} = {P_T}(({{f}}^{(k)} + {{b}}_2^{(k - 1)}) - \tau {{{H}}^*}{{\omega }}^{(i)})$

3 实验结果与分析

为了验证所提方法的有效性, 评估了通过不同惩罚项所获得的图像复原性能, 特别是, 对采用TV模型、高阶混合惩罚项^[8]、Hessian惩罚项的正则化进行比较. 考虑所提方法中的正则化参数 $\mu $ , 它起到平衡图像轮廓和细节的作用. 通常来说, 增加其值往往会获得更清晰的结果, 但也会放大噪音. 在这种情况下, 我们采用有效的策略, 通过实验性地调节参数值, 以获得最佳的视觉效果.

3.1 仿真验证

为了说明提出的方法在共聚焦图像上去模糊的复原效果, 本文首先在仿真数据上进行验证. 仿真运算的配置为Matlab r2014b, 2.60 GHz Intel Core i5-3230M CPU以及4 GB内存. 对于合成的共聚焦数据集, 参照了近期发布的软件confocalGN^[11], 模糊函数采用为方差为1的高斯点扩散函数. 由于泊松噪声是依赖于像素的灰度值的, 所以噪声水平可以由图形的最大灰度值来决定. 泊松噪声可以通过matlab的内置函数“poissrnd”来实现. 对于仿真数据, 采用2幅较为代表性的图片(“woman”以及“cameraman”)进行观察. 去模糊结果如图2、图3和图4所示, 其中, 图3是图2(c)、图2(d)、图2(e)对应于图2(a)中红框区域的放大显示.

此外, 为了对仿真结果作定量分析以估计算法的重建质量, 本文采用峰值信噪比(PSNR)以及均方误差(MSE)来度量, 图2、图4中(c)、(d)、(e)的指标对比如表1所示. 显然, 在两张经典图片上, 相比于TV模型、高阶混合惩罚项方法这两种参照, 本文所提方法复原结果的PSNR更高, MSE更低, 充分体现了其有效性.

3.2 实验设置

建立的实验平台^[12]如图5所示, 它是用于荧光图像重建的激光扫描共聚焦显微镜, 由激光光源、扫描装置、共聚焦装置和计算机控制的图像处理系统组成. 该平台被用于比较所提出的Hessian正则化模型和TV模型、高阶混合惩罚项方法.

本文选择用平台所观察的扫描断层图像作为原始图像, 采集的图像一般会被泊松噪声和散焦模糊所污染, 针对散焦模糊中的点扩散函数, 我们采用荧光小球来进行标定. 其具体步骤为: (1)通过实验获得不同尺寸的荧光微珠图像; (2)采用高斯点扩散函数近似法^[13]拟合荧光珠的分布; (3)估算并取点扩散函数的平均值. 我们注意到, 平均预处理不仅可以降低测量误差, 而且可以减轻噪声对点扩散函数标定的影响.

3.3 图像复原结果

基于标定的点扩散函数, 分别使用TV模型、高阶混合惩罚项方法和所提出的Hessian正则化模型来处理它, 并比较图像复原结果, 如图6所示. 其中, 图6(a)与图6(e)是不同视野内所采集到的样本图像, 图6(b)、图6(c)、图6(d)和图6(f)、图6(g)、图6(h)分别表示使用TV模型、高阶混合惩罚项方法和所提出的Hessian模型对样本1和2的去模糊结果.

从图6明显可以看出, 本文所提出的方法对模糊共聚焦图像进行了有效的复原, 图像质量有了显著的提高. 同时, 通过与其他复原方法作对比可知, 应用本文方法得到的复原结果中含有较少的振铃、噪声等负面效应, 视觉效果相当清晰, 这是明显优于其他方法所得复原结果的.

综合以上实验结果可得, 在泊松噪声污染的模糊共聚焦图像复原方面, 本文提出的方法具有突出优势.

4 结论

本文针对泊松噪声下共聚焦图像的去卷积问题, 提出了一种基于Hessian矩阵范数的正则化方法. 与经典的TV正则化方法相比, 本文所提方法能够在保留图像细节的同时, 消除噪声效应和阶梯效应, 得到更加优良的复原结果, 实验结果证实了其有效性.

而且, 该方法不仅适用于共聚焦显微镜, 也能扩展到其他的快速生物成像系统. 将来, 还可探讨深度学习方法在模型中迭代参数选择的应用.