[发明专利]一种基于迁移学习的超声或CT医学影像三维重建方法

权利要求书

1.一种基于迁移学习的超声或CT医学影像三维重建方法，其特征在于，该方法输入一个超声或者CT影像序列，其影像分辨率为M×N，100≤M≤1500，100≤N≤1500，三维重建的过程具体包括以下步骤：

步骤1：构建数据集

(a)构建自然图像数据集D

选取一个自然图像网站，要求具有图像序列及对应的摄像机内部参数，从所述自然图像网站下载a个图像序列及序列对应的内部参数,1≤a≤20,对于每个图像序列，每相邻3帧图像记为图像b、图像c和图像d，将图像b和图像d按照颜色通道进行拼接，得到图像τ，由图像c与图像τ构成一个数据元素，图像c为自然目标图像，图像c的采样视点作为目标视点，图像b、图像c和图像d的内部参数均为e_t(t＝1，2，3，4)，其中e₁为水平焦距，e₂为垂直焦距，e₃及e₄是主点坐标的两个分量；如果同一图像序列中最后剩余图像不足3帧，则舍弃；利用所有序列构建数据集D，数据集D有f个元素，而且3000≤f≤20000；

(b)构建超声影像数据集E

采样g个超声影像序列，其中1≤g≤20，对于每个序列，每相邻3帧影像记为影像i、影像j和影像k，将影像i和影像k按照颜色通道进行拼接得到影像π，由影像j与影像π构成一个数据元素，影像j为超声目标影像,影像j的采样视点作为目标视点，如果同一影像序列中最后剩余影像不足3帧，则舍弃，利用所有序列构建数据集E，数据集E有F个元素，而且1000≤F≤20000；

(c)构建CT影像数据集G

采样h个CT影像序列，其中1≤h≤20，对于每个序列，每相邻3帧记为影像l、影像m和影像n，将影像l和影像n按照颜色通道进行拼接得到影像σ，由影像m与影像σ构成一个数据元素，影像m为CT目标影像,影像m的采样视点作为目标视点，如果同一影像序列中最后剩余影像不足3帧，则舍弃，利用所有序列构建数据集G，数据集G有ξ个元素，而且1000≤ξ≤20000；步骤2：构建神经网络

神经网络处理的图像或影像的分辨率均为p×o，p为宽度，o为高度，以像素为单位，100≤o≤2000，100≤p≤2000；

(1)网络A的结构

张量H作为输入，尺度为α×o×p×3，张量I作为输出,尺度为α×o×p×1，α为批次数量；

网络A由编码器和解码器组成，对于张量H，依次经过编码和解码处理后，获得输出张量I；

编码器由5个残差单元组成，第1至5个单元分别包括2，3，4，6，3个残差模块，每个残差模块进行3次卷积，卷积核的形状均为3×3，卷积核的个数分别是64，64，128，256，512，其中，第一个残差单元之后，包含一个最大池化层；

解码器由6个解码单元组成，每个解码单元均包括反卷积和卷积两步处理，反卷积和卷积处理的卷积核形状、个数都相同，第1至6个解码单元卷积核的形状均为3×3，卷积核的个数分别是512，256，128，64，32，16，编码器与解码器的网络层之间进行跨层连接，跨层连接的对应关系为：1与4、2与3、3与2、4与1；

(2)网络B的结构

张量J和张量K作为输入，尺度分别为α×o×p×3和α×o×p×6，张量L和张量O作为输出，尺度分别为α×2×6和α×4×1，α为批次数量；

网络B由模块P和模块Q构成，共有11层卷积单元,首先，将张量J和张量K按照最后一个通道进行拼接，获得尺度为α×o×p×9的张量，对于所述张量经过模块P和模块Q处理后，分别获得输出张量L和张量O；

模块Q与模块P共享前4层卷积单元，前4层卷积单元结构为：前两层单元中卷积核尺度分别为7×7、5×5，第3层到第4层的卷积核尺度均为3×3，1至4层的卷积核个数依次为16、32、64、128；

对于模块P，除了共享4层外，占据网络B的第5层到第7层卷积单元，卷积核尺度均为3×3,卷积核个数均为256，第7层的处理结果利用12个3×3的卷积核进行卷积处理后，将12个结果顺次排成2行，得到张量L的结果；

对于模块Q，除了共享网络B的1至4层外，还占据网络B的第8至11层卷积单元，网络B的第2层输出作为网络B的第8层输入，第8至11层卷积单元中卷积核的形状均为3×3，卷积核个数均为256，第11层的结果利用4个3×3的卷积核进行卷积处理后，从4个通道得到张量O的结果；

(3)网络C的结构

张量R和张量S作为网络输入，尺度均为α×o×p×3，张量T作为网络输出，尺度为α×o×p×2，α为批次数量；

网络C设计为编码和解码结构，首先，将张量R和张量S按照最后一个通道进行拼接，获得尺度为α×o×p×6的张量，对于所述张量经过编码和解码处理后，获得输出张量T；

对于编码结构，由6层编码单元组成，每层编码单元包含1次卷积处理，1次批归一化处理和1次激活处理，其中第1层编码单元采用7×7卷积核，其它层编码单元均采用3×3卷积核，第1和3层编码单元的卷积步长为1，其它层卷积步长均为2，对于每层编码单元，均采用Relu函数激活，1-6层编码单元的卷积核个数分别为16、32、64、128、256、512；

对于解码结构，由6层解码单元组成，每层解码单元由反卷积单元、连接处理单元和卷积单元构成，其中反卷积单元包括反卷积处理与Relu激活处理，1-6层反卷积核的大小均为3x3，对于第1-2层解码单元，反卷积步长为1，3-6层解码单元的反卷积步长为2，1-6层反卷积核的个数依次为512、256、128、64、32、16，连接处理单元将编码单元和对应解码单元的反卷积结果连接后，输入到卷积单元，1-5层卷积单元的卷积核大小为3x3，第6层卷积单元的卷积核大小为7x7，1-6层卷积单元的卷积步长均为2，将第6层的卷积结果经过2个3x3的卷积处理后，得到结果T；

(4)网络μ的结构

张量Ω作为网络输入，尺度为α×o×p×3，张量作为网络输出，尺度为α×o×p×1，α为批次数量；

网络μ由编码器和解码器组成，对于张量Ω，依次经过编码和解码处理后，获得输出张量

编码器由14层编码单元组成，每个编码单元包含1次卷积处理，1次批归一化处理和1次激活处理，第1和2层的编码单元均采用7×7大小卷积核结构，第3和4层的编码单元采用5×5大小卷积核结构，其余编码单元设为3×3的卷积核，每个编码单元的卷积步长设计为：第1、3、5、7、9、11和13层的步长设为2，其它层的步长设为1，对于每个编码单元都采用Relu函数激活处理，编码结构中，第1层至第8层的卷积核个数分别设为32、32、64、64、128、128、256、256，其余层的卷积核个数均为512；

解码器由7层解码单元组成，每层解码单元由反卷积单元、连接处理单元和卷积单元构成，在任意一层的解码单元中，反卷积单元包括反卷积处理与Relu激活处理，各层反卷积核的大小均为3x3，步长均为2，1至7层反卷积核的个数分别为512、512、256、128、64、32、16，连接处理单元将编码单元和对应层的反卷积特征连接后，输入到下一卷积单元进行卷积和Relu激活处理，在卷积单元中，1至7层卷积核的个数分别为512、512、256、128、64、32、16，在各卷积单元中，卷积核大小均为3x3，步长均为1，4至7层解码单元的输出分别乘以权重后得到输出结果张量

步骤3：神经网络的训练

分别将数据集D、数据集E和数据集G中样本按照9:1划分为训练集和测试集，训练集中数据用于训练，测试集数据用于测试，在下列各步骤训练时，分别从对应的数据集中获取训练数据，统一缩放到分辨率p×o，输入到对应网络中，迭代优化，通过不断修改网络模型参数，使得每批次的损失达到最小；

在训练过程中，各损失的计算方法：

内部参数监督合成损失：在自然图像的网络模型参数训练中，将网络A的输出张量I作为深度，将网络B的输出结果L与训练数据的内部参数标签e_t(t＝1，2，3，4)分别作为位姿参数和摄像机内部参数，根据计算机视觉算法，利用图像b和图像d分别合成图像c视点处的两个图像，利用图像c分别与所述的两个图像，按照逐像素、逐颜色通道强度差之和计算得到；

无监督合成损失：在超声或者CT影像的网络模型参数训练中，将网络μ的输出张量作为深度，将网络B的输出张量L和张量O分别作为位姿参数和摄像机内部参数，根据计算机视觉算法，利用目标影像的两相邻影像分别合成目标影像视点处的影像，利用目标影像分别与所述目标影像视点处的影像，按照逐像素、逐颜色通道强度差之和计算得到；

内部参数误差损失：利用网络B的输出结果O与训练数据的内部参数标签e_t(t＝1，2，3，4)按照各分量差的绝对值之和计算得到；

空间结构误差损失：在超声或者CT影像的网络模型参数训练中，将网络μ的输出张量作为深度，将网络B的输出张量L和张量O分别作为位姿参数和摄像机内部参数，根据计算机视觉算法，以目标影像的视点为摄像机坐标系的原点，对目标影像进行重建，采用RANSAC算法对重建点进行空间结构的拟合，利用目标影像的每个重建点到空间几何结构的欧式距离计算得到；

变换合成损失：在超声或者CT影像的网络模型参数训练中，将网络μ的输出张量作为深度，将网络B的输出张量L和张量O分别作为位姿参数和摄像机内部参数，根据计算机视觉算法，利用目标影像的两个相邻影像合成目标影像视点处的两个影像的过程中，对于所述合成的两个影像中的每个影像，在得到每个像素位置后，将每个像素的坐标加上网络C输出的每个像素位移结果，得到每个像素新的位置，构成合成结果影像，利用该合成结果影像与影像j之间的逐像素、逐颜色通道强度差之和计算得到；

(1)在数据集D上，分别对网络A及网络B的模块P训练80000次

每次从数据集D中取出训练数据，统一缩放到分辨率p×o，将图像c输入网络A，将图像c及图像τ输入网络B，对网络B的模块P进行训练，每批次的训练损失由内部参数监督合成损失计算得到；

(2)在数据集D上，对网络B的模块Q训练80000次

每次从数据集D中取出训练数据，统一缩放到分辨率p×o，将图像c输入网络A，将图像c及图像τ输入网络B，对网络B的模块Q进行训练，每批次的训练损失由内部参数监督合成损失和内部参数误差损失之和计算得到；

(3)在数据集E上，对网络μ及网络B的模块Q训练80000次，进行特征迁移

每次从数据集E中取出超声训练数据，统一缩放到分辨率p×o，将影像j输入网络μ，将影像j及影像π输入网络B，对网络B的模块Q进行训练，每批次的训练损失计算如下：

z＝v+W+χ (1)

其中v为无监督合成损失，W为空间结构误差损失，常深度损失χ利用网络μ的输出结果的均方差计算得到；

(4)在数据集E上，按照如下步骤对网络B两个模块训练80000次

每次从数据集E中取出超声训练数据，统一缩放到分辨率p×o，将影像j输入网络μ，将影像j及影像π输入网络B，训练过程中调整网络B的两个模块参数，迭代优化，使得每批次的每幅影像的损失达到最小，每批次的训练损失由无监督合成损失、空间结构误差损失和常深度损失之和组成，常深度损失利用网络μ的输出结果的均方差计算得到；

(5)在数据集E上，对网络C和网络B的两个模块训练80000次

每次从数据集E中取出超声影像训练数据，统一缩放到分辨率p×o，将影像j输入网络μ，将影像j及影像π输入网络B，将网络μ的输出作为深度，将网络B的输出作为位姿参数和摄像机内部参数，分别根据影像i和影像k合成影像j视点处的两张影像，将所述的两张影像输入网络C，通过不断修改网络C和网络B的参数，迭代优化，使得每批次的每幅影像的损失达到最小，每批次的损失计算为变换合成损失、空间结构误差损失和常深度损失之和，其中常深度损失利用网络μ的输出结果的均方差计算得到；

(6)在数据集E上，对网络C和网络B的两个模块训练50000次，得到模型ρ

每次从数据集E中取出超声影像训练数据，统一缩放到分辨率p×o，将影像j输入网络μ，将影像j及影像π输入网络B，将网络μ的输出作为深度，将网络B的输出作为位姿参数和摄像机内部参数，分别根据影像i和影像k合成影像j视点处的两张影像，将所述的两张影像输入网络C，通过不断修改网络C和网络B的参数，迭代优化，使得每批次的每幅影像的损失达到最小，迭代后得到最优的网络模型参数ρ，每批次的损失计算为变换合成损失、空间结构误差损失之和；

(7)在数据集G上，对网络C和网络B训练80000次

每次从数据集G中取出CT影像训练数据，统一缩放到分辨率p×o，将影像m输入网络μ，将影像m及影像σ输入网络B，将网络μ的输出作为深度，将网络B的输出作为位姿参数和摄像机内部参数，分别根据影像l和影像n合成影像m视点处的两张影像，将所述的两张影像输入网络C，通过不断修改网络C和网络B的参数，迭代优化，使得每批次的每幅影像的损失达到最小，每批次的损失计算为变换合成损失、空间结构误差损失、常深度损失和摄象机平移运动损失Y之和，其中，常深度损失利用网络μ的输出结果的均方差计算得到，Y由网络B的输出位姿参数，根据摄像机平移运动的约束计算得到；

(8)在数据集G上，对网络C和网络B训练50000次，得到模型ρ′

每次从数据集G中取出CT影像训练数据，统一缩放到分辨率p×o，将影像m输入网络μ，将影像m及影像σ输入网络B，将网络μ的输出作为深度，将网络B的输出作为位姿参数和摄像机内部参数，分别根据影像l和影像n合成影像m视点处的两张影像，将所述的两张影像输入网络C，通过不断修改网络C和网络B的参数，迭代优化，使得每批次的每幅影像的损失达到最小，迭代后得到最优的网络模型参数ρ′，每批次的损失计算为变换合成损失、空间结构误差损失和摄象机平移运动损失Y之和，Y由网络B的输出位姿参数，根据摄像机平移运动的约束计算得到；

步骤4：超声或者CT影像三维重建

利用自采样的一个超声或者CT序列影像，将每一帧影像统一缩放到分辨率p×o，使用模型参数ρ或者模型参数ρ′进行预测，对于超声序列影像，将影像j输入网络μ，将影像j及影像π输入网络B，对于CT序列影像，将影像m输入网络μ，将影像m及影像σ输入网络B，将网络μ的输出作为深度，将网络B的输出作为位姿参数和摄像机内部参数，按照下列步骤选取关键帧，序列中第一帧作为当前关键帧，依次将序列影像中的每一帧作为目标帧，根据当前关键帧，利用摄像机位姿参数和内部参数，合成目标帧视点处的影像，利用所述合成影像与目标帧之间逐像素逐颜色通道强度差之和的大小计算误差λ，再根据目标帧的相邻帧，利用摄像机位姿参数和内部参数，合成目标帧视点处的影像，利用所述合成影像与目标帧之间逐像素逐颜色通道强度差之和的大小计算误差γ，进一步利用公式(2)计算合成误差比Z，当Z大于阈值η时，1η2,将当前关键帧更新为此时的目标帧；

对任一目标帧，将其分辨率缩放到M×N，根据摄像机内部参数，依据计算机视觉的重建算法，计算每帧影像每个像素的摄象机坐标系中的三维坐标，进一步，将第一帧的视点作为世界坐标系的原点，再结合所有关键帧的位姿参数，利用三维空间几何变换，计算得到该序列每一帧影像每个像素的世界坐标系中的三维坐标。