1 引言

激光共焦扫描显微镜是一种无损观测微观结构的高端大型科学仪器,广泛应用于生物和医学领域^[1]。这种高端激光共焦扫描显微镜主要被Zeiss、Nikon、Olympus、Leica等国外著名品牌垄断,相关关键图像重构算法并未对外公布^[2]。在这样的情况下,自主研发国产激光共焦扫描显微镜,尤其是能够实现高速扫描成像的高端激光共焦扫描显微镜显得尤为重要。

激光共焦扫描显微系统中常用的扫描成像方式分为单程式隔行扫描和往复式逐行扫描两种^[3]。相比于单程式隔行扫描,往复式逐行扫描方式具有成像速度快、图像数据利用率高等优点,是提高激光共焦高速扫描显微系统成像帧速和成像分辨率的理想扫描方式^[4]。重构图像时,为使图像不错位,图像列间相邻像素必须与样品上实际空间点位置对应一致。然而,在对样品进行往复式逐行扫描的过程中,受振镜运动特性以及图像重建过程中相关采样信号匹配程度的影响,振镜扫描正行程与反行程中的数据采集点对应的实际空间点位置无法准确匹配,即数据采集点重组成的图像列间相邻像素数据无法准确匹配。如果采用采集数据直接排列的方式,重构图像会明显存在锯齿样式错位,这种错位严重影响图像质量,大大降低了图像分辨率。

本文在最大程度降低每帧图像扫描时间的基础上,针对高速振镜扫描时图像出现的错位现象,提出了基于形态学梯度的图像错位评价算法,并利用实验验证了该算法的可靠性。

2 高速扫描振镜扫描方式及错位分析

2.1 高速振镜成像轨迹

激光共焦高速扫描显微系统采用二维振镜配合进行平面扫描^[5],其常用的两种振镜搭配方式为X、Y方向均为检流计振镜,或者X方向是扫描速度较快的共振振镜、Y方向为慢速检流计振镜。检流计振镜由磁驱动器控制,由外界提供的模拟控制信号驱动振镜进行步进式摆动,其扫描线速度约为2 kHz;与检流计振镜不同,共振振镜的镜面连接一个可产生机械谐振的调谐棒,当接受一个固定电压时,调谐棒会带着反射镜以固定的振幅、固定的频率往复转动,输入的电压只与转动的幅度相关,其扫描线速度可以达到12 kHz^[6]。下文提及的检流计振镜和共振振镜型号分别为Cambridge Technology公司出产的 Galvanometer Optical Scanner 6230H和CRS 8 kHz Resonant Scanner。

在激光共焦扫描显微成像系统中,振镜扫描方式决定着每帧图像的成像速度^[7]。以线扫描为基础的往复式逐行扫描方式比单程式隔行扫描方式的成像速度快,数据利用率高,其X方向的振镜以固定频率来回摆动,Y方向的振镜在一帧的扫描时间内保持单一方向摆动。所使用的往复式逐行扫描轨迹如图1所示,这种扫描方式对X方向振镜的摆动速度要求较高。为达到更高的系统成像速度,X方向振镜摆动速度需为Y方向振镜摆动速度的好几倍,因此用双检流计振镜方案扫描一幅128 pixel×128 pixel的图像时帧速仅为10 frame/s;该方法的成像速度慢,但是图像错位等失真程度较低。用共振振镜(X方向)搭配检流计振镜(Y方向)方案扫描一幅512 pixel×32 pixel的图像时,帧速可以高达300 frame/s;该方法帧速高且利用了两个不同运动特性的振镜,但图像错位等失真程度极其严重,因此高成像质量需要更为复杂的图像重构技术支撑。

图 1. 高速振镜逐行扫描成像轨迹图

Fig. 1. Traces of high speed mirror reciprocating progressive scanning

下载图片 查看所有图片

2.2 高速振镜错位分析

针对双检流计振镜扫描系统,成像错位的原因主要是检流计振镜对其驱动电压的响应有一定的延迟^[8]。在激光共焦显微成像系统启动时,数据采集指令与控制振镜的驱动电压指令同时发出,这导致数据采集卡在连续采集样本灰度数据初期,由于振镜系统响应延迟,数据信号起始端存入了重复性的无效灰度数据。在图像重构时,首行使用了无效灰度数据进行排列,致使后续每行对应的灰度数据均混入前一行的灰度信息,这是导致双检流计振镜扫描系统图像错位的主要原因。

对于共振-检流计振镜扫描系统,驱动板反馈位置信号与角速度信号不匹配是造成图像错位的主要原因。在共振振镜驱动板输出的信号中,振镜反馈位置信号(即TTL Sync信号)提供用来确定振镜方向的信息,由于其下降沿陡峭,系统用其作为开始采集指令的触发信号,信号的高电平部分代表振镜处于正向行程运动阶段,低电平代表反向行程运动阶段;振镜运动角速度信号是由其角度位置信号求导得出的余弦函数,共振振镜具有与其扫描谐振频率相等的正弦角位移,因此振镜运动角速度信号也作为精确计算振镜运动实际频率的依据。在实际重构图像过程中,每行图像灰度数据必须完全对应振镜运动角速度信号正周期或负周期。然而,由扫描系统硬件导致的相关反馈信号无法完全匹配,如图2(a)所示,致使由TTL同步脉冲信号触发开始采集指令后,由获得的灰度数据直接排列图像时,每行灰度数据无法对应相应的角速度,行间像素排列错乱,产生了图像错位。模拟仿真结果如图2所示。TTL Sync信号与角速度反馈信号不完全匹配的情况如图2(a)所示;对应信号下重构产生错位的艾里斑图像如图2(b)所示,图像中艾里斑呈现锯齿样式,特征区域边缘模糊不清,图像分辨率下降;理想信号反馈图如图2(c)所示;当TTL Sync信号高电平部分完全对应角速度反馈信号正周期时,其对应的重建艾里斑图像如图2(d)所示,从图中可以看出特征区域边缘平滑,图像分辨率较高。因此,综合以上对不同振镜搭配系统的成像错位分析,通过确定有效灰度数据起始位置可以显著校正图像错位。

图 2. 振镜驱动板反馈信号与图像质量关系图。(a)错位信号反馈图;(b)图2(a)对应的重建错位艾里斑图像;(c)理想信号反馈图;(d)图2(c)对应的重建艾里斑图像

Fig. 2. Relationship between feedback signals of mirror driving board and image quality. (a) Dislocation signal feedback; (b) reconstructed dislocation Airy spot corresponding to Fig. 2(a); (c) ideal signal feedback; (d) reconstructed Airy spot corresponding to Fig. 2(c)

下载图片 查看所有图片

3 高速振镜扫描图像错位校正

由成像错位分析可知,图像出现错位的主要原因是扫描振镜系统对驱动电压响应延迟以及其反馈位置信号与角速度信号不匹配,从而导致在利用同步采集的灰度数据重建图像时,无法确定重构图像的有效灰度数据起始位置,重构出来的图像行间数据错乱,因此确定图像有效灰度数据的起始位置(即采集数据错位量)是校正算法的关键。激光共焦高速扫描显微系统的图像错位校正算法结构框图如图3所示。

图 3. 错位校正算法结构框图

Fig. 3. Block diagram of dislocation correction algorithm

下载图片 查看所有图片

处理系统成像错位问题的过程为:首先设置错位量ΔS的变化阈值,根据不同的ΔS构造图像;经基本图像轨迹正弦畸变校正^[9]后,对每一帧图像运用基于形态学梯度的图像错位评价算法得到图像错位评价系数;在实际应用中,选取的ΔS阈值过大会使目标图像完全分离,严重干扰错位评价算法的运算结果。针对这种情况,结合单目标约束优化粒子群算法,搜索图像错位评价系数最小值及其对应的最优错位量ΔS_T,最后应用ΔS_T重建图像,完成图像错位校正。

3.1 图像构建与ΔS的关系

在图像重建中,ΔS直接影响图像每行端点数据的确定。根据光学相关理论公式及共振振镜说明书手册(双检流计振镜扫描系统轨迹参考文献[ 10]),得出采样点在扫描平面的实际横坐标X(i)与采集数据点数i的关系为

X(i)=ftan-2Amaxcos2πfgfs(i-1),i=1,2,…,NS,(1)

式中:f为激光共焦扫描显微镜物镜焦距;A_max为振镜运动的最大角度;f_g为振镜运动频率;f_s为数据采集频率;N_S为一帧图像灰度数据的个数。根据ΔS确定不同的采样点序列i,i∈[ΔS,N_S+ΔS-1],计算得到实际扫描平面的横坐标位置X(i)。将X(i)中邻域极小值作为图像每行左端点数据,确定其所有左端点对应的采集点序列I_S(i_r),i_r=1,2,…,N_r,其中N_r为每帧图像的行数。针对因数据采集频率与振镜运动频率无法完全对应导致每行采集数据点不相等的情况,根据N_l=[f_s/f_g]进行计算,并相应提取出每个左端点I_S(i_r)后N_l-1个数据作为每行排列的像素点,其中N_l为每行像素数据的个数,[]为向上取整符号。在重建过程中,选取不同的ΔS造成首行图像灰度数据不完整,这时需要舍弃首行灰度数据,直接使用第二行端点坐标位置I_S(2)对应的灰度值作为图像起始灰度数据进行图像重建。

3.2 基于形态学梯度的图像错位评价算法

ΔS的取值不同,重建得到的图像错位程度也不同,因此需要获得不同ΔS下的图像错位程度评价系数。这里采用基于形态学梯度的图像错位评价算法对每帧图像进行评价。该算法先利用具有一定特征的结构元素探测一帧图像,进行变换提取梯度,然后利用错位评价函数进行目标图像边缘提取检测,得出对应ΔS下的图像错位评价系数。

3.2.1 结构元素

结构元素的形状和大小决定了形态学运算中所提取的图像信息。形态学有膨胀和腐蚀两种基本运算,灰度图像的形态学定义为^[11]

fd(x,y)=f(x,y)⊕b(m,n)=max{f(x-m,y-n)+b(m,n)|[(x-m),(y-n)∈Df];(m,n)∈Db},(2)fe(x,y)=f(x,y)Θb(m,n)=min{f(x+m,y+n)-b(m,n)|[(x-m),(y-n)∈Df];(m,n)∈Db},(3)

式中:f(x,y)为原始灰度图像,其定义域为D_f;f_d(x,y)和f_e(x,y)分别表示膨胀和腐蚀运算后得到的图像;b(m,n)表示结构元素,其定义域为 Db12。(2)、(3)式分别表达了灰度膨胀和腐蚀的几何意义,即通过将结构元素的原点平移到与信号重合,膨胀运算可得到信号上每一点与结构元素对应各点和的图像灰度区域最大值,腐蚀运算可得到信号在结构元素大小范围内与结构元素各点差的图像灰度区域最小值。在对结构元素的实际选取中,既要考虑能否计算得到显微镜对不同生物细胞的成像上各个方向的灰度梯度,也要避免实验光路模块、电路模块以及环境因素引入各种各样的图像噪声对错位评价函数造成的影响。采用的多结构元素可定义为

b1=010010000,b2=000110000,b3=100010000,b4=001010000,b5=000010010,b6=000011000,b7=000010001,b8=000010100。(4)

应用以上8个方向的结构元素对图像进行变换并提取梯度g_i,将所得梯度进行累加得到图像梯度g=∪ig_i。通过使用上述的结构元素,可以有效地检测出图像错位部分,同时避免孤立噪声点的影响。

3.2.2 错位评价函数

激光共焦扫描显微镜分辨率的评判标准主要来自于图像的边缘和轮廓,因此可通过形态学梯度计算梯度值来提取图像的边缘错位信息^[13]。图像中荧光激发被染料标记的部分相比于未标记部分,边缘区域具有较大的梯度,因此采用的形态学梯度边缘检测算子为

g(x,y)=Fd(x,y)-Fe(x,y),(5)

式中:F_d(x,y)为形态学外部梯度;F_e(x,y)为形态学内部梯度。F_d(x,y)和F_e(x,y)可分别表示为

Fd(x,y)=fd(x,y)-f(x,y),(6)Fe(x,y)=f(x,y)-fe(x,y)。(7)

为衡量图像的错位程度,将形态学梯度边缘检测算子累加,以此作为形态学梯度的图像错位程度评价系数,其表达式为

EMGD=∑xM∑yN|g(x,y)|,(8)

式中:M、N为应用错位校正算法的图像区域尺寸,这里对整幅图像像素数据进行算法处理。错位评价系数E_MGD越小,表示图像错位程度越小;反之,图像错位越严重。基于形态学梯度的图像错位评价算法流程如图4所示。

图 4. 基于形态学梯度的图像错位评价算法流程图

Fig. 4. Flow chart of imaging dislocation evaluation algorithm based on morphological gradient

下载图片 查看所有图片

通过基于形态学梯度的图像错位评价算法可以对每帧图像的错位程度进行准确评估,此时需确定错位评价系数最小的一帧图像及此时的ΔS_T,以完成错位校正。但是,实验表明每次采用相同的ΔS_T进行图像错位校准时,效果并不好,这是由于振镜单次启动延迟时间及相关反馈信号错位程度存在细微的差别,相同的参数不能保证每次获得理想无错位的图像。因此,进一步结合单约束粒子群算法进行ΔS_T的自适应搜索。

3.3 单目标约束优化粒子群算法

根据图像错位评价算法可以计算出选取ΔS时图像的错位评价系数,这时需要用搜索算法确定最小的错位评价系数ΔS_T。经多次实验发现错位量和错位评价系数的关系曲线的单调性明显。因此采用单目标约束优化粒子群算法搜索ΔS_T,算法流程如图5所示。

图 5. 单目标约束优化粒子群算法流程图

Fig. 5. Flow chart of single-objective constrained optimization particle swarm algorithm

下载图片 查看所有图片

搜索算法首先初始化一批粒子,对于其中的每个粒子,根据错位量搜索阈值及速度设置区间,随机生成一个当前错位量和当前速度,根据当前错位量对应的目标函数(即错位评价函数)确定一个评价系数,即E_MGD,每个粒子对应的最小评价系数记为个体最优值E_PMX,整个粒子群中最小的E_PMX记为全局最优值E_GMX,每个粒子根据E_PMX与E_GMX不断变化着错位量和速度的大小,最终得到目标最优错位量。本文算法中,适应度函数设为目标函数,且目标函数的评价系数与错位量一一对应。粒子速度和错位量的更新方程分别为

vi(t+1)=wvi(t)+c1r1(t)[pi(t)-Si(t)]+c2r2(t)[pg(t)-Si(t)],(9)Si(t+1)=Si(t)+vi(t+1),(10)

式中:S_i(t)为粒子当前的错位量;S_i(t+1)为更新一次后粒子的错位量;v_i(t)为当前粒子的移动速度,v_i(t+1)为对应S_i(t+1)的粒子移动速度^[14];c₁和c₂为加速因子,其中c₁为认知因子,c₂为社会因子,为了在搜索中协调粒子群自身认知部分以及社会部分,将c₁和c₂都设定为2^[15];p_i为个体最优错位量;p_g(t)为全局最优错位量^[14];r₁、r₂为[0,1]范围内变换的随机数;w为惯性权重。采用线性调整w策略,其方程为

w=w1-(w1-w2)tT,(11)

式中:w₁为惯性权重上限,w₂为惯性权重下限,经多次实验验证,将w₁与w₂的值分别设为0.9与0.4;t和T分别为当前迭代次数和最大迭代次数。在算法实际应用过程中,由于振镜扫描速度极快,ΔS阈值选取不当会造成重构图像中特征区域完全分离,对错位评价算法造成很大的干扰。针对此情况,引入罚函数作为错位函数的附加约束条件。具体做法是:选取一帧图像最强灰度信号点,标定该点为目标图像特征点。对图像上该点G(x,y)和邻域点G(x,y+1)的灰度值作差,以此差值构造特殊的罚函数,其目标函数可表示为

F(x)=EMGD+P(x,y),(x,y)∈Ω,(12)

式中:P(x,y)=|G(x,y)-G(x,y+1)|为惩罚项,Ω代表整幅图像。在校正错位图像过程中运用单约束粒子群算法的目的在于自动准确地搜索出ΔS_T。

4 错位校正实验

为验证本文算法,系统实验在自主搭建的激光共焦高速扫描显微系统中进行。实验中分别使用双检流计振镜搭配方式以及共振-检流计振镜搭配方式扫描图像,以此验证算法的准确性及自适应特性。

4.1 双检流计振镜扫描算法验证

图6为在激光共焦高速扫描显微系统采用双检流计振镜作为扫描模块条件下,系统采集的老鼠脑神经细胞样本以及通过算法校正错位后的对比图像。图像大小为128 pixel×128 pixel,细胞尺寸约为10 μm,扫描图像的实际空间尺寸约为83.78 μm×83.78 μm,每个像素对应的样品尺寸约为0.655 μm。

图 6. 双检流计振镜扫描成像错位校正。(a)未经算法校正;(b)经算法校正后;(c)图6(a)错位图像中红色线段部分灰度值与像素的关系曲线;(d)图6(b)校正后图像中红色线段部分灰度值与像素的关系曲线

Fig. 6. Imaging dislocation correction in double-glvanometer scanning. (a) Before correction; (b) after correction; (c) gray value and pixel relationship curve of the red line in dislocation image corresponding to Fig. 6(a); (d) gray value and pixel relationship curve of the red line in correction image corresponding to Fig. 6(b)

下载图片 查看所有图片

系统采用双检流计振镜作为扫描模块时,产生图像错位的主要原因是振镜电机响应时间延迟,应用错位校正算法解决了振镜运动与数据采集指令不同时启动产生的图像错位问题。图6(a)所示为ΔS=0(即对采集数据直接排列时)的图像,其对应的E_MGD=7188498;图6(b)为经算法校正后的图像,此时ΔS_T=4.27,E_MGD=7062851。图6(c)、6(d)分别由图6(a)、(b)中红线标出的细胞区域相邻两行灰度数据作图得到,经错位校正算法处理后,细胞宽度由11.89 μm校正为7.20 μm,横向分辨率提高约39.44%。

双检流计振镜系统的成像帧速较慢,因此所得数据直接排列的图像错位现象并不严重,如图6(a)所示。由图6(c)、6(d)可以看出,相邻两行灰度数据的错位幅度适中,没有出现目标图像奇偶行像素完全分离的情况,此时搜索算法中的罚函数对搜索结果的影响较小。ΔS_T的取值受振镜扫描频率的影响较大,由于检流计振镜扫描频率较低,故ΔS_T取值较小。

利用经单目标约束优化粒子群算法搜索的结果绘制曲线,如图7所示,图中离散的点代表了总数为10的粒子在10次迭代中,每次迭代对应ΔS与E_MGD之间的关系,其中ΔS的初始化阈值范围为[-10,10],速度取值范围为[0,5]。算法中粒子的收敛较快,并大量集中在最小评价系数附近,搜索结果可靠。

图 7. 双检流计振镜扫图模式下ΔS与EMGD的关系曲线

Fig. 7. Relationship between ΔS and EMGD in double-galvanometer scanning

下载图片 查看所有图片

4.2 共振-检流计振镜扫描算法验证

图8为在激光共焦高速扫描显微系统中采用共振-检流计振镜作为扫描模块条件下,系统采集的405 nm荧光染料小球以及通过错位校正算法处理后的对比图像,图像大小为630 pixel×630 pixel,荧光小球尺寸约为4 μm,扫描图像实际空间尺寸约为139.1 μm×139.1 μm,每个像素对应的样品尺寸约为0.22 μm。

图 8. 共振-检流计振镜扫描成像错位校正图。(a)未经算法校正;(b)经算法校正后;(c)图8(a)错位图像中红色线段部分灰度值与像素的关系曲线;(d)图8(b)校正后图像中红色线段部分灰度值与像素的关系曲线

Fig. 8. Imaging dislocation correction in resonant-galvanometer scanning. (a) Before correction; (b) after correction; (c) gray value and pixel relationship curve of the red line in dislocation image corresponding to Fig. 8(a); (d) gray value and pixel relationship curve of the red line in correction image corresponding to Fig. 8(b)

下载图片 查看所有图片

与采用双检流计振镜作为系统扫描模块元器件不同,共振-检流计振镜搭配方式的图像产生错位的主要原因为反馈位置信号与角速度信号不匹配,应用错位校正算法解决了图像列间相邻像素数据无法准确匹配的问题。图8(a)所示为ΔS=0下采集数据直接排列的图像,此时对应的E_MGD=181124389;图8(b)为经算法校正后的图像,其对应ΔS_T=10.51,E_MGD=180987951。图8(c)、(d)由红线标出的细胞区域相邻两行灰度数据绘制得到,经过错位校正后,红线标出的荧光小球宽度由8.15 μm校正为2.54 μm,横向分辨率提高约68.83%。

共振-检流计振镜搭配方式较双检流计振镜模式系统的成像帧速高,但所得数据直接排列的图像错位现象严重。可以观察到图8(a)中目标图像完全分离,图8(c)中相邻两行灰度数据错位幅度大,此时搜索算法中罚函数对搜索结果的影响较大。由于共振-检流计振镜扫描频率高,因此ΔS_T取值较大。

利用单目标约束优化粒子群算法的搜索结果绘制的曲线如图9所示,图中离散点代表在共振-检流

图 9. 共振-检流计振镜扫图模式下ΔS与EMGD的关系曲线

Fig. 9. Relationship between ΔS and EMGD in resonant-galvanometer scanning

下载图片 查看所有图片

计扫描振镜成像条件下,总数为10的粒子在10次迭代中ΔS与E_MGD的关系,其中ΔS的初始化阈值为[-15,15],速度取值范围为[0,5]。可以看出,引入罚函数的单目标约束粒子群算法的曲线单调性较好,算法中粒子收敛较快,并大量分布在最小评价系数周围,搜索结果可靠。该研究为实际应用中根据扫描振镜配置选取ΔS阈值提供了参考。

5 结论

基于形态学梯度的图像错位评价函数可以有效反映图像错位程度,结合单目标约束优化粒子群算法可以准确搜索最佳错位量。实验结果表明,所提算法满足高帧频共聚焦显微成像下的图像错位校正,并可有效提高横向分辨率。基于形态学梯度的错位校正算法为下一步高帧频图像的分割提供了基础,并为利用共聚焦显微镜实现高帧频超分辨成像提供了前提条件。