CN109816609B

CN109816609B - 一种基于傅里叶变换的数字pcr图像还原方法及其应用

Info

Publication number: CN109816609B
Application number: CN201910099525.8A
Authority: CN
Inventors: 高鹏; 夏江
Original assignee: Pilot Gene Technologies Hangzhou Co ltd
Current assignee: Pilot Gene Technologies Hangzhou Co ltd
Priority date: 2019-01-31
Filing date: 2019-01-31
Publication date: 2021-12-21
Anticipated expiration: 2039-01-31
Also published as: CN109816609A

Abstract

本发明涉及图像处理技术领域，具体的涉及一种基于傅里叶变换的数字PCR图像还原方法及其应用。所述基于傅里叶变换的数字PCR图像还原方法包括以下步骤：S1.获取待处理的图像；S2.图像的预处理；S3.设计滤波器；S4.频域图像滤波；S5.计算矩阵校正系数；S6.输出校正后图像。本发明公开的基于傅里叶变换的数字PCR图像还原方法，其解决了光照不均一导致的图片中心区域像素高、边缘像素低的问题，为后续依据荧光信号强弱判断液滴的阳性或阴性提供准确的数据。

Description

一种基于傅里叶变换的数字PCR图像还原方法及其应用

技术领域

本发明涉及图像处理技术领域，具体的涉及一种基于傅里叶变换的数字PCR图像还原方法及其应用。

背景技术

数字PCR技术是指DNA或RNA样本被稀释、均分到独立的微反应单元，通过特定基因片段标记的荧光探针进行标记，进行分子模板PCR扩增，扩增完成后依据荧光信号强弱判断微反应单元的阳性或阴性的方法。荧光信号强弱是用光学成像系统检测。光学镜头在成像时，会发生光照不均一的现象，导致视场中心区域成像正常，视场边缘成像亮度暗。这种现象是由透镜的特征决定的，即球面透镜的中心区域透过率高，距离中心越远透镜率越低，导致成像系统拍摄的图像中心区域亮，视场边缘区域暗，严重影响图像质量。图1是数字PCR光学成像系统示意图，LED光源穿过透镜后形成平行光线，滤光片的作用是允许特定波段的光穿过，在芯片出反射后在CCD处成像。

现有的特殊处理的透镜的制作工艺比较复杂，需要更改现有的LED光源，制作成本加重；另外，图像的荧光值强度与液滴阴阳性、试剂浓度、芯片材质存在关系，多项式函数关系是理想条件；由于人为因素使透镜安装倾斜，导致光路中心偏移了图片的中心，函数的对应关系发生变换，现有的技术只能人工重新校对透镜的位置，重新拍照和计算，造成时间和资源的浪费。

现有的解决思路通常有两种：

方案一：使用特殊处理的透镜(类似于非球面)，改进透镜光学性能使透光率更均一。

方案二：假设导致视场边缘灰度值偏低的原因只与成像系统有关(与液滴的阴阳性、试剂浓度、芯片材质等无关系)。荧光值强度与距离图像中心远近符合一定的函数关系。制作一个标准芯片，拍摄标准芯片图片，计算不同距离(与中心位置)的像素变化，并用多项式函数拟合这种变化，将同一个成像系统拍照的图片带入到多项式中解决光能损失引起的光路不均一问题。图2显示不同距离下获取的补偿系数(标准片中心像素值和其他位置像素值的比值)。

发明内容

本发明公开了一种基于傅里叶变换的数字PCR图像还原方法，其解决了现有技术中的不能很好解决透镜成像系统的光照不均一的问题，为后续依据荧光信号强弱判断液滴的阳性或阴性提供准确的数据。

术语解释：

图像傅里叶变换：将图像从空间域转换到频率域(即图像的灰度分布函数变换为图像的频率分布函数)，其逆变换是将图像从频率域转换到空间域。

数字PCR(即Digital PCR)技术：是一种核酸分子绝对定量技术，准确地、可信地检测稀有基因的突变的技术。

滤波器，是指对波进行过滤的器件，其可以对电源线中特定频率的频点或该频点以外的频率进行有效滤除，得到一个特定频率的电源信号，或消除一个特定频率后的电源信号。

本发明的具体方案如下：

本发明一方面公开了一种基于傅里叶变换的数字PCR图像还原方法，包括以下步骤：

S1.获取待处理的图像；

S2.图像的预处理；

S3.设计滤波器；

S4.频域图像滤波；

S5.计算矩阵校正系数；

S6.输出校正后图像。

应该理解，本发明不限于上述步骤，还可以包含其他的步骤，例如在步骤S1之前、步骤S1和S2之间、步骤S2和S3之间、步骤S3和S4之间、步骤S4和S5之间、步骤S5和S6之间、步骤S6之后，还包含其他额外的步骤，而不超出本发明的保护范围。

在发明的一优选实施例中，在S1中，通过CCD感光元件获取芯片区域的数字图像。

在本发明的一优选实施例中，在S2中，通过高斯滤波或者直方图均衡化方法进行图像的预处理。

在本发明的一优选实施例中，所述步骤S3包括：

S31.第一部分滤波器的设计：

其中C₁＝0.1，W₁＝5，HEIGHT和WIDTH分别表示图像的高度和宽度，F1(i，j)表示图像位于点(i，j)的滤波值；

S32.第二部分滤波器的设计：

其中C₂＝5，W₂＝0.1，HEIGHT和WIDTH分别表示图像的高度和宽度，F2(i，j)表示图像位于点(i，j)的滤波值；

S32.总的滤波器的设计：

其中H₁＝H₂＝1，L₁＝L₂＝0.1，F(i，j)表示图像位于点(i，j)的滤波值。

应该理解，在上述步骤之前、之间、之后，还包含其他额外的步骤，而不超出本发明的保护范围。

在发明的一优选实施例中，所述步骤S4包括：

S41.图像的傅里叶变换。

在发明的一优选实施例中，所述步骤S41包括：

S411.计算傅里叶变换的扩展宽度Col_min＝2^p+3^q+5^r，其中p，q，r为大于0的整数，同时Col_min为大于WIDTH的最小整数值；扩展长度Row_min＝2^p+3^q+5^r，其中p，q，r为大于0的整数，同时Row_min为大于HEIGHT的最小整数值；HEIGHT和WIDTH表示图像的高度和宽度；

S412.创建两个大小为Row_min*Col_min的图片，记作Re和Im；设P^back(m，n)表示步骤S2中预处理之后的图像，将P^back(m，n)对应位置的像素值赋予Re，Re其它位置的像素值赋予0，Im所有位置的像素值赋予0；

S413.将P^back(m，n)带入傅里叶变换公式

其中P^back(m，n)表示步骤S2中预处理之后的图像，FFT(k，l)为傅里叶变换后的结果，傅里叶转换的结果是复数，是实数图像和虚数图像的叠加，用Re和Im表示；

S414.计算幅度图像

其中0≤m≤Row_min，0≤n≤Col_min；对图像Mag(m，n)取对数，记为logM；

S415.以(0，0，cx，cy)、(cx，0，cx，cy)、(0，cy，cx，cy)和(cx，cy，cx，cy)划分幅度图像Mag，其中

将左上角区域和右下角区域交换，右上角区域和左下角区域交换组成幅度图像Magswap；

S416、滤波处理：将幅度图像Mag_swap带入S3设计的滤波器F中，得到滤波后的图像B_filer＝Mag_swap*F，其中F为步骤S3得到的滤波器，F需要按步骤S411-S412进行尺度扩展；

S417、图像进行归一化处理：通过公式计算，所述公式为

其中0≤y≤Row_min，0≤x≤Col_min，

其中P_min和P_max为图像B_filer(x，y)的最小值和最大值；

S418.B_filer(x，y)图像进行逆变换得到逆变换图像P_inverse，同时尺度缩小为HEIGHT*WIDTH；

S419、图像进行归一化处理：通过公式计算，所述公式为

其中0≤y≤HEIGHT，0≤x≤WIDTH，其中P_min和P_max为图像P_inverse(x，y)的最小值和最大值。

在本发明的一实施例中，所述步骤S5包括：

S51.设傅里叶逆变换图像为I_ifft(i，j)，对逆变换图像为I_ifft(i，j)进行高斯滤波和均值滤波，其中I_ifft(i，j)表示逆变换图像(i，j)处的数值；

S52.任选逆变换图像上一点(i，j)，以(i，j)为中心，R为半径，计算包含区域内所有像素的平均值Mean_ifft(i，j)，其中，R表示像素个数；

S53.遍历所有的像素位置，计算出最大平均值

S54.任选逆变换图像上一点(i，j)，以(i，j)为中心，R为半径，计算包含区域内所有像素的平均值Mean_ifft(i，j)；

S55.遍历所有的图像像素位置，计算出最小平均值

S56.计算出标准值

S57.计算矩阵校正系数

其中0≤i≤HEIGHT，0≤j≤WIDTH，HEIGHT和WIDTH表示图像的高度和宽度，I_ifft(i，j)表示逆变换图像(i，j)处的数值。

在本发明的一具体实施例中，在S52中，R＝100。

在本发明的一优选实施例中，所述步骤S6包括：

S61.输入图像为P(i，j)，其中P(i，j)表示图像(i，j)处的像素值，校正后图像为Corr(i，j)＝P(i，j)*W(i，j)，其中0≤i≤HEIGHT，0≤j≤WIDTH，HEIGHT和WIDTH表示图像的高度和宽度，W(i，j)表示(i，j)处校正矩阵系数；

S62.输出校正后图像Corr(i，j)。

本发明第二个方面公开了上述的方法在图像处理技术领域中的应用。

与现有技术相比，本发明具有以下有益效果：

(1)本发明公开的技术方案利用图像处理的方法解决了光路均一性的问题，无需更换现有的成像系统的透镜，大大节约了成本且缩短了开发的周期。

(2)本发明设计的滤波器只保留了频域的低频信息，减少液滴阴阳性和试剂浓度对光路均一校正的影响。

(3)及时人为因素使透镜安装倾斜，本发明的技术方案仍可以用图像进行光路均一性校正，无需人工重新校对透镜的位置，避免了时间和资源的浪费。

附图说明

图1是数字PCR光学成像系统示意图；

图2是不同距离下获取的补偿系数的示意图；

图3是本发明实施例中基于傅里叶变换的数字PCR图像还原方法的流程示意图；

图4是本发明实施例中中心化处理的示意图；

图5是本发明实施例中校正前后图像横向投影像素均值与坐标对应关系的示意图；

图6是本发明实施例中数字PCR的平面图；

图7为图6光路均一校正之后的图像；

图8是本发明实施例中数字PCR的三维示意图；

图9是本发明实施例中校正前图像X-Z平面投影；

图10是本发明实施例中校正后图像X-Z平面投影。

具体实施方式

下面结合附图和实施例对本发明的技术方案进行详细描述，但并不因此将本发明限制在所述的实施例范围之中。

一、在本发明的实施方式中，公开了一种基于傅里叶变换的数字PCR图像还原方法，包括以下步骤：

S1.获取待处理的图像；

S2.图像的预处理；

S3.设计滤波器；

S4.频域图像滤波；

S5.计算矩阵校正系数；

S6.输出校正后图像。

具体的技术方案流程图如图3所示。

二、在本发明的实施方式中，在S1中，通过CCD感光元件获取芯片区域的数字图像。

三、在本发明的实施方式中，在S2中，通过高斯滤波或者直方图均衡化方法进行图像的预处理。

四、在本发明的实施方式中，步骤S3包括：

S31.第一部分滤波器的设计：

S32.第二部分滤波器的设计：

S32.总的滤波器的设计：

五、在本发明的实施方式中，所述步骤S4包括：

S41.图像的傅里叶变换。

六、在本发明的实施方式中，所述步骤S41包括：

S413.将P^back(m，n)带入傅里叶变换公式

其中Pback(m，n)表示步骤S2中预处理之后的图像，FFT(k，l)为傅里叶变换后的结果，傅里叶转换的结果是复数，是实数图像和虚数图像的叠加，用Re和Im表示；

S414.计算幅度图像

将左上角区域和右下角区域交换，右上角区域和左下角区域交换组成幅度图像Mag_swap，如图4所示；

S416、滤波处理：将幅度图像Mag_swap带入S3设计的滤波器F中，得到滤波后的图像B_filer＝Mag_swap*F，其中F为步骤S3得到的滤波器，F需要按步骤S411-S412进行尺度扩展，和得到图像Re的方法一样；

S417、图像进行归一化处理：通过公式计算，所述公式为

其中0≤y≤Row_min，0≤x≤Col_min，其中P_min和P_max为图像B_filer(x，y)的最小值和最大值；

S419、图像进行归一化处理：通过公式计算，所述公式为

七、在本发明的实施方式中，所述步骤S5包括：

S53.遍历所有的像素位置，计算出最大平均值

S55.遍历所有的图像像素位置，计算出最小平均值

S56.计算出标准值

S57.计算矩阵校正系数

在一个具体实施方式里面，在S52中，R＝100

八、在本发明的实施方式中，所述步骤S6包括：

S62.输出校正后图像Corr(i，j)。

具体的，选择一张图片(HEIGHT*WIDTH)横向投影，

表示纵坐标(每一列像素之和的均值)，j表示横坐标0≤j≤WIDTH。统计校正前后的像素变化绘制图如图5所示，校正前像素的范围为[47,80]；校正后像素的范围在[57,65],校正后的像素更平缓，均一性更好。

选择一张图片(HEIGHT*WIDTH)如图6所示，图7为光路均一校正之后的图像，整体图像的均一性得到了很好的改善。以图像的宽为X轴，图像的高为Y轴，对应位置的像素值为Z轴，绘制一幅三维图像如图8所示。图9为校正前图像X-Z平面投影，像素的范围为[13,48]，像素在一些区域突变，中心区域的像素明显高于边界区域，光路不均一导致了这种情况的发生；图10为校正后图像X-Z平面投影，像素的范围为[40,57]，像素缓慢变化，中间区域与边界区域像素近似，很好解决了光路不均一导致的中间区域像素高、边界区域像素低的问题。

上述实施例为本发明较佳的实施方式，但本发明的实施方式并不受上述实施例的限制，其他的任何未背离本发明的精神实质与原理下所作的改变、修饰、替代、组合、简化，均应为等效的置换方式，都包含在本发明的保护范围之内。