本发明属于数字全息成像技术与微流控芯片技术领域，涉及一种使用微流控芯片进行细胞全息重建观测的流路结构去除方法。
背景技术：
微流控芯片实验室是本世纪一项重要的科学技术。微流控芯片实验室简称微流控芯片或芯片实验室，指的是在一块几平方厘米的芯片上构建的化学或生物实验室。它把化学和生物等领域中所涉及的样品制备、反应、分离、检测，细胞培养、分选、裂解等基本操作单元集成到一块很小的芯片上，由微通道形成网络，以可控流体贯穿整个系统，用以实现常规化学或生物实验室的各种功能。
从物理上说，微流控芯片是一种操控微小体积的流体在微小通道或构建中流动的系统，其中通道和构件的尺度为几十到几百微米，承载流体的量通常仅为10-9至10-18L。
从工作原理来讲，微流控芯片实验室各个操作单元通过微通道网络内流体的流动相互联系。设计适当结构的微通道网络，使流体在该微通道网络内按一定的方式流动，微流控芯片就可以从整体上实现特定的功能。
但就观察手段而言，对常规微流控芯片内细胞而言，仍无法摆脱荧光标记、染料染色等有损伤标记方法。且对于尺寸较大的芯片与通道较长的流路，光学显微镜难以在成像视野与成像精度间权衡，尤其对于追踪观察的情况需要实验人员不断调整视野，不仅操作繁琐而且极易造成目标丢失。因此如果能够保证成像精度与细胞活性的前提下扩大观察视野，那么无论是对于简化实验流程还是便于后续处理分析都具有重大意义。
数字全息是一种精度可达微米级别的成像方法，由CCD或CMOS等光学元器件直接记录获取样本图像，光路极其精简，易于集成和小型化。同时其具有全视场、非接触、无损伤、实时性、定量化的优点。加之无需对样本染色即可清晰成像，故特别适合于活体生物样品的定量三维重建和快速跟踪，得以在生物医学应用领域尤其是细胞培养观测中发展迅速。此外由于全息图视野直接对应于光学元器件的尺寸，所以全息图具有得天独厚的大视野成像的优点。在同时能够对微米级别物体准确成像的基础上，全息图视野下所记录的物体信息远丰富于常规光学显微镜。
但是由于全息成像的记录过程是由光源发出的光线照射到样本上，样本表面产生的散射光(即物光)和参考光在光敏电子成像器件表面发生干涉，从而形成全息图。而芯片内流路结构也会与光产生强烈的衍射作用，同时由于具有明显的深度结构，所以在进行常规重建时会带来难以克服的共轭像，甚至会严重影响内部细胞成像。
技术实现要素：
有鉴于此，本发明的目的在于提供一种使用微流控芯片进行细胞全息重建观测的流路结构去除方法，实现对微流控芯片内细胞实时追踪，简化观察步骤，提高实验人员观察分析效率。
为达到上述目的，本发明提供如下技术方案：
使用微流控芯片进行细胞全息重建观测的流路结构去除方法，该方法包含如下步骤：
S1：读入待重建的细胞全息图，进行预处理，将细胞全息图从RGB彩色空间映射为灰度图像，对应像素映射关系：
Gray(a,b)＝0.229×R(a,b)+0.587×G(a,b)+0.114×B(a,b)
其中，a、b为行列坐标，R、G、B分别为三原色通道信息；
S2：设置图像分块尺寸，将全息图进行半窗长滑动分块处理；
S3：对于某一分块处理后全息图像块H在记录距离、重建波长下进行卷积重建，得到复矩阵；
S4：对重建所得复矩阵取模，得实数矩阵，经过灰度拉伸与归一化取整得到初步重建图像块，取各初步重建图块交叠部分对应像素点灰度取均值作为整合重建图中的像素灰度完成所有重建图块的整合；
S5：对于整幅整合重建图再次分块，使用最大类间方差设置灰度阈值得到整合重建图块内流路结构可能的成像区域；
S6：结合梯度判断与形态学滤波手对流路结构进行进一步处理，排除掉观察用目标物体，获得流路结构的成像区域；
S7：将成像区域所对应的全息图重建所得的复矩阵中流路结构对应区域的数值保留，将其他位置的数据实部、虚部分别替换为矩阵各元素实部、虚部的均值，构建复数矩阵；
S8：按照卷积重建步骤反向执行，反衍射出复数矩阵对应的衍射信息矩阵holo；
S9：将成像区域所对应的全息图中的流路结构的衍射信息进行抵消，对抵消后的全息图进行卷积重建，得到纯净的细胞物体重建图；
S10：将细胞物体重建图按照对应全息图块分块顺序重新排列，得到芯片流路内所有细胞分布的完整重建图。
进一步，所述卷积重建具体为：
将全息图中衍射信息再现后对应光波表示为：
其中，
式中，j表示虚数单位，λ为重建光波长，θ为衍射角，zi为观察面与记录面距离，xi为重建图中像点横坐标，yi为重建图中像点纵坐标，x为全息图中像点横坐标，y为全息图中像点纵坐标，C(x,y)为是重建光波在传感器平面上的分布，H(x,y)为图像传感器所记录的全息图光强分布；
令C(x,y)的共轭为R(x,y)：
R(x,y)＝exp[-jk(xsinα+ysinβ)]
式中，α、β分别为平面参考光波与空间yoz平面和xoz平面的夹角，k为波数，k＝2π/λ；
根据线性系统理论，将再现图像复振幅表示为：
其中，k＝2π/λ，g为自由空间脉冲响应；
由卷积理论，使用傅里叶变换计算得到再现图像：
Z＝FT-1[FT(H·C)·FT(g)]
式中，FT()表示傅里叶变换，FT-1()表示逆傅里叶变换。
进一步，步骤S2具体包含：
S21：设图像分块边长为L，整幅图像的行列值为M和N，对应整除取图像的行、列所能排布图像块的个数分别为Row和Col；
Row＝M|L
Col＝N|L
S22：取整幅图像左上Row×L行与Col×L列区域为有效重建区域，进行图像裁剪；
S23：以L为分块边长对全息图有效区域进行滑动分块，横向、纵向滑动步长为L/2，直至各图像块平铺排布覆盖整个有效区域，且相邻块间均有四分之一面积区域交叠。
进一步，步骤S3中所述复矩阵为：
U(xi,yi)＝R(xi,yi)+O(xi,yi)×i
其中，R与O分别为复矩阵的实部矩阵与虚部矩阵系数，i表示虚数单位，(xi,yi)为重建图在再现图像平面上的位置分布。
进一步，步骤S4包含如下步骤：
S41：将复矩阵U(xi,yi)取模，获得实数矩阵I(xi,yi)；
S42：对实数矩阵I(xi,yi)进行归一化和取整处理，获得元素数值为整数且拉伸分布在整个图像灰度级[0,level]区间的实整矩阵Is(xi,yi)，
其中，max和min分别代表矩阵I(xi,yi)中的最大、最小元素数值，level取255，round的取整规则为四舍五入取整，Is(xi,yi)为重建图像块。
S43：对于各重建图块交叠部分对应像素点灰度取均值作为整合重建图中像素灰度，完成重建图块整合。
进一步，步骤S5包含如下步骤：
S51：将整合重建图记为Rec，通过图像分块边长L，再次对整合重建图Rec进行分块处理，获得平铺无交叠的整合重建图块rec；
S52：对于整合重建图块rec进行基于最大类间方差的灰度阈值计算，
对于整合重建图块rec，记T为前景与背景的分割阈值，前景点数占图像比例为w0，平均灰度为u0；背景点数占图像比例为w1，平均灰度为u1，图像的总平均灰度为u，g为前景和背景图象的方差，则有：
u＝w0×u0+w1×u1
g＝w0×(u0-u)2+w1×(u1-u)2
则：
当方差g最大时，认为此时前景和背景差异最大，分割阈值T即是最佳灰度阈值；
S53：根据所得到的分割阈值T，对整合重建图块rec进行二值化处理，得二值图像recBW，
其中，x,y为整合重建图块的位置坐标，rec(x,y)表示坐标x,y处的灰度值。
进一步，步骤S6包含如下步骤：
S61：使用Sobel算子对二值图像recBW中流路结构可能区域所有边缘点坐标对应的整合重建图块rec中的像素点进行梯度值计算，
Gx＝[2rec(x,y+1)-2rec(x,y-1)]+[rec(x-1,y+1)-rec(x-1,y-1)]+[rec(x+1,y+1)-rec(x+1,y-1)]
Gy＝[2rec(x+1,y)-2rec(x-1,y)]+[rec(x+1,y-1)+rec(x-1,y-1)]+[rec(x+1,y+1)-rec(x-1,y+1)]
该点对应梯度值为：
S62：计算recBW中所有边缘点对应整合重建图块rec中像素点梯度的均值Gavr；
其中，num为整合重建图块rec中连通域边缘点总数，edge整合重建图块rec中连通域边缘点合集；
S63：设recBW中共初检测出k个连通域，计算其各连通域边缘对应rec中像素灰度梯度均值Gt，t∈1,2,…,k；
S64：进行形态学滤波，根据实际观察细胞尺寸设置滤波范围，实物面积对应在重建图上像素点个数n，经过1.5倍放大，去除所含像素点数少于1.5n的连通域，将此时的二值矩阵为M。
进一步，步骤SS7具体包含，
S71：计算重建所得复矩阵的实部R与虚部系数O的均值，分别记为real和imag，
S72：根据二值矩阵M的逻辑值分布，构建针对流路成像区域的复数矩阵S(xi,yi)，
其中，M(xi,yi)为二值矩阵M的第i行第i列元素。
进一步，步骤S9包含如下步骤：
S91：计算反衍射矩阵holo元素均值holoavr，
S92：从全息图像块H中抵消掉衍射环矩阵holo所含信息，得到仅含有细胞衍射信息的全息图像块H*，
H*(x,y)＝H(x,y)-2[holo(x,y)-holoavr]
S93：对矩阵H*进行卷积重建，得重建复矩阵U*，将重建复矩阵U*取模，得到U*对应的实数矩阵I′，并对实数矩阵I′进行归一化和取整处理。
本发明的有益效果在于：本发明方法所解决的根源问题在于，在使用全息成像技术对微流控芯片流路中的细胞进行重建再现时，由于微流控芯片自身的流路结构等相较细胞尺寸较大，进而重建后会产生非常明显的共轭像，以致严重干扰影响到观测的主体——流路中的细胞。本发明对于实现微流控芯片内细胞实时追踪，简化观察步骤，提高实验人员观察分析效率具有重大意义。
附图说明
为了使本发明的目的、技术方案和有益效果更加清楚，本发明提供如下附图进行说明：
图1为本发明流程图。
具体实施方式
下面将结合附图，对本发明的优选实施例进行详细的描述。
如图1所示，本发明方法包含如下步骤：
S1：读入待重建全息图，进行预处理。
具体包括：
将全息图从RGB彩色空间映射为灰度图像，对应像素映射关系：
Gray(a,b)＝0.229×R(a,b)+0.587×G(a,b)+0.114×B(a,b)
其中a、b为行列坐标，R、G、B分别为三原色(红、绿、蓝)通道信息。
S2：设置分块尺寸，将全息图进行半窗长滑动分块处理。
具体包括：
S21：设图像块分块边长为L，整幅图像行列值为M和N，单位为像素(pixel)。对应整除取其行、列所能排布图像块的个数Row和Col：
Row＝M|L
Col＝N|L
S22：取整幅图像左上Row×L行与Col×L列区域为有效重建区域，完成图像裁剪。
S23：以L为分块边长对全息图有效区域进行滑动分块，横向、纵向滑动步长为L/2，直至各图像块平铺排布覆盖整个有效区域，且相邻块间均恰有四分之一面积区域交叠。由于对于每一图像的处理均相同，故在下述步骤重以H统一代指某一图像块。
S3：对于某一分割后全息图像块H在记录距离、重建波长下进行卷积重建，重建结果得到复矩阵U。
U(xi,yi)＝R(xi,yi)+O(xi,yi)×i
其中R与O分别为复矩阵的实部矩阵与虚部矩阵系数，i在此表示虚数单位，(xi,yi)为重建图在再现像平面上的位置分布。
S4：对重建所得复数矩阵取模，得实数矩阵，经过灰度拉伸与归一化取整得到初步重建图像块。对于各重建图块交叠部分对应像素点灰度取均值作为整合重建图中像素灰度完成重建图块整合。
具体包括：
S41：将S3中重建得到的复数矩阵U取模，得实数矩阵I。
S42：对实数矩阵进行归一化和取整处理，得元素数值为整数且拉伸分布在整个图像灰度级[0,level]区间的实整矩阵Is
其中max和min分别代表矩阵I中的最大、最小元素数值。level通常为255。round的取整规则使用四舍五入取整。至此Is为重建图像块。
S43：对于各重建图块交叠部分对应像素点灰度取均值作为整合重建图中像素灰度完成重建图块整合，记得到的整幅有效区域重建图为Rec。
S5：对于整幅重建图再次分块，使用最大类间方差设置灰度阈值得到块内流路结构可能成像区域。
具体包括：
S51：对于整幅重建图Rec，仍沿用设置的分块边长L，再次进行分块处理。
与S2不同的，本次所分图像块平铺无交叠，由于下述步骤对于重建图像块做相同处理，故均以rec代指某一所分得的重建图像块。
S52：对于rec图像块像素进行基于最大类间方差(也称大津法，简称OTSU)的灰度阈值计算，阈值灰度记为T。
具体地，对于重建图块rec，记T为前景与背景的分割阈值，前景点数占图像比例为w0，平均灰度为u0；背景点数占图像比例为w1，平均灰度为u1，图像的总平均灰度为u，前景和背景图象的方差，则有：
u＝w0×u0+w1×u1
g＝w0×(u0-u)2+w1×(u1-u)2
联立上面两式可得：
g＝w0×w1×(u0-u1)2
或
当方差g最大时，认为此时前景和背景差异最大，灰度T即是最佳灰度阈值。
S53：根据所得到的灰度阈值T，对rec图像块进行二值化处理，得二值图像recBW，具体地：
S6：使用梯度信息并结合形态学滤波手段，对于S5各重建块中粗提取出的流路结构做进一步，排除掉被选中的细胞等观察用目标物体。
具体包括：
S61：使用Sobel算子对S5中提取出的recBW中流路结构可能区域所有边缘点(集合记为edge)坐标对应的重建图rec中的像素点进行梯度值计算。
对于坐标位置(x,y)的像素点，其横向梯度Gx与纵向梯度Gy计算方法如下：
Gx＝[2rec(x,y+1)-2rec(x,y-1)]+[rec(x-1,y+1)-rec(x-1,y-1)]+[rec(x+1,y+1)-rec(x+1,y-1)]
Gy＝[2rec(x+1,y)-2rec(x-1,y)]+[rec(x+1,y-1)+rec(x-1,y-1)]+[rec(x+1,y+1)-rec(x-1,y+1)]
该点对应梯度大小：
S62：计算recBW中所有边缘点对应rec像素点梯度的均值Gavr。
其中，num为rec中连通域边缘点总数。
S63：设recBW中共初检测出k个连通域，计算其各连通域边缘对应rec中像素灰度梯度均值Gt，t∈1,2,…,k。
对于Gt值低于Gavr的强度的连通域，其对应物体与背景对比度不足，排除为流路结构成像的可能，认为是由于光场分布不均或流路结构的共轭像所致。
S64：进行形态学滤波，根据实际观察细胞尺寸设置滤波范围，即实物面积对应在重建图上像素点个数n，经过1.5倍放大，去除所含像素点数少于1.5n的连通域。记此时二值矩阵为M。
S7：经过S6中梯度判断和形态学滤波处理后，仍保留下的区域(M中真值)认为是微流控芯片流路结构成像区域。对于rec区域对应的全息图H重建所得的复矩阵U，将流路结构对应区域的数值加以保留，而将其他位置处的数据实部、虚部分别替换为U矩阵各元素实部、虚部的均值。
S71：计算重建所得复矩阵U的实部R与虚部系数O的均值，分别记为real和imag。
S72：根据M矩阵逻辑值分布，构建针对流路成像区域的复数矩阵S。
S8：按照卷积重建步骤反向执行，反衍射出复数矩阵S对应的衍射信息矩阵holo。
S9：将重建图块rec对应原始全息图H中流路结构的衍射信息抵消，此时全息图中将仅(主要)包含细胞物体的衍射信息。再对抵消后的全息图重建，得到纯净的细胞物体重建图。
具体包括：
S91：计算反衍射矩阵holo元素均值holoavr。
S92：从全息图像块H中抵消掉衍射环矩阵holo所含信息，得到仅含有细胞衍射信息的全息图像块H*。
H*(x,y)＝H(x,y)-2[holo(x,y)-holoavr]
S93：对矩阵H*进行卷积重建，得重建复矩阵U*。与S4相同的，得到U*对应的实数矩阵I′，并对其进行归一化和取整处理。
S10：将各去除了芯片流路结构的重建图块按照对应全息图块分块顺序重新排列，得到芯片流路内所有细胞分布的完整重建图。
S0：特别的，对于S3、S8、S9中的全息图重建(逆)方法，采用卷积重建法加以实现，具体计算流程如下：
卷积重建法是基于线性系统理论和瑞利-索默菲积分公式而构建的，能够准确反映光的空间衍射分布。由瑞利-索默菲衍射积分公式，全息图中衍射信息再现后对应光波：
其中，
j在此表示虚数单位，λ为重建光波长。
H(x,y)为图像传感器所记录的全息图光强分布，也即是所采集记录得到的全息图。
C(x,y)是重建光波在传感器平面上的分布，取为R(x,y)的共轭：
R(x,y)＝exp[-jk(xsinα+ysinβ)]
式中α、β分别为平面参考光波与空间yoz平面和xoz平面的夹角。
而根据线性系统理论，再现像复振幅可以表示为如下卷积形式：
其中，k＝2π/λ，g为自由空间脉冲响应。
再由卷积理论，可使用傅里叶变换计算得到再现像，也即是本发明所使用的重建图：
Z＝FT-1[FT(H·C)·FT(g)]
最后说明的是，以上优选实施例仅用以说明发明的技术方案而非限制，尽管通过上述优选实施例已经对本发明进行了详细的描述，但本领域技术人员应当理解，可以在形式上和细节上对其作出各种各样的改变，而不偏离本发明权利要求书所限定的范围。