基于重映射的SPRI视频自动找点方法

袁彬峰，汪之又，黄彩霞，罗 宇

(1.国防科技大学 计算机学院，湖南 长沙 410073；2.长沙学院 电子信息与电气工程学院，湖南 长沙 410022；3.长沙学院 计算机工程与应用数学学院，湖南 长沙 410022)

0 引 言

表面等离子共振成像(surface plasmon resonance imaging，SPRI)是基于表面等离子共振(surface plasmon resonance，SPR)的光学成像技术，可基于生物芯片表面折射率的变化测定待测溶液中特定物质的浓度[1]。其测定原理可以被描述为：在SPRI实验中放置在生物芯片表面的蛋白分子与待测溶液蛋白分子之间的结合会引起生物芯片表面折射率的变化[2]并被相机拍摄记录。因此判断待测溶液中是否含有某种物质，只需要查找相机记录的SPRI视频中是否存在亮度变化区域，并定位该区域对应放置的蛋白分子。

由于表面等离子共振成像(SPRI)检测技术具有无需标记、灵敏度高、专属性强、消耗样品少的特点[3]，目前已广泛用于食品安全[4]、环境监测[5]、药物研究[6]等领域。当前SPRI技术在结合微点阵技术后，可同时、原位、实时观测成千上万个无标记样点[1]。若利用传统的人工识别蛋白点方法工作量大且准确度低。目前已有成熟的商业化软件用于半自动定位生物芯片上的蛋白点，例如Array Pro，GenePix Pro，MicroVigene，Dapple and Mapix等等[7]。这些商业软件中通常用椭圆检测算法直接提取SPRI视频图像中的蛋白点[8]，或者在人工干预下先将SPRI视频图像按网格划分，然后在每个网格内检测蛋白点[9]。这些商业化软件工具使用方便但对输入数据的格式和使用范围有严格要求，且不能全自动识别与定位蛋白点。商业软件中使用的定位算法处理步骤繁琐，时间开销大，不支持在性能受限的嵌入式设备上运行。

在SPRI实验中发现从生物芯片制备工具中导出的阵列打印效果图为生物芯片的扫描图，并且利用CCD相机拍摄的SPRI视频是对生物芯片上蛋白点的放置区域的记录，因此阵列打印效果图中放置蛋白阵列区域与CCD相机拍摄的SPRI视频是对生物芯片上放置蛋白点区域的不同成像，阵列打印效果图与SPRI视频每帧图像中的蛋白点之间可能存在空间映射关系。结合重映射的思想，该文提出了基于重映射的自动找点方法。该方法需要先定位阵列打印效果图上的蛋白点，然后将打印效果图中蛋白点的坐标映射为SPRI视频中目标蛋白点的坐标。由于阵列打印效果图成像清晰，相较于成像图噪音干扰较少，且只需通过一次映射即可得到目标蛋白点坐标。因此基于重映射的SPRI视频自动找点算法运行时间短，且无需人工干预，支持在嵌入式设备上运行。在测试实验中，时长为50分钟，帧率为1帧每秒，排列大小为48×48的SPRI视频需要214秒定位并提取每个蛋白点区域的灰度值。

1 实验条件与软件环境

生物芯片的制备过程可以被简述为：将与待测物质发生反应的蛋白溶液从试剂盒中吸取，并按阵列形式固定在特定金属芯片上[10-11]。完成生物芯片的制备后对生物芯片进行扫描，可得到如图1(a)所示的芯片的阵列打印效果图。在SPRI实验前需要将生物芯片固定在流通池的下方，并将生物芯片上的蛋白阵列区域放置在CCD相机的成像区域内，图1(a)中蛋白阵列外围正方形黑框代表成像区域。实验开始后利用单色光源照射生物芯片，激发表面等离子共振现象并将待测溶液通入流通池。在激发表面等离子共振现象后，若通入溶液与芯片上放置蛋白分子发生反应则会引起生物芯片折射率变化并被CCD相机记录。图1(a)中的黑色正方形成像区域在CCD记录成像后会被拉伸、压缩为图1(b)所示的矩形图像。

图1 实现原理

在SPRI实验中，由CCD相机拍摄图像(图1(b))通常为灰度图且含有大量噪音。图1(b)中噪声主要分为两种：空间噪声和时间噪声。空间噪声是指由于光源分布不均匀造成的成像面上光强分布不均匀的情况。成像视频常常出现中间亮高，四周亮度低的状况[12]。时间噪声包括高频随机噪声和长期基线漂移。随机噪声主要是由于光源波动、光电探测器热燥引起，基线漂移通常由环境温度变化，光源波动和探测器热噪声等引起[12]。以上噪音都会影响基于图1(b)进行自动定位算法的准确度，导致自动定位算法的准确度与CCD相机的成像质量相关，当成像质量较低时定位蛋白点的准确度降低。

通过在多次SPRI实验中对比阵列打印效果图与成像视频中蛋白点的成像后发现，生物芯片上同一个蛋白点在不同成像图上的圆心坐标，宽和高之间存在映射关系。具体来说，成像视频中蛋白点相对于打印效果图中蛋白点的宽度会被拉伸，高度会被压缩，圆心坐标之间存在线性变换。具体的映射参数取决于当前CCD相机的拍摄参数。

该文提出的自动找点算法使用C++在Ubuntu上利用Qtcreator开发，开发过程中使用opencv3.6依赖库[13-15]，用于截取图像、边缘检测等图像处理操作。为验证自动定位算法的性能，在树莓派4b(Cortex-A72 1.5 GHz处理器，4 GB)基于Debian的Linux操作系统Raspbian lite(2019-07-10)上运行测试。

2 方法设计

该文提出了基于阵列打印效果图的自动找点算法，该算法不直接对拍摄视频(图1(b))进行处理，而是先对阵列打印效果图进行识别，然后将在效果图上的蛋白点坐标映射为成像视频中蛋白点的坐标。基于以上思路，基于重映射的自动识别方法的处理过程可分为六步，处理流程如图2所示。第一步截取阵列打印效果图中的成像区域。第二步对截取的图像进行边缘检测。第三步将上一步中提取的边缘拟合为椭圆形状并过滤噪音。第四步将第三步中拟合的椭圆组织为椭圆阵列。第五步需要输入CCD成像图，并将上一步中椭圆坐标映射为CCD成像图中目标蛋白点坐标。第六步修正坐标映射后的椭圆点坐标，使映射后的椭圆与CCD图像中的蛋白点重合。该方法中需要输入的阵列打印效果图与CCD成像图分别对应图1(a)与图1(b)。

图2 算法处理流程

2.1 截 取

在这一步中需要截取阵列打印效果图中的视频成像区域即图1(a)中的黑色边框区域。由于阵列打印效果图是对生物芯片的扫描，且在制备生物芯片时必须确定成像区域的起始点与截止点坐标，因此可以使用成像区域的起始点与截止点坐标截取阵列打印效果图中的视频成像阵列区域。

2.2 边缘检测与椭圆拟合

在这一步中需要对阵列打印效果图的截图进行边缘检测，然后将提取的蛋白点边缘拟合为椭圆形状。常规边缘检测可以使用边缘算子或其他方法[16-17]。但对于边缘不清晰的蛋白点成像，可以使用改进的Canny边缘算子进行边缘检测[18-20]。最后可使用opencv库中的fitelipse方法将之前检测的蛋白点边缘拟合为椭圆。

2.3 阵列识别

阵列识别的目标是输出椭圆矩阵，矩阵中每一个元素代表一个椭圆，记录了对应椭圆的大小、圆心坐标等信息。为达到该目标，需要先猜测阵列大小，然后修正上一步中拟合的椭圆，去除非目标椭圆，并将目标椭圆按阵列方式组织。最后检查当前阵列是否与蛋白点的排列大小相同，若不相同则需要通过阵列修正操作补全缺失椭圆。

2.3.1 构建阵列与过滤

为构建阵列，首先需要将上一步中拟合椭圆点按圆心纵坐标分组。可依据两个椭圆点圆心的纵坐标之间的差值来判断两个椭圆是否处于同一行。若差值接近于0，则说明这两个椭圆点位于同一行。基于以上思路可将上一步中拟合的椭圆点按行分组，最后分组的数量即为阵列的行数。

完成按行分组后，则需要对每组内的椭圆进行初步过滤。初步过滤中需要过滤的错误椭圆可以被分为两类：第一类是拟合阵列打印图中小黑点噪音得到的椭圆，这类椭圆面积远小于正常椭圆且接近于0，可以通过设置面积的阈值进行过滤。第二类是蛋白点边缘不清晰导致拟合一个蛋白点的边缘得到的多个椭圆。为过滤这类椭圆，需要对组内椭圆按圆心的横坐标排序。若两个椭圆圆心横坐标相近，则可判定这两个椭圆为同一个蛋白点的不同拟合结果，则只需要保留其中一个。

完成以上构建操作后可以得到椭圆阵列，阵列中每个元素表示一个椭圆。但该阵列的列数还未确定，每行中的椭圆数量可能不一致。为求椭圆阵列的列数需要统计阵列中每行的椭圆数量并排序，选择重复次数最多的数量作为当前阵列的列数。

椭圆阵列构建完成后，还需要对阵列内椭圆进行第二次过滤。第二次过滤的目标是去除拟合错误导致的偏大椭圆。这类椭圆是因为在拟合时错误将多个蛋白点拟合为一个椭圆点造成的，对于这类椭圆点可以通过设立椭圆的宽度的上限[imageWidth/cols]和高度上限[imageHeight/rows]进行过滤。imageWidth与imageHeight为CCD视频图像的宽度与高度，rows与cols为之前步骤中确立的椭圆阵列的行数与列数。

2.3.2 阵列修正

构建椭圆阵列后需要检查当前阵列大小与2.3.1节中计算得到的阵列模式是否一致。若发现某行缺失椭圆，则可以通过椭圆圆心之间的距离差定位缺失椭圆点编号。在正常情况下，各个椭圆圆心之间的距离是固定的，但会在某个小范围内波动。如果两个椭圆点中心之间的差距过大，则可以判定两个椭圆点之间可能存在遗漏椭圆点。如图3所示，ΔX5与ΔX8明显大于ΔX1，ΔX2等正常间隔之间的距离。由此可以判断5号与6号椭圆之间，8号椭圆与9号椭圆之间存在缺失椭圆。

图3 补全椭圆

定位缺失椭圆点的编号后，可以依据缺失椭圆点左右两侧的椭圆点补全缺失椭圆点的大小和圆心坐标。如图3所示，5号椭圆与6号椭圆之间缺失椭圆点A的宽度等于5号椭圆点与6号椭圆宽度的平均值，同理可计算缺失椭圆点A的高度和圆心Y轴坐标。缺失椭圆点A圆心的X轴坐标等于5号椭圆圆心的X轴坐标加上ΔX1，ΔX2，ΔX3，ΔX4的平均值。补全A号椭圆后，按照之前定义的步骤依次补全B和C号椭圆。

2.4 坐标映射

其中，offsetx、offsety已在图1中标注，表示蛋白阵列中左上角蛋白点距离起始点在横轴和纵轴上的差值。ratex、ratey表示CCD成像图相对于打印效果图在横轴和纵轴上的压缩比。以上参数的取值和CCD相机参数相关，在SPRI实验中保证拍摄参数不变则以上参数不变。

2.5 自动校准

映射后的椭圆点可能与CCD成像图中的蛋白点轮廓不能完全吻合。为此需要对映射后的椭圆点进行自动校正。在实验中CCD相机采集灰度图像，因此在成像图中目标蛋白点呈现白色。在理想情况下，映射后的椭圆点与蛋白点轮廓重合，椭圆点所在区域内为白色，椭圆点的重心与圆心重合。当映射效果不理想时，重心与圆心之间的差值代表应该修正的方向。此时用重心替换椭圆点圆心便可得到修正后的椭圆点。椭圆点的重心计算公式如下：

其中，xc与yc表示重心的横坐标与纵坐标，M00表示图像的0阶矩，M10与M01表示图像的一阶矩，其计算方式如下所示：

其中，V(i,j)表示图像V像素点(i,j)处的灰度值。

在自动校准后，可以将灰度图转换为二值化图，然后使用椭圆点区域内白色像素点数量与该区域内像素点数量的比值衡量椭圆点的定位准确度。在理想情况下白色像素点数量应该约等于该区域内像素点数量，比值接近于1说明定位准确度高，比值接近于0说明当前准确度低。

3 实验结果与性能评价

通过实验验证，使用文中的自动找点方法可以准确定位视频成像图中的每个蛋白点。图4左侧是对48×48排列模式的成像视频使用本算法识别后的局部截图。图中白圈代表预测的蛋白点位置，椭圆形白斑为蛋白点的成像。对比后可以发现白圈基本能圈中每个蛋白点。

在SPRI实验中，由于相机的曝光与聚焦设置会导致CCD相机在实验开始的一段时间内拍摄图像接近于空白。此时若使用传统的自动定位方法将不能定位蛋白点。但由于本方法是对打印效果图进行分析，然后将分析结果映射变换到CCD成像图中。通过这种模型，只要保证映射关系准确，可以消除视频图像质量低、杂音多的影响。例如在图4右图中成像图亮度太高使蛋白点轮廓不清晰的前提下，也能定位每个蛋白点的位置。

图4 找点结果

为验证自动找点算法的性能，该文提出的自动找点算法在树莓派4b的硬件设备上运行测试。对于时长为50分钟，帧率为1帧每秒，排列大小为48×48的SPRI视频，需要花费214秒定位所有蛋白点的位置并提取每帧图像中每个蛋白点区域的平均灰度值。

使用自动找点算法定位的蛋白点为目标蛋白点的成像。图5是利用自动找点算法识别的一个蛋白点在SPRI实验中反应强度的变化曲线。图中横轴是时间编号(图像帧的编号)，纵轴是该蛋白点的反应强度(蛋白点成像区域内的平均值灰度值减去本次实验的空白信号)。在一次完整SPRI实验中会依次通入清洗液，生理盐水，甘油，清洗液，生理盐水，测试溶液。不同的溶液会导致该蛋白点的放置区域内出现不同的反应强度。通入清洗液时会破坏蛋白质分子之间的结合导致蛋白点放置区域内的折射率降低，反应强度降低。生理盐水是缓冲剂，通入生理盐水时当前蛋白点的反应强度代表本次实验的空白信号。通入甘油时蛋白点的折射率达到最高，通入甘油用于测试和校准实验。在实验最后会通入样品溶液，若发生蛋白分子间的结合反应则会导致折射率明显提高，使用此时的强度值与通入生理盐水时强度值之间的差值对照查表，可得样品溶液中与该蛋白点发生反应物质的浓度含量。最终该蛋白点在一轮实验中的强度值变化如图5所示。

图5 局部强度变化曲线

4 结束语

该文对如何定位SPRI视频中的蛋白点提出了一种新的解决思路。该算法不直接对含有较多噪音的SPRI视频进行分析，而是通过生物芯片的打印效果图间接提取生物蛋白点坐标。因此能适应各种蛋白点的排列模式，适用于噪音较多并且模式固定的SPRI视频分析，并且支持移植到嵌入式设备内运行。但是该算法的时间开销存在优化空间，这也是今后的研究方向。