CN111951315B - 皮肤光学相干层析图像配准的方法 - Google Patents

Abstract

一种皮肤光学相干层析图像的配准方法，首先对采集的两组快速扫描方向正交的光学相干层析三维数据集进行预处理，再利用多层分辨率法与梯度优化法对两组数据集进行先纵向后整体的两步配准，消除运动伪影。完成配准后，将两组三维数据集进行加权融合，并沿纵向累加得到二维图像。本发明可以消除受运动影响严重的皮肤光学相干层析图像中的运动伪影，并通过图像融合提高了图像质量。

Description

皮肤光学相干层析图像配准的方法

技术领域

本发明涉及光学相干层析成像(Optical Coherence Tomography，简称OCT)，特别是一种皮肤光学相干层析图像配准的方法。

背景技术

光学相干层析成像(Optical Coherence Tomography，以下简称OCT)是一种非侵入、高分辨率、可在体检测生物组织内部微结构的生物医学光学成像技术。1991年，美国麻省理工学院的J.G.Fujimoto和D.Huang等人首先提出了此概念，并对视网膜和冠状动脉进行了离体成像。OCT可分为时域OCT(Time Domain OCT，TDOCT)和频域OCT(Fourier DomainOCT，FDOCT)。频域OCT相比时域OCT，成像速度更快且信噪比更高。目前OCT技术已被广泛应用于眼科、皮肤科、心血管等领域的临床诊断和研究。

光学相干层析微血管造影技术(Optical Coherence Tomography Angiography,OCTA)是OCT技术的功能化拓展，不但可以获得样品的结构差异，而且可以得到层析的血管图像。OCTA利用血红细胞等粒子运动引起的干涉信号变化，计算得到组织中的三维血液灌注图。相比于FA与ICGA，OCTA是一种可以提供更加完善的组织结构和血管结构信息的三维无标记、安全、快速的成像手段。OCTA最初主要用于视网膜微血管成像，随后在皮肤微血管成像方面也得到了应用。

OCTA数据采集时需要在同一位置连续采集具有一定时间间隔的B-scan，通过比较B-scan之间的去相关信号以重建微血管图像。数据采集过程中，不可避免地会受到样品运动、环境噪声、系统振动等影响，使数据不连续，并表现为血管图像的扭曲和断裂，在图像中产生运动伪影。样品运动造成的伪影会严重降低OCTA的成像质量，并影响对组织微血管的定量研究。样品运动有多种来源，如呼吸和心搏、被测部位不自主的抖动等。即使很小的运动，如心搏在皮肤中造成的微米级运动，都会增加背景噪声，最终影响血管成像效果。数据采集的时间间隔越长，运动伪影越严重。因此在OCTA的血管提取算法中，运动伪影的消除是一项关键内容。伪影消除的主要方法是图像配准，即利用变换、插值、位移等方法将图像中受运动影响产生的扭曲、断裂等校正，得到无运动伪影的图像。

现有的OCT伪影消除方法主要基于硬件方法和软件方法。

硬件方法中，一种方式是使用高速OCT以减少采集时间进而减少采样中的运动(参见在先技术[1]Z.Zhi et al.,“4D optical coherence tomography-based micro-angiography achieved by 1.6-MHz FDML swept source,”Opt.Lett.40(8),1779–1782(2015).)；另一种方式是使用外加的仪器探测采集中的运动以使系统重新扫描以消除运动影响(参见在先技术[2]B.Braaf et al.,“Real-time eye motion correction in phase-resolved OCT angiography with tracking SLO,”Biomed.Opt.Express 4(1),51–65(2013).)。但硬件方法中，使用高速OCT的方式增加了系统成本，使用外加仪器的方式增加了系统复杂度和成本。

软件方法中，一种方式是利用互相关算法对连续B-scan进行处理实现图像三维配准(参见在先技术[3]Zawadzki R J,Fuller AR,Choi S S,et al.Correction of motionartifacts and scanning beam distortions in 3D ophthalmic optical coherencetomography imaging[C]//Ophthalmic Technologies XVII.International Society forOptics and Photonics,2007,6426:642607.)；另一种方式是通过表面分割，从OCT结构图中提取了组织特征，并进行几何变换以校正横向扭曲(参见在先技术[4]Antony B,Abramoff M D,Tang L,et al.Automated 3-D method for the correction of axialartifacts in spectral-domain optical coherence tomography images[J].Biomedical optics express,2011,2(8):2403-2416.)。但这几种方式均难以满足微血管造影的配准精度需求。

发明内容

本发明的目的是为了克服上述在先技术的不足，提供一种三维皮肤OCT图像配准方法。对采集的两组快速扫描方向正交的OCT三维数据集进行预处理，再利用多层分辨率法与梯度优化法对两组数据集进行先纵向后整体的两步配准，消除运动伪影。完成配准后，将两组三维数据集进行加权融合，并沿纵向累加得到二维图像。本发明可以消除受运动影响严重的皮肤光学相干层析图像中的运动伪影，并通过图像融合提高了图像质量。

本发明的技术解决方案如下：

一种皮肤光学相干层析图像配准的方法，其特点在于该方法包括以下步骤：

①定义一个扫描坐标系，该坐标系为空间三维直角坐标系，x和y为水平维度，z为竖直维度；

②定义光学相干层析成像使用栅线扫描(raster scan)方式采集三维数据集时，采集B-scan的方向称为快速扫描方向；

③使用光学相干层析成像系统在待测皮肤的成像区域使用栅线扫描方式采集一组三维数据集命名为XRAW；XRAW中所有体素定义在①中所述扫描坐标系中坐标(x_i,y_j,z_k)上，x_i,y_j,z_k分别表示x、y、z轴上的第i、j、k个正整数点，x_i＝1,2,3…w_x，y_j＝1,2,3…h_x，z_k＝1,2,3…d_x，w_x、h_x、d_x分别为XRAW在x、y、z三轴上的采样点数；快速扫描方向沿x轴方向；

④使用③中同一光学相干层析成像系统在待测皮肤与③中同一成像区域使用栅线扫描方式采集另一组三维数据集命名为YRAW；YRAW中所有体素定义在①中所述扫描坐标系中坐标(x_i,y_j,z_k)上，x_i,y_j,z_k分别表示x、y、z轴上的第i、j、k个正整数点，x_i＝1,2,3…h_y，y_j＝1,2,3…w_y，z_k＝1,2,3…d_y，h_y、w_y、d_y分别为YRAW在x、y、z三轴上的采样点数；快速扫描方向沿y轴方向；h_y＝w_x，w_y＝h_x，d_x＝d_y；

⑤对三维数据集XRAW和YRAW做预处理，得到三维数据集XFAST和YFAST，具体步骤如下：

5.1)对XRAW和YRAW分别沿z方向做一维中值滤波，得到两组三维数据集XP1和YP1；

5.2)对XP1和YP1分别做对数变换，得到XP2和YP2；

5.3)对XP2和YP2做阈值处理：计算XP2和YP2中所有体素强度的直方图，取直方图中最大值对应的区间的中值作为阈值下限thmin，取所有体素强度的最大值作为阈值上限thmax，对XP2和YP2分别进行阈值处理和归一化得到两组三维数据集，命名为XP3，YP3：

5.4)对XP3和YP3分别进行照明校正：对XP3和YP3分别做沿z方向的累加，得到二维数组Fx和Fy；选择Fx中大于75％体素的强度的强度值作为v_ref,x，选择Fy中大于75％体素的强度的强度值作为v_ref,y，用下式计算照明校正强度二维图Bx和By：

Bx_i,j＝v_ref,x-Fx_i,j

By_i,j＝v_ref,y-Fy_i,j

定义阈值v_thx和v_thy，认为XP3中强度大于v_thx的体素为信号体素，YP3中强度大于v_thy的体素为信号体素；定义三维数据集Mx和My：

对XP3和YP3做下式处理，得到结果XP4和YP4：

5.5)对XP4和YP4分别做一次沿z方向系数为2的下采样；将下采样后的数据集化为均值为0，方差为1的数据集，结果分别命名为XFAST和YFAST。

⑥定义XFAST为0级分辨率数据集，对XFAST做N次沿x、y、z方向的系数为2的下采样，下采样一次的结果为1级分辨率数据集，命名为XDOWN₁，以此类推，下采样N次的结果为N级分辨率数据集，命名为XDOWN_N；

定义YFAST为0级分辨率数据集，对YFAST做N次沿x、y、z方向的系数为2的下采样，下采样一次的结果为1级分辨率数据集，命名为YDOWN₁，以此类推，下采样N次的结果为N级分辨率数据集，命名为YDOWN_N；

⑦定义采集XRAW时，在y_j＝j处采集的B-scan的采样时刻为tx_j，其中y_j＝1处采集的B-scan的采样时刻tx₁＝0；定义采集YRAW时在x_i＝i处采集的B-scan的采样时刻为ty_i；ty_i大于tx_j；

⑧定义在采样时刻tx₁＝0时，按③中所述的采样点采得的三维数据集为理想采样数据集，理想采样数据集的B-scan称为理想采样B-scan；从y_j＝j(j＞1)处的理想采样B-scan运动到tx_j(j＞1)时刻的实际采样的B-scan在z方向的运动量为D_z(tx_j)，从x_i＝i处的理想采样B-scan运动到ty_i时刻的实际采样点在z方向的运动量为D_z(ty_i)；

⑨对两组N级分辨率数据集XDOWN_N、YDOWN_N进行纵向配准，得到N级分辨率数据集的纵向运动量，将N级分辨率数据集的纵向运动量进行插值得到N-1级纵向运动量初值，对N-1级分辨率数据集XDOWN_N-1、YDOWN_N-1进行纵向配准，以此类推，直到对XFAST、YFAST纵向配准，得到0级分辨率数据集的纵向运动量，利用0级分辨率数据集的纵向运动量D_z(tx_j)和D_z(ty_i)，使用插值方法对XFAST和YFAST进行变换：变换前XFAST_i,j,k＝V(XFAST,x_i,y_j,z_k)，表示三维数据集XFAST在第i列、第j行、第k页的值为XFAST在三维坐标系中(x_i,y_j,z_k)点的值，变换后得到新的三维数据集XCOR，XCOR_i,j,k＝V(XFAST,x_i,y_j,z_k+D_z(tx_j))，变换前YFAST_i,j,k＝V(YFAST,x_i,y_j,z_k)，表示三维数据集YFAST在第i列、第j行、第k页的值为YFAST在三维坐标系中(x_i,y_j,z_k)点的值，变换后得到新的三维数据集YCOR，YCOR_i,j,k＝V(YFAST,x_i,y_j,z_k+D_z(ty_i))；

所述的纵向配准的具体步骤如下：

9.1)将进行纵向配准的两组三维数据集分别命名为X₀和Y₀；定义数据集X₀在y_j＝j处的纵向运动量D_z(tx_j)由平移运动量和旋转运动量组成；平移运动量为从y_j＝j(j＞1)处的理想采样B-scan运动到tx_j(j＞1)时刻的X₀的B-scan在z方向上的平移运动量为D_zd(tx_j)；旋转运动量为从y_j＝j(j＞1)处的理想采样B-scan运动到tx_j(j＞1)时刻的X₀的B-scan在B-scan平面内的旋转运动量D_zt(tx_j)；D_z(tx_j)＝D_zd(tx_j)+D_zt(tx_j)；定义数据集Y₀在x_i＝i处的纵向运动量D_z(ty_i)由平移运动量和旋转运动量组成；平移运动量为从x_i＝i处的理想采样B-scan运动到ty_i时刻的Y₀的B-scan在z方向上平移运动量D_zd(ty_i)；旋转运动量为从x_i＝i处的理想采样B-scan运动到ty_i时刻的Y₀的B-scan在B-scan平面内的旋转运动量D_zt(ty_i)；D_z(ty_i)＝D_zd(ty_i)+D_zt(ty_i)；

9.2)使用插值方法对X₀变换，变换前X_0i,j,k＝V(X₀,x_i,y_j,z_k)，变换后得到新的三维数据集X_g，X_gi,j,k＝V(X₀,x_i,y_j,z_k+D_z(tx_j))；使用插值方法对Y₀变换，变换前为Y_0i,j,k＝V(Y₀,x_i,y_j,z_k)，变换后得到新的三维数据集Y_g，Y_gi,j,k＝V(Y₀,x_i,y_j,z_k+D_z(ty_i))；

9.3)将变换后的数据集做差，得到残差：

R_i,j,k(X_g,Y_g,D_z)＝X_{g i,j,k}-Y_{g i,j,k}

将一种罚函数L(x)应用于残差，定义相似度项S(D_z)：

定义正则项E(D_z)：

其中tx_l和ty_l分别为采集X₀和Y₀时第l个采样时刻；

将相似度项和正则项组合，得到目标函数：

F(D_z)＝S(D_z)+αE(D_z)，

其中α为正则项相对相似度项重要程度的权值，计算目标函数对于D_z的一阶导数G(D_z)；

9.4)使用迭代优化算法，利用G(D_z)得到纵向运动量D_z，公式如下：

D_z＝argminF(D_z)

其中argmin为求使F(D_z)最小的D_z。

⑩定义XCOR为0级分辨率数据集，对XCOR做N次沿x、y、z方向的系数为2的下采样，下采样一次的结果为1级分辨率数据集，命名为XCORDOWN₁，以此类推，下采样N次的结果为N级分辨率数据集，命名为XCORDOWN_N；

定义YCOR为0级分辨率数据集，对YCOR做N次沿x、y、z方向的系数为2的下采样，下采样一次的结果为1级分辨率数据集，命名为YCORDOWN₁，以此类推，N次的结果为N级分辨率数据集，命名为YCORDOWN_N；

根据⑧所述，从y_j＝j(j＞1)处的理想采样B-scan运动到tx_j(j＞1)时刻的实际采样的B-scan在x、y、z方向的运动量分别为D_x(tx_j)、D_y(tx_j)、D_z(tx_j)，从x_i＝i处的理想采样B-scan运动到ty_i时刻的实际采样点在x、y、z方向的运动量分别为D_x(ty_i)、D_y(ty_i)、D_z(ty_i)；

对两组N级分辨率数据集XCORDOWN_N、YCORDOWN_N进行整体配准，得到N级分辨率数据集的整体运动量，将N级分辨率数据集的整体运动量进行插值得到N-1级整体运动量初值，对N-1级分辨率数据集XCORDOWN_N-1、XCORDOWN_N-1进行整体配准，以此类推，直到对XCOR、YCOR整体配准，得到0级分辨率数据集的整体运动量，利用0级整体运动量D_x(tx_j)，D_y(tx_j)，D_z(tx_j)和D_x(ty_i)，D_y(ty_i)，D_z(ty_i)，使用插值方法对XCOR和YCOR进行变换：变换前XCOR_i,j,k＝V(XCOR,x_i,y_j,z_k)，表示三维数据集XCOR在第i列、第j行、第k页的值为XCOR在三维坐标系中(x_i,y_j,z_k)点的值，变换后得到新的三维数据集XREG，XREG_i,j,k＝V(XCOR,x_i+D_x(tx_j),y_j+D_y(tx_j),z_k+D_z(tx_j))，变换前YCOR_i,j,k＝V(YCOR,x_i,y_j,z_k)，表示三维数据集YCOR在第i列、第j行、第k页的值为YCOR在三维坐标系中(x_i,y_j,z_k)点的值，变换后得到新的三维数据集YREG，YREG_i,j,k＝V(YCOR,x_i+D_x(ty_i),y_j+D_y(ty_i),z_k+D_z(ty_i))；

所述的整体配准具体步骤如下：

12.1)将进行整体配准的两组三维数据集分别命名为X₀和Y₀；定义数据集X₀在y_j＝j处在x、y、z三个方向的运动量分别为D_x(tx_j)，D_y(tx_j)，D_z(tx_j)，运动量为从y_j＝j(j＞1)处的理想采样B-scan运动到tx_j(j＞1)时刻X₀中B-scan在x、y、z方向上的平移运动量；定义数据集Y₀在x_i＝i处x、y、z三个方向的运动量分别为D_x(ty_i)，D_y(ty_i)，D_z(ty_i)，运动量为从x_i＝i处理想采样B-scan运动到ty_i时刻Y₀中B-scan在x、y、z方向上的平移运动量；

12.2)使用插值方法对X₀变换，变换前X_0i,j,k＝V(X₀,x_i,y_j,z_k)，变换后得到新的三维数据集X_g，为X_{g i,j,k}＝V(X₀,x_i+D_x(tx_j),y_j+D_y(tx_j),z_k+D_z(tx_j))；使用插值方法对Y₀变换，变换前为Y_0i,j,k＝V(Y₀,x_i,y_j,z_k)，变换后得到新的三维数据集Y_g，Y_{g i,j,k}＝V(Y₀,x_i+D_x(ty_i),y_j+D_y(ty_i),z_k+D_z(ty_i))；

12.3)将变换后的数据集做差，得到残差：

R_i,j,k(X_g,Y_g,D_x,D_y,D_z)＝X_{g i,j,k}-Y_{g i,j,k}

将一种罚函数L(x)应用于残差，定义相似度项S(D_x,D_y,D_z)：

定义正则项E(D_z)：

其中tx_l和ty_l分别为采集X₀和Y₀时第l个采样时刻；

将相似度项和正则项组合，得到目标函数：

F(D_x,D_y,D_z)＝S(D_x,D_y,D_z)+αE(D_x,D_y,D_z)

其中α为正则项相对相似度项重要程度的权值，计算目标函数对于D_x,D_y,D_z的一阶导数G(D_x,D_y,D_z)；

12.4)使用迭代优化算法，利用G(D_x,D_y,D_z)得到整体运动量D_x,D_y,D_z：

D_x,D_y,D_z＝argminF(D_x,D_y,D_z)

其中argmin为求使F(D_x,D_y,D_z)最小的D_x,D_y,D_z。

将两组三维数据集XREG和YREG加权融合，得到V_merged，将V_merged沿z向累加，得到En-face投影图，具体步骤如下：

计算y_j+D_y(tx_j)的概率密度分布和x_i+D_x(ty_i)的概率密度分布，XREG的y_j＝j处B-scan中的A-line的权值为y_j+D_y(tx_j)的概率密度的倒数；YREG的x_i＝i处B-scan中的A-line的权值为对应的x_i+D_x(ty_i)的概率密度的倒数，XREG和YREG的每个对应A-line的权值做归一化，得到融合时的权值。

本发明与现有技术相比有益的效果是：

1.与在先技术[1]相比，本发明不需要高采样速度的OCT系统。

2.与在先技术[2]相比，本发明不需要额外的运动探测仪器。

3.与在先技术[3]和[4]相比，本发明可以将受运动影响严重，含有大量运动伪影的数据集中的运动校正，实现配准。

4.本发明可以消除受运动影响严重的皮肤光学相干层析图像中的运动伪影，并通过图像融合提高了图像质量。

具体实施方式

下面结合实施例对本发明作进一步说明，但不应以此实施例限制本发明的保护范围。

1.配准原理

数据采集时，扫描由二维振镜按栅线扫描方式实现。扫描过程中，二维振镜两个轴的扫描方向分别称为快速扫描方向和慢速扫描方向。沿快速扫描方向完成一次连续扫描，得到一帧B-scan。每完成一次快速扫描，快速扫描振镜转回起始位置，慢速扫描振镜前进一步，直至OCT系统采集一组完整的三维数据。对同一采样区域，连续采集两组三维数据，两组数据的快速扫描方向正交。将得到的两组三维数据集分别命名为XFAST和YFAST。在以振镜扫描的x、y维度所组成的扫描坐标系中，两次栅线扫描的对应采样点位置重合，即所有采样点均进行了两次采样。

两个数据集各自定义在一个三维直角坐标系中。坐标系中的点定义为X(i,j,k)和Y(i,j,k)，i、j、k分别为x、y、z三个维度中的索引。每个数据集在水平方向上有h×w个A-scan，每个A-scan在竖直方向由d个体素组成。

和/>是XFAST和YFAST中每一个A-scan的采样时刻。使用下面的函数定义一个数据集中的每一个体素值

V＝A(x,y,z,t) (1)

其中x、y、z是扫描坐标系的三个坐标，t是采样时刻。每个水平方向坐标(x,y)处，沿竖直方向z＝1,2…d的体素构成了一个A-scan。每个体素的值除了与采样坐标相关，也与采样时刻相关。因为当OCT与被测组织之间发生相对运动时，体素值也会随之改变。基于以上定义，数据集X和Y可定义为：

x_i,Y_j,z_k是扫描坐标系中x、y、z方向第i、j、k个正整数坐标点。

将开始采集第一个数据集的第一个A-line的时刻被测皮肤与系统的对准位置作为理想采样位置，在此对准下对待测皮肤采样得到的体素为理想采样点。由于数据采集时的样品运动或环境扰动，理想采样点并不在实际采样坐标(x_i,y_j,z_k)。设理想采样坐标为(x_pi,y_pj,z_pk)，x_pi,y_pj,z_pk分别表示在x、y、z三个维度上与实际采样坐标(x_i,y_j,z_k)相对应的理想采样坐标，这些理想采样坐标不完全位于正整数坐标点上。运动是随时间变化的，因此某一时刻t，在物理坐标位置(x_i,y_j,z_k)处采得体素值为：

其中D_x(t)、D_y(t)、D_z(t)分别是t时刻采集的A-scan中的体素相对理想采样坐标在x、y、z三个方向上的运动偏移量。

若要得到在理想采样点(x_pi,y_pj,z_pk)的体素值，需使用运动偏移量在V中进行插值，定义插值函数I(V,x_i,y_j,z_k)，则：

V(x_pi,y_pj,z_pk)＝I(V,x_i,y_j,z_k)＝V(x_i+D_x(t),y_j+D_y(t),z_k+D_z(t)) (5)

只有求得D_x、D_y、D_z才能得到理想采样点的体素值。由于XFAST和YFAST是对同一位置进行采样的数据集，因此在XFAST和YFAST中对同一理想坐标(x_pi,y_pj,z_pk)处进行插值得到的体素值应高度相似。定义两者之间的残差R：

当两者相似度高时残差项最小。定义XFAST和YFAST两数据集之间的相似度项S表示两者插值后的相似程度：

S(X,Y,D_x,D_y,D_z)＝∑_i∑_j∑_kL(R_i,j,,k(X,Y,D_x,D_y,D_z)) (7)

其中L是损失函数。当插值后的两数据集之间相似度最高时，损失函数的值应最小，即：

D_x,D_y,D_z＝argminS(X,Y,D_x,D_y,D_z) (8)

其中argmin表示使S(X,Y,D_x,D_y,D_z)最小的D_x,D_y,D_z。在相似度之外，还应对运动的变化速度进行惩罚。在实际中，运动通常为平缓运动，在相对小的时间间隔内不会发生大的运动量改变，因此定义正则项：

其中L(x)也是损失函数。取出运动量对时间间隔的导数即运动速度，对其应用损失函数，并将所有时间间隔对应的速度累加，得到正则项。当求得的运动速度越低时，该项越小。

将式(7)和式(8)组合，得到基本的优化目标函数：

D_x,D_y,D_z＝argmin(S(D_x,D_y,D_z)+αE(D_x,D_y,D_z)) (10)

其中α为权值，表示相似度项和正则项之间的相对重要程度。得到目标函数后，利用优化算法，通过对目标函数的参数进行优化，使目标函数值最小，以求解D_x,D_y,D_z，实现运动偏移量的校正。

2.图像预处理

为了提高配准算法处理步骤的运行速度和准确度，需要对采集的数据进行预处理，以减小两组数据之间的不一致，降低噪声，并且提高配准算法的运行速度。图像预处理主要包括去噪及阈值化、照明校正、下采样和标准化。

2.1去噪及阈值化

为了降低散斑噪声，对两组数据集进行空域中值滤波。对每一个A-scan沿深度方向进行一维中值滤波，滤波窗格根据数据集的信噪比进行选择。

中值滤波后，对两组数据进行对数化处理，以降低原始数据的动态范围，再通过阈值化处理消除背景噪声并归一化。由于在采集的数据中，背景体素数量最多，因此背景体素中占比最高的灰度与背景体素灰度平均值近似，而灰度低于背景体素灰度平均值的体素不包含有效信号。因此利用灰度直方图对两组数据集进行阈值化处理。首先求得两组数据集中所有体素灰度的直方图，直方图中包含128个组。将体素分布最多的灰度设为最低阈值thmin，两组数据集中最大灰度值设为最高阈值thmax，对每一数据集应用阈值处理：

通过阈值化处理消除了背景噪声，同时将所有体素归一化到0和1之间。

2.2照明校正

阈值化处理后，对两个数据集进行照明校正。在数据采集的过程中，受样品运动影响，OCT系统与被测组织的对准会发生改变，导致对同一采样位置的照明不同，引起干涉信号强度的变化，进而影响图像配准。为补偿照明的变化量，首先计算每个数据集的En-face投影图，将其定义为：

假设照明的变化是缓慢的，因此照明改变是低频的。对投影图进行二维高斯低通滤波，标准差设为7.5像素。将滤波结果二维图记为设用于校正照明的二维图为C_i,j，参考值为v_ref，得到照明校正图：

参考值v_ref为二维图中灰度分布占75％的值，这使得照明校正图中75％的像素为正。将照明图应用于对应的数据集V中每个A-scan有效信号的体素上。定义一个掩模M：

其中v_ret是区分信号与背景的阈值，认为大于该值的体素是信号。将照明场应用于数据集：

4.3下采样和标准化

将经过照明校正的数据集沿深度方向进行下采样，以减少数据量。先沿深度方向做高斯滤波，然后去除其中的偶数项。深度方向的下采样可以进一步降低噪声，但也会降低纵向分辨率，影响纵向校正精度。

最后，对经过以上处理的数据集进行标准化处理，使两个数据集各自的均值为0，方差为1。

3.图像配准算法

完成预处理后，进行两组数据集的图像配准。通过寻找适当参数最小化目标函数，计算每个数据集采样时被测组织的运动量，再利用插值实现对图像的配准。基本目标函数如1节所述，为相似度项和正则项的加权和。

3.1相似度测量

尽管在预处理中进行了照明校正和去噪，但由于数据集中的散斑噪声为非高斯分布，且椒盐噪声难以完全去除，照明不均也只是得到一定程度的补偿，数据中仍然存在剩余噪声。这可能导致两组数据集中同一物理位置的体素强度不等，从而引起相似度测量不准确。为解决以上问题，使用基于伪胡贝尔范数作为相似度项的损失函数。伪胡贝尔范数是L1范数的拓展。L1范数如下式：

L1(x)＝|x| (16)

使用L1范数作为损失函数，可以较好地解决数据不一致和噪声等导致的残差较大问题，具有更高的鲁棒性。然而，由于L1范数非连续可导，若使用其作为损失函数，则目标函数同样非连续可导，不适用于基于二阶梯度的非线性优化。因此，此处使用伪胡贝尔范数：

其中ε_H为一个很小的正常数，该值越接近0，式(4.17)越近似于L1范数，当x远大于ε_H，斜率近似为1，表示损失函数值随x线性变化。ε_H通常选为0.00001。

3.2正则化

除了相似度项，基本的目标函数中的另一项是正则项。相似度项的目的是使优化结果中两组数据集相似，正则项是为了保证得到的运动量与实际相符。根据实际情况，被测样品在相对短时间内不会发生大幅度的运动，即运动速度较小。因此，运动量对时间的导数越大，正则项的损失函数值越大。

通过计算在每个采样时刻的运动量对时间的导数以得到正则项。在实际中，使用有限差分法求得运动量对时间的导数。以D_x为例，使用前向差分，其在t_l采样时刻的运动量对时间的导数为：

其中t_l+1是下一个采样的时刻。同理，可得到D_y、D_z对时间的导数。

正则项为所有时刻的运动量对时间间隔的损失函数的和，即：

在实际情况中，通常在较长采样时间内运动变化缓慢，在短时间内出现运动速度快而幅度大的抖动。因此使用L0.5范数作为运动损失函数具有更好的适应性。在一维运动中：

L0.5范数可以成比例地惩罚缓慢的运动，也会较少地惩罚大幅的运动。但是L0.5范数并非连续可导，也不符合目标函数需要二阶连续可导的条件。因此使用L0.5范数的近似表示：

其中ε_0.5是一个正的常量，决定了上式与L0.5范数的近似程度。ε_0.越小，与L0.5范数越接近。但是ε_0.5越接近0，损失函数的导数越近似于x＝0处的导数，因此选择ε_0.5＝0.1。

3.3运动量参数化

在每个采样时刻分配x、y、z三个维度的运动量参数，表示在该采样时刻时，样品在三个维度上的运动量。使用优化算法对这些参数进行优化，以使前述目标函数的值最小。

3.4数据集插值

通过使用插值函数I得到t时刻理想物理采样位置的A-scan。插值结果由t时刻真实采得的A-scan和运动量(D_x(t),D_y(t),D_z(t))决定。使用三维三阶样条法进行插值。使用差分求导法得到每个插值点的一阶导数。当被插值点位于原数据集定义的坐标网格之外时，使用数据集内距离被插值点最近的点作为插值点。

3.5优化算法

由于目标函数的非线性和非凸性，使用迭代算法对其进行求解。目标函数中含有大量参数，而插值时需计算每个参数，因此计算复杂度较高。为降低计算复杂度，使用limited-memory Broyden-Pletcher-Goldfarb-Shanno(L-BFGS)方法对目标函数进行优化。L-BFGS方法是一种基于梯度的拟牛顿优化方法，在优化时需要将目标函数每一项的梯度数学表达式也作为该方法的输入。相较于牛顿优化方法，L-BFGS算法计算量更少，因此可以更快地进行优化，适用于含有大量参数的优化问题。

3.6多级分辨率优化

由于目标函数的非凸性，若对所有参数进行迭代数值优化只能得到目标函数的局部极小值。因此，使用多级分辨率优化以得到全局最小值。

多级分辨率优化将一个复杂的优化问题化为一系列相对简单的优化问题。将简单问题的结果作为下一级更复杂问题的初始值。由于简单问题已经得到了接近最优解的值，因此该值可以作为下一级更复杂问题的良好初始值。从最简单的问题开始优化，逐步增加问题复杂度，直到将原始问题优化。

通过对原始三维数据集进行下采样，得到一系列分辨率更低的三维数据集，从而建立数据集金字塔。通过对数据集高斯滤波后取每个维度的奇数体素，对原始数据集进行连续下采样。将所得的每一个下采样数据集进行标准化，使其均值为0，方差为1。

对时间模型使用相似的方法进行下采样，以与数据集的数据相对应。将低级下采样数据集优化的结果转化为更高级的数据集的初始值时，使用双线性插值以匹配更高分辨率的数据集。

3.7分步优化

同一B-scan内的A-scan之间的采样时间间隔较短，因此认为同一B-scan内的A-scan之间无相对运动，将B-scan作为刚体进行运动校正。为了提高鲁棒性，使用了多分辨率与多步优化结合的两步优化方法。第一步为粗配准，只校正纵向运动。对每个B-scan的纵向维度进行优化，将得到的纵向运动量应用于原数据集上，校正数据集中的纵向运动，将其作为第二步全局配准的初始数据集。在第二步中，对每个B-scan所有维度的运动量进行校正，实现精细配准。在两步配准过程中，均使用多级分辨率方法，逐级优化以得到使目标函数值为全局最小值的三个维度的运动量。

第一步配准基本将被测部位的纵向运动校正，第二步配准时在初始值基本无纵向运动的情况下，在所有维度上进行精细配准。如前所述，迭代优化方法一般仅能得到最近的局部最优解，为减轻该问题，使用了多步配准方法。当第一步配准已将纵向运动量基本校正后，第二步配准的初始数据集就更加接近全局最优配准结果，减少同时对纵向和横向维度配准时仅得到局部最优解的情况。

除了仅沿纵向的相对运动之外，还存在样品光入射角度相对被测部位发生扭转导致的运动伪影。因此在第一步对准中，不仅需要校正纵向运动，也需要对扭转进行校正。在纵向运动量参数中加入一个校正每个B-scan扭转的扭转参量，通过将该值与B-scan中的每一个A-scan的索引相乘，表示扭转引起的运动量。同样，因为扭转与仅沿一个维度的运动类似，不会在短时间内大幅度改变，所以在正则化项中也要惩罚其相对时间的变化量。

因此在第一步对准中，相似度项中的残差在(x_i,Y_j,z_k)处的表达式为：

其中D_zt为在采样时刻(x_i,y_j,z_k)的扭曲率。

4.图像融合

完成精细配准后，将配准后的两组数据集进行融合，以得到信噪比提高的数据集。假设两组数据集中对应位置均被配准，且两组数据集中的散斑噪声不相关，则图像融合会提高信噪比。

运动可能导致被测部位的某个位置在一个数据集中未被采集，而在另一数据集中则被采集到。在这种情况下，尽管原数据集中不包含该位置的数据，但是配准过程中会将其相邻部分的相似数据插值至此处以使两数据集相似。这样导致某实际采样的A-scan被重复多次插值于其他位置。若原数据集中某一位置的数据被多次插值于其他位置，则其对应的配准后的理想采样点可能为无效点。图像融合时应该降低这些点的权重。

使用帕尔森密度估计得到原数据集中每个B-scan被用于插值的采样密度。由于配准中将B-scan视为刚体进行优化，若某处B-scan被多次插值到慢速扫描方向上的其他位置，说明该B-scan被用来补充未采到的部分数据。统计用来插值的原始B-scan慢速扫描坐标，使用窗长度为0.05的帕尔森窗核函数对各数据集中慢速扫描方向上被用于插值的B-scan进行采样密度估计。

得到采样密度后，将两个数据集中对应B-scan的采样密度比值的反比作为权重，将每组权重归一化后分配给两组数据集。将两组数据集加权求和，得到融合结果。

Claims (5)

1.一种皮肤光学相干层析图像的配准方法，其特征在于，该方法包括以下步骤：

①定义一个扫描坐标系，该扫描坐标系为空间三维直角坐标系，x和y为水平维度，z为竖直维度；

②定义光学相干层析成像使用栅线扫描方式采集三维数据集时，采集B-scan的方向称为快速扫描方向；

③使用光学相干层析成像系统在待测皮肤的成像区域使用栅线扫描方式采集一组三维数据集命名为XRAW；XRAW中所有体素定义在①中所述扫描坐标系中坐标(x_i，y_j，z_k)上，x_i，y_j，z_k分别表示x、y、z轴上的第i、j、k个正整数点，x_i＝1，2，3...w_x，y_j＝1，2，3...h_x，z_k＝1，2，3…d_x，w_x、h_x、d_x分别为XRAW在x、y、z三轴上的采样点数；快速扫描方向沿x轴方向；

④使用③中同一光学相干层析成像系统在待测皮肤与③中同一成像区域使用栅线扫描方式采集另一组三维数据集命名为YRAW；YRAW中所有体素定义在①中所述扫描坐标系中坐标(x_i，y_j，z_k)上，x_i，y_j，z_k分别表示x、y、z轴上的第i、j、k个正整数点，x_i＝1，2，3...h_y，y_j＝1，2，3...w_y，z_k＝1，2，3…d_y，h_y、w_y、d_y分别为YRAW在x、y、z三轴上的采样点数；快速扫描方向沿y轴方向；h_y＝w_x，w_y＝h_x，d_x＝d_y；

⑤对三维数据集XRAW和YRAW做预处理，得到三维数据集XFAST和YFAST；

⑥定义XFAST为0级分辨率数据集，对XFAST做N次沿x、y、z方向的系数为2的下采样，下采样一次的结果为1级分辨率数据集，命名为XDOWN₁，以此类推，下采样N次的结果为N级分辨率数据集，命名为XDOWN_N；

定义YFAST为0级分辨率数据集，对YFAST做N次沿x、y、z方向的系数为2的下采样，下采样一次的结果为1级分辨率数据集，命名为YDOWN₁，以此类推，下采样N次的结果为N级分辨率数据集，命名为YDOWN_N；

⑦定义采集XRAW时，在y_j＝j处采集的B-scan的采样时刻为tx_j，其中y_j＝1处采集的B-scan的采样时刻tx₁＝0；定义采集YRAW时在x_i＝i处采集的B-scan的采样时刻为ty_i，且ty_i＞tx_j；

⑧定义在采样时刻tx₁＝0时，按③中所述的采样点采得的三维数据集为理想采样数据集，理想采样数据集的B-scan称为理想采样B-scan；从y_j＝j(j＞1)处的理想采样B-scan运动到tx_j(j＞1)时刻的实际采样的B-scan在z方向的运动量为D_z(tx_j)，从x_i＝i处的理想采样B-scan运动到ty_i时刻的实际采样点在z方向的运动量为D_z(ty_i)；

⑨对两组N级分辨率数据集XDOWN_N、YDOWN_N进行纵向配准，得到N级分辨率数据集的纵向运动量，将N级分辨率数据集的纵向运动量进行插值得到N-1级纵向运动量初值，对N-1级分辨率数据集XDOWN_N-1、YDOWN_N-1进行纵向配准，以此类推，直到对XFAST、YFAST纵向配准，得到0级分辨率数据集的纵向运动量，利用0级分辨率数据集的纵向运动量D_z(tx_j)和D_z(ty_i)，使用插值方法对XFAST和YFAST进行变换：变换前XFAST_i，j，k＝V(XFAST，x_i，y_j，z_k)，表示三维数据集XFAST在第i列、第j行、第k页的值为XFAST在三维坐标系中(x_i，y_j，z_k)点的值，变换后得到新的三维数据集XCOR，XCOR_i，j，k＝V(XFAST，x_i，y_j，z_k+D_z(tx_j))，变换前YFAST_i，j，k＝V(YFAST，x_i，y_j，z_k)，表示三维数据集YFAST在第i列、第j行、第k页的值为YFAST在三维坐标系中(x_i，y_j，z_k)点的值，变换后得到新的三维数据集YCOR，YCOR_i，j，k＝V(YFAST，x_i，y_j，z_k+D_z(ty_i))；

⑩定义XCOR为0级分辨率数据集，对XCOR做N次沿x、y、z方向的系数为2的下采样，下采样一次的结果为1级分辨率数据集，命名为XCORDOWN₁，以此类推，下采样N次的结果为N级分辨率数据集，命名为XCORDOWN_N；

定义YCOR为0级分辨率数据集，对YCOR做N次沿x、y、z方向的系数为2的下采样，下采样一次的结果为1级分辨率数据集，命名为YCORDOWN₁，以此类推，N次的结果为N级分辨率数据集，命名为YCORDOWN_N；

根据⑧所述，从y_j＝j(j＞1)处的理想采样B-scan运动到tx_j(j＞1)时刻的实际采样的B-scan在x、y、z方向的运动量分别为D_x(tx_j)、D_y(tx_j)、D_z(tx_j)，从x_i＝i处的理想采样B-scan运动到ty_i时刻的实际采样点在x、y、z方向的运动量分别为D_x(ty_i)、D_y(ty_i)、D_z(ty_i)；

对两组N级分辨率数据集XCORDOWN_N、YCORDOWN_N进行整体配准，得到N级分辨率数据集的整体运动量，将N级分辨率数据集的整体运动量进行插值得到N-1级整体运动量初值，对N-1级分辨率数据集XCORDOWN_N-1、XCORDOWN_N-1进行整体配准，以此类推，直到对XCOR、YCOR整体配准，得到0级分辨率数据集的整体运动量，利用0级整体运动量D_x(tx_j)，D_y(tx_j)，D_z(tx_j)和D_x(ty_i)，D_y(ty_i)，D_z(ty_i)，使用插值方法对XCOR和YCOR进行变换：变换前XCOR_i，j，k＝V(XCOR，x_i，y_j，z_k)，表示三维数据集XCOR在第i列、第j行、第k页的值为XCOR在三维坐标系中(x_i，y_j，z_k)点的值，变换后得到新的三维数据集XREG，XREG_i，j，k＝V(XCOR，x_i+D_x(tx_j)，y_j+D_y(tx_j)，z_k+D_z(tx_j))，变换前YCOR_i，j，k＝V(YCOR，x_i，y_j，z_k)，表示三维数据集YCOR在第i列、第j行、第k页的值为YCOR在三维坐标系中(x_i，y_j，z_k)点的值，变换后得到新的三维数据集YREG，YREG_i，j，k＝V(YCOR，x_i+D_x(ty_i)，y_j+D_y(ty_i)，z_k+D_z(ty_i))；

将两组三维数据集XREG和YREG加权融合，得到V_merged，将V_merged沿z向累加，得到En-face投影图。

2.根据权利要求1所述皮肤光学相干层析图像的配准方法，其特征在于，所述的步骤⑤中预处理包括以下步骤：

①对XRAW和YRAW分别沿z方向做一维中值滤波，得到两组三维数据集XP1和YP1；

②对XP1和YP1分别做对数变换，得到XP2和YP2；

③对XP2和YP2做阈值处理：计算XP2和YP2中所有体素强度的直方图，取直方图中最大值对应的区间的中值作为阈值下限thmin，取所有体素强度的最大值作为阈值上限thmax，对XP2和YP2分别进行阈值处理和归一化得到两组三维数据集，命名为XP3，YP3：

④对XP3和YP3分别进行照明校正：对XP3和YP3分别做沿z方向的累加，得到二维数组Fx和Fy；选择Fx中大于75％体素的强度的强度值作为v_ref，x，选择Fy中大于75％体素的强度的强度值作为v_ref，y，用下式计算照明校正强度二维图Bx和By：

Bx_i，j＝v_ref，x-Fx_i，j

By_i，j＝v_ref，y-Fy_i，j

定义阈值v_thx和v_thy，认为XP3中强度大于v_thx的体素为信号体素，YP3中强度大于v_thy的体素为信号体素；定义三维数据集Mx和My：

对XP3和YP3做下式处理，得到结果XP4和YP4：

⑤对XP4和YP4分别做一次沿z方向系数为2的下采样；将下采样后的数据集化为均值为0，方差为1的数据集，结果分别命名为XFAST和YFAST。

3.根据权利要求1所述皮肤光学相干层析图像的配准方法，其特征在于，所述的纵向配准包括以下步骤：

①将进行纵向配准的两组三维数据集分别命名为X₀和Y₀；定义数据集X₀在y_j＝j处的纵向运动量D_z(tx_j)由平移运动量和旋转运动量组成；平移运动量为从y_j＝j(j＞1)处的理想采样B-scan运动到tx_j(j＞1)时刻的X₀的B-scan在z方向上的平移运动量为D_zd(tx_j)；旋转运动量为从y_j＝j(j＞1)处的理想采样B-scan运动到tx_j(j＞1)时刻的X₀的B-scan在B-scan平面内的旋转运动量D_zt(tx_j)；D_z(tx_j)＝D_zd(tx_j)+D_zt(tx_j)；定义数据集Y₀在x_i＝i处的纵向运动量D_z(ty_i)由平移运动量和旋转运动量组成；平移运动量为从x_i＝i处的理想采样B-scan运动到ty_i时刻的Y₀的B-scan在z方向上平移运动量D_zd(ty_i)；旋转运动量为从x_i＝i处的理想采样B-scan运动到ty_i时刻的Y₀的B-scan在B-scan平面内的旋转运动量D_zt(ty_i)；D_z(ty_i)＝D_zd(ty_i)+D_zt(ty_i)；

②使用插值方法对X₀变换，变换前X_{0 i，j，k}＝V(X₀，x_i，y_j，z_k)，变换后得到新的三维数据集X_g，X_{g i，j，}k＝V(X₀，x_i，y_j，z_k+D_z(tx_j))；使用插值方法对Y₀变换，变换前为Y_{0 i，j，k}＝V(Y₀，x_i，y_j，z_k)，变换后得到新的三维数据集Y_g，Y_{g i，j，k}＝V(Y₀，x_i，y_j，z_k+D_z(ty_i))；

③将变换后的数据集做差，得到残差：

R_i，j，k(X_g，Y_g，D_z)＝X_{g i，j，k}-Y_{g i，j，k}

将一种罚函数L(x)应用于残差，定义相似度项S(D_z)：

定义正则项E(D_z)：

其中tx_l和ty_l分别为采集X₀和Y₀时第1个采样时刻；

将相似度项和正则项组合，得到目标函数：

F(D_z)＝S(D_z)+αE(D_z)，

其中α为正则项相对相似度项重要程度的权值，计算目标函数对于D_z的一阶导数G(D_z)；

④使用迭代优化算法，利用G(D_z)得到纵向运动量D_z，公式如下：

D_z＝argminF(D_z)

其中argmin为求使F(D_z)最小的D_z。

4.根据权利要求1所述皮肤光学相干层析图像的配准方法，其特征在于，所述的整体配准包括以下步骤：

①将进行整体配准的两组三维数据集分别命名为X₀和Y₀；定义数据集X₀在y_j＝j处在x、y、z三个方向的运动量分别为D_x(tx_j)，D_y(tx_j)，D_z(tx_j)，运动量为从y_j＝j(j＞1)处的理想采样B-scan运动到tx_j(j＞1)时刻X₀中B-scan在x、y、z方向上的平移运动量；定义数据集Y₀在x_i＝i处x、y、z三个方向的运动量分别为D_x(ty_i)，D_y(ty_i)，D_z(ty_i)，运动量为从x_i＝i处理想采样B-scan运动到ty_i时刻Y₀中B-scan在x、y、z方向上的平移运动量；

②使用插值方法对X₀变换，变换前x_{0 i，j，k}＝V(X₀，x_i，y_j，z_k)，变换后得到新的三维数据集X_g，为X_{g i，j，k}＝V(X₀，x_i+D_x(tx_j)，y_j+D_y(tx_j)，z_k+D_z(tx_j))；使用插值方法对Y₀变换，变换前为Y_{0 i，j，k}＝V(Y₀，x_i，y_j，z_k)，变换后得到新的三维数据集Y_g，Y_{g i，j，k}＝V(Y₀，x_i+D_x(ty_i)，y_j+D_y(ty_i)，z_k+D_z(ty_i))；

③将变换后的数据集做差，得到残差：

R_i，j，k(X_g，Y_g，D_x，D_y，D_z)＝X_{g i，j，k}-Y_{g i，j，k}

将一种罚函数L(x)应用于残差，定义相似度项S(D_x，D_y，D_z)：

定义正则项E(D_z)：

其中tx_l和ty_l分别为采集X₀和Y₀时第1个采样时刻；

将相似度项和正则项组合，得到目标函数：

F(D_x，D_y，D_z)＝S(D_x，D_y，D_z)+αE(D_x，D_y，D_z)

其中α为正则项相对相似度项重要程度的权值，计算目标函数对于D_x，D_y，D_z的一阶导数G(D_x，D_y，D_z)；

④使用迭代优化算法，利用G(D_x，D_y，D_z)得到整体运动量D_x，D_y，D_z：

D_x，D_y，D_z＝argminF(D_x，D_y，D_z)

其中argmin为求使F(D_x，D_y，D_z)最小的D_x，D_y，D_z。

5.根据权利要求1所述皮肤光学相干层析图像的配准方法，其特征在于，所述的加权融合包括以下步骤：

计算y_j+D_y(tx_j)的概率密度分布和x_i+D_x(ty_i)的概率密度分布，XREG的y_j＝j处B-scan中的A-line的权值为y_j+D_y(tx_j)的概率密度的倒数；YREG的x_i＝i处B-scan中的A-line的权值为对应的x_i+D_x(ty_i)的概率密度的倒数，XREG和YREG的每个对应A-line的权值做归一化，得到融合时的权值。