[發(fā)明專利]一種基于知識蒸餾的超聲或CT醫(yī)學影像的三維重建方法

有效

• 申請?zhí)枺?/td>: 202011621388.9
• 申請日: 2020-12-30
• 公開（公告）號: CN112734906B
• 公開（公告）日: 2022-08-19
• 發(fā)明（設計）人: 全紅艷;錢笑笑
• 申請（專利權）人: 華東師范大學
• 主分類號: G06T17/00
• 分類號: G06T17/00;G06N3/04;G06N3/08
• 代理公司: 上海藍迪專利商標事務所(普通合伙) 31215
• 代理人: 徐筱梅;張翔

1.一種基于知識蒸餾的超聲或CT醫(yī)學影像的三維重建方法，其特征在于，該方法輸入一個超聲或者CT影像序列，其影像分辨率為M×N，100≤M≤1500，100≤N≤1500，三維重建的過程具體包括以下步驟：

步驟1：構建數(shù)據(jù)集

(a)構建自然圖像數(shù)據(jù)集D

選取一個自然圖像網(wǎng)站，要求具有圖像序列及對應的攝像機內部參數(shù)，從所述自然圖像網(wǎng)站下載a個圖像序列及序列對應的內部參數(shù),1≤a≤20,對于每個圖像序列，每相鄰3幀圖像記為圖像b、圖像c和圖像d，將圖像b和圖像d按照顏色通道進行拼接，得到圖像τ，由圖像c與圖像τ構成一個數(shù)據(jù)元素，圖像c為自然目標圖像，圖像c的采樣視點作為目標視點，圖像b、圖像c和圖像d的內部參數(shù)均為e_t(t＝1，2，3，4)，其中e₁為水平焦距，e₂為垂直焦距，e₃及e₄是主點坐標的兩個分量；如果同一圖像序列中最后剩余圖像不足3幀，則舍棄；利用所有序列構建數(shù)據(jù)集D，數(shù)據(jù)集D有f個元素，而且3000≤f≤20000；

(b)構建超聲影像數(shù)據(jù)集E

采樣g個超聲影像序列，其中1≤g≤20，對于每個序列，每相鄰3幀影像記為影像i、影像j和影像k，將影像i和影像k按照顏色通道進行拼接得到影像π，由影像j與影像π構成一個數(shù)據(jù)元素，影像j為超聲目標影像,影像j的采樣視點作為目標視點，如果同一影像序列中最后剩余影像不足3幀，則舍棄，利用所有序列構建數(shù)據(jù)集E，數(shù)據(jù)集E有F個元素，而且1000≤F≤20000；

(c)構建CT影像數(shù)據(jù)集G

采樣h個CT影像序列，其中1≤h≤20，對于每個序列，每相鄰3幀記為影像l、影像m和影像n，將影像l和影像n按照顏色通道進行拼接得到影像σ，由影像m與影像σ構成一個數(shù)據(jù)元素，影像m為CT目標影像,影像m的采樣視點作為目標視點，如果同一影像序列中最后剩余影像不足3幀，則舍棄，利用所有序列構建數(shù)據(jù)集G，數(shù)據(jù)集G有ξ個元素，而且1000≤ξ≤20000；步驟2：構建神經(jīng)網(wǎng)絡

神經(jīng)網(wǎng)絡處理的圖像或影像的分辨率均為p×o，p為寬度，o為高度，以像素為單位，100≤o≤2000，100≤p≤2000；

(1)網(wǎng)絡A的結構

張量H作為輸入，尺度為α×o×p×3，張量I作為輸出,尺度為α×o×p×1，α為批次數(shù)量；

網(wǎng)絡A由編碼器和解碼器組成，對于張量H，依次經(jīng)過編碼和解碼處理后，獲得輸出張量I；

編碼器由5個殘差單元組成，第1至5個單元分別包括2，3，4，6，3個殘差模塊，每個殘差模塊進行3次卷積，卷積核的形狀均為3×3，卷積核的個數(shù)分別是64，64，128，256，512，其中，第一個殘差單元之后，包含一個最大池化層；

解碼器由6個解碼單元組成，每個解碼單元均包括反卷積和卷積兩步處理，反卷積和卷積處理的卷積核形狀、個數(shù)都相同，第1至6個解碼單元卷積核的形狀均為3×3，卷積核的個數(shù)分別是512，256，128，64，32，16，編碼器與解碼器的網(wǎng)絡層之間進行跨層連接，跨層連接的對應關系為：1與4、2與3、3與2、4與1；

(2)網(wǎng)絡B的結構

張量J和張量K作為輸入，尺度分別為α×o×p×3和α×o×p×6，張量L、張量O和張量作為輸出，尺度分別為α×2×6，α×4×1和α×1×1，α為批次數(shù)量；

網(wǎng)絡B由模塊P、模塊Q和模塊μ構成，共有16層卷積單元，首先，將張量J和張量K按照最后一個通道進行拼接，獲得尺度為α×o×p×9的張量，對于所述張量經(jīng)過模塊P、模塊Q和模塊μ處理后，分別獲得輸出張量L、張量O和張量

模塊Q、模塊P和模塊μ共享前4層卷積單元，這4層結構為：前兩層單元中卷積核尺度分別為7×7、5×5，第3層到第4層的卷積核尺度均為3×3，這4層結構的卷積核個數(shù)依次為16、32、64、128；

對于模塊P，除了共享4層外，占據(jù)網(wǎng)絡B的第5層到第7層卷積單元，卷積核尺度均為3×3，卷積核個數(shù)均為256，第7層的處理結果利用12個3×3的卷積核進行卷積處理后，從12個通道得到張量L的結果；

對于模塊Q，除了共享網(wǎng)絡B的1至4層外，還占據(jù)網(wǎng)絡B的第8至11層卷積單元，網(wǎng)絡B的第2層輸出作為網(wǎng)絡B的第8層輸入，第8至11層卷積單元中卷積核的形狀均為3×3，卷積核個數(shù)均為256，第11層的結果利用4個3×3的卷積核進行卷積處理后，從4個通道得到張量O的結果；

對于模塊μ，除了共享網(wǎng)絡B的1至4層外，還占據(jù)網(wǎng)絡B的第12層至第16層卷積單元，網(wǎng)絡B的第2層輸出作為第12層輸入，第12層至第16層中卷積核的形狀均為3×3，卷積核個數(shù)均為256，第16層的結果利用1個3×3的卷積核進行卷積處理后，得到張量的結果；

(3)網(wǎng)絡C的結構

張量R和張量S作為網(wǎng)絡輸入，尺度均為α×o×p×3，張量T作為網(wǎng)絡輸出，尺度為α×o×p×2，α為批次數(shù)量；

網(wǎng)絡C設計為編碼和解碼結構，首先，將張量R和張量S按照最后一個通道進行拼接，獲得尺度為α×o×p×6的張量，對于所述張量經(jīng)過編碼和解碼處理后，獲得輸出張量T；

對于編碼結構，由6層編碼單元組成，每層編碼單元包含1次卷積處理，1次批歸一化處理和1次激活處理，其中第1層編碼單元采用7×7卷積核，其它層編碼單元均采用3×3卷積核，第1和3層編碼單元的卷積步長為1，其它層卷積步長均為2，對于每層編碼單元，均采用Relu函數(shù)激活，1-6層編碼單元的卷積核個數(shù)分別為16、32、64、128、256、512；

對于解碼結構，由6層解碼單元組成，每層解碼單元由反卷積單元、連接處理單元和卷積單元構成，其中反卷積單元包括反卷積處理與Relu激活處理，1-6層反卷積核的大小均為3x3，對于第1-2層解碼單元，反卷積步長為1，3-6層解碼單元的反卷積步長為2，1-6層反卷積核的個數(shù)依次為512、256、128、64、32、16，連接處理單元將編碼單元和對應解碼單元的反卷積結果連接后，輸入到卷積單元，1-5層卷積單元的卷積核大小為3x3，第6層卷積單元的卷積核大小為7x7，1-6層卷積單元的卷積步長均為2，將第6層的卷積結果經(jīng)過2個3x3的卷積處理后，得到結果T；

步驟3：神經(jīng)網(wǎng)絡的訓練

分別將數(shù)據(jù)集D、數(shù)據(jù)集E和數(shù)據(jù)集G中樣本按照9:1劃分為訓練集和測試集，訓練集中數(shù)據(jù)用于訓練，測試集數(shù)據(jù)用于測試，在下列各步驟訓練時，分別從對應的數(shù)據(jù)集中獲取訓練數(shù)據(jù)，統(tǒng)一縮放到分辨率p×o，輸入到對應網(wǎng)絡中，迭代優(yōu)化，通過不斷修改網(wǎng)絡模型參數(shù)，使得每批次的損失達到最小；

在訓練過程中，各損失的計算方法：

內部參數(shù)監(jiān)督合成損失：在自然圖像的網(wǎng)絡模型參數(shù)訓練中，將網(wǎng)絡A的輸出張量I作為深度，將網(wǎng)絡B的輸出結果L與訓練數(shù)據(jù)的內部參數(shù)標簽e_t(t＝1，2，3，4)分別作為位姿參數(shù)和攝像機內部參數(shù)，根據(jù)計算機視覺算法，利用圖像b和圖像d分別合成圖像c視點處的兩個圖像，利用圖像c分別與所述的兩個圖像，按照逐像素、逐顏色通道強度差之和計算得到；

無監(jiān)督合成損失：在超聲或者CT影像的網(wǎng)絡模型參數(shù)訓練中，將網(wǎng)絡B的模塊μ的輸出張量作為深度，將網(wǎng)絡B的輸出張量L和張量O分別作為位姿參數(shù)和攝像機內部參數(shù)，根據(jù)計算機視覺算法，利用目標影像的兩相鄰影像分別合成目標影像視點處的影像，利用目標影像分別與所述目標影像視點處的影像，按照逐像素、逐顏色通道強度差之和計算得到；

內部參數(shù)誤差損失：利用網(wǎng)絡B的輸出結果O與訓練數(shù)據(jù)的內部參數(shù)標簽e_t(t＝1，2，3，4)按照各分量差的絕對值之和計算得到；

空間結構誤差損失：在超聲或者CT影像的網(wǎng)絡模型參數(shù)訓練中，將網(wǎng)絡B的模塊μ的輸出張量作為深度，將網(wǎng)絡B的輸出張量L和張量O分別作為位姿參數(shù)和攝像機內部參數(shù)，根據(jù)計算機視覺算法，以目標影像的視點為攝像機坐標系的原點，對目標影像進行重建，采用RANSAC算法對重建點進行空間結構的擬合，利用目標影像的每個重建點到空間幾何結構的歐式距離計算得到；

變換合成損失：在超聲或者CT影像的網(wǎng)絡模型參數(shù)訓練中，將網(wǎng)絡B的模塊μ的輸出張量作為深度，將網(wǎng)絡B的輸出張量L和張量O分別作為位姿參數(shù)和攝像機內部參數(shù)，根據(jù)計算機視覺算法，利用目標影像的兩個相鄰影像合成目標影像視點處的兩個影像的過程中，對于所述合成的兩個影像中的每個影像，在得到每個像素位置后，將每個像素的坐標加上網(wǎng)絡C輸出的每個像素位移結果，得到每個像素新的位置，構成合成結果影像，利用該合成結果影像與影像j之間的逐像素、逐顏色通道強度差之和計算；

(1)在數(shù)據(jù)集D上，分別對網(wǎng)絡A及網(wǎng)絡B的模塊P訓練80000次

每次從數(shù)據(jù)集D中取出訓練數(shù)據(jù)，統(tǒng)一縮放到分辨率p×o，將圖像c輸入網(wǎng)絡A，將圖像c及圖像τ輸入網(wǎng)絡B，對網(wǎng)絡B的模塊P進行訓練，每批次的訓練損失由內部參數(shù)監(jiān)督合成損失計算得到；

(2)在數(shù)據(jù)集D上，對網(wǎng)絡B的模塊Q訓練80000次

每次從數(shù)據(jù)集D中取出訓練數(shù)據(jù)，統(tǒng)一縮放到分辨率p×o，將圖像c輸入網(wǎng)絡A，將圖像c及圖像τ輸入網(wǎng)絡B，對網(wǎng)絡B的模塊Q進行訓練，每批次的訓練損失由內部參數(shù)監(jiān)督合成損失和內部參數(shù)誤差損失之和計算得到；

(3)在數(shù)據(jù)集E上，對網(wǎng)絡B的模塊Q和模塊μ進行訓練80000次，進行知識蒸餾

每次從數(shù)據(jù)集E中取出超聲訓練數(shù)據(jù)，統(tǒng)一縮放到分辨率p×o，將影像j及影像π輸入網(wǎng)絡B，對網(wǎng)絡B的模塊Q和模塊μ進行訓練，每批次的訓練損失計算如下：

z＝v+W+χ (1)

其中，v為無監(jiān)督合成損失，W為空間結構誤差損失，常深度損失χ利用模塊μ的輸出結果的均方差計算得到；

(4)在數(shù)據(jù)集E上，按照如下步驟對網(wǎng)絡B的三個模塊訓練80000次

每次從數(shù)據(jù)集E中取出超聲訓練數(shù)據(jù)，統(tǒng)一縮放到分辨率p×o，將影像j及影像π輸入網(wǎng)絡B，訓練時，通過不斷修改網(wǎng)絡B的三個模塊的參數(shù)，迭代優(yōu)化，使得每批次的每幅影像的損失達到最小，每批次的訓練損失由無監(jiān)督合成損失、空間結構誤差損失和常深度損失之和組成，常深度損失利用網(wǎng)絡B的模塊μ的輸出結果的均方差計算得到；

(5)在數(shù)據(jù)集E上，對網(wǎng)絡C和網(wǎng)絡B的三個模塊訓練80000次

每次從數(shù)據(jù)集E中取出超聲影像訓練數(shù)據(jù)，統(tǒng)一縮放到分辨率p×o，將影像j及影像π輸入網(wǎng)絡B，將網(wǎng)絡B的模塊μ的輸出張量作為深度，將網(wǎng)絡B輸出張量L和張量O分別作為位姿參數(shù)和攝像機內部參數(shù)，分別根據(jù)影像i和影像k合成影像j視點處的兩張影像，將所述的兩張影像輸入網(wǎng)絡C，通過不斷修改網(wǎng)絡C和網(wǎng)絡B的參數(shù)，迭代優(yōu)化，使得每批次的每幅影像的損失達到最小，每批次的損失計算為變換合成損失、空間結構誤差損失和常深度損失之和，其中常深度損失利用網(wǎng)絡B的模塊μ的輸出結果的均方差計算得到；

(6)在數(shù)據(jù)集E上，對網(wǎng)絡C和網(wǎng)絡B的三個模塊訓練50000次，得到模型ρ

訓練時，每次從數(shù)據(jù)集E中取出超聲影像訓練數(shù)據(jù)，統(tǒng)一縮放到分辨率p×o，將影像j及影像π輸入網(wǎng)絡B，將網(wǎng)絡B的模塊μ的輸出張量作為深度，將網(wǎng)絡B輸出張量L和張量O分別作為位姿參數(shù)和攝像機內部參數(shù)，分別根據(jù)影像i和影像k合成影像j視點處的兩張影像，將所述的兩張影像輸入網(wǎng)絡C，通過不斷修改網(wǎng)絡C和網(wǎng)絡B的參數(shù)，迭代優(yōu)化，使得每批次的每幅影像的損失達到最小，迭代后得到最優(yōu)的網(wǎng)絡模型參數(shù)ρ，每批次的損失計算為變換合成損失和空間結構誤差損失之和；

(7)在數(shù)據(jù)集G上，對網(wǎng)絡C和網(wǎng)絡B的三個模塊訓練80000次

訓練時，每次從數(shù)據(jù)集G中取出CT影像訓練數(shù)據(jù)，統(tǒng)一縮放到分辨率p×o，將影像m及影像σ輸入網(wǎng)絡B，將網(wǎng)絡B的模塊μ的輸出張量作為深度，將網(wǎng)絡B輸出張量L和張量O分別作為位姿參數(shù)和攝像機內部參數(shù)，分別根據(jù)影像l和影像n合成影像m視點處的兩張影像，將所述的兩張影像輸入網(wǎng)絡C，通過不斷修改網(wǎng)絡C和網(wǎng)絡B的參數(shù)，迭代優(yōu)化，使得每批次的每幅影像的損失達到最小，每批次的損失計算為變換合成損失、空間結構誤差損失、常深度損失和攝象機平移運動損失Y之和，其中，常深度損失利用網(wǎng)絡B的模塊μ的輸出結果的均方差計算得到，Y由網(wǎng)絡B的輸出位姿參數(shù)，根據(jù)攝像機平移運動的約束計算得到；

(8)在數(shù)據(jù)集G上，對網(wǎng)絡C和網(wǎng)絡B的三個模塊訓練50000次，得到模型ρ′

每次從數(shù)據(jù)集G中取出CT影像訓練數(shù)據(jù)，統(tǒng)一縮放到分辨率p×o，將影像m及影像σ輸入網(wǎng)絡B，將網(wǎng)絡B的模塊μ的輸出張量作為深度，將網(wǎng)絡B輸出張量L和張量O分別作為位姿參數(shù)和攝像機內部參數(shù)，分別根據(jù)影像l和影像n合成影像m視點處的兩張影像，將所述的兩張影像輸入網(wǎng)絡C，通過不斷修改網(wǎng)絡C和網(wǎng)絡B的參數(shù)，迭代優(yōu)化，使得每批次的每幅影像的損失達到最小，迭代后得到最優(yōu)的網(wǎng)絡模型參數(shù)ρ′，每批次的損失計算為變換合成損失、空間結構誤差損失和攝象機平移運動損失Y之和，Y由網(wǎng)絡B的輸出位姿參數(shù)，根據(jù)攝像機平移運動的約束計算得到；

步驟4：超聲或者CT影像三維重建

利用自采樣的一個超聲或者CT序列影像，將每一幀影像統(tǒng)一縮放到分辨率p×o，使用模型參數(shù)ρ或者模型參數(shù)ρ′進行預測，對于超聲序列影像，將影像j及影像π輸入網(wǎng)絡B，對于CT序列影像，將影像m及影像σ輸入網(wǎng)絡B，將網(wǎng)絡B的模塊μ的輸出張量作為深度，將網(wǎng)絡B的輸出張量L和張量O分別作為位姿參數(shù)和攝像機內部參數(shù)，按照下列步驟選取關鍵幀，序列中第一幀作為當前關鍵幀，依次將序列影像中的每一幀作為目標幀，根據(jù)當前關鍵幀，利用攝像機位姿參數(shù)和內部參數(shù)，合成目標幀視點處的影像，利用所述合成影像與目標幀之間逐像素逐顏色通道強度差之和的大小計算誤差λ，再根據(jù)目標幀的相鄰幀，利用攝像機位姿參數(shù)和內部參數(shù)，合成目標幀視點處的影像，利用所述合成影像與目標幀之間逐像素逐顏色通道強度差之和的大小計算誤差γ，進一步利用公式(2)計算合成誤差比Z，當Z大于閾值η時，1η2,將當前關鍵幀更新為此時的目標幀；

對任一目標幀，將其分辨率縮放到M×N，根據(jù)攝像機內部參數(shù)，依據(jù)計算機視覺的重建算法，計算每幀影像每個像素的攝象機坐標系中的三維坐標，進一步，將第一幀的視點作為世界坐標系的原點，再結合所有關鍵幀的位姿參數(shù)，利用三維空間幾何變換，計算得到該序列每一幀影像每個像素的世界坐標系中的三維坐標。