陳繼剛，王曉康，康永興，關(guān)亞彬，董學(xué)剛，張子路

(1. 燕山大學(xué)機(jī)械工程學(xué)院 河北 秦皇島 066004；2. 東北大學(xué)機(jī)械工程與自動(dòng)化學(xué)院 沈陽 110819)

骨組織是人體重要的組織之一，屬于獨(dú)特的多孔結(jié)構(gòu)[1]。在骨組織工程領(lǐng)域，三維骨多孔結(jié)構(gòu)的數(shù)字建模是諸多研究的基礎(chǔ)技術(shù)[2]。當(dāng)前，對(duì)于具有仿生特征的三維骨多孔結(jié)構(gòu)數(shù)字模型，大多還是利用骨切片圖像(以下稱：切片圖像)進(jìn)行三維重構(gòu)得到的。

文獻(xiàn)[3-8]提出兩點(diǎn)相關(guān)函數(shù)法[3]、分形蒙特卡洛法[4]、多點(diǎn)統(tǒng)計(jì)法[5]、基于深度學(xué)習(xí)的方法[6-8]等切片圖像訓(xùn)練方法，同時(shí)基于切片圖像(包括重構(gòu)生成的圖像)建模得到滿意的三維骨多孔結(jié)構(gòu)數(shù)字模型。由于成像設(shè)備精度、存儲(chǔ)設(shè)備容量等原因，切片圖像在獲取過程中會(huì)出現(xiàn)數(shù)據(jù)丟失、圖像受損或圖像尺寸異常等問題，導(dǎo)致只能得到具有局部特征信息的圖像[9]。另外，雖然計(jì)算機(jī)斷層掃描技術(shù)[10]和掃描電子顯微鏡[11]等成像技術(shù)取得了進(jìn)步，但在實(shí)際情況下，獲取大區(qū)域、完整的切片圖像還是比較困難的[12]。本文把受損的切片圖像以及小尺寸切片圖像稱為局部骨切片圖像(以下稱：局部切片圖像)，即只包含局部骨多孔特征的切片圖像。由于局部切片圖像中的數(shù)據(jù)是不完整的，所以其無法在統(tǒng)計(jì)意義和實(shí)際意義上表達(dá)三維結(jié)構(gòu)的主要特征。因此，為了得到準(zhǔn)確的三維骨多孔結(jié)構(gòu)數(shù)字模型，對(duì)局部切片圖像進(jìn)行修復(fù)并重構(gòu)，得到完整且相似的圖像是很關(guān)鍵的。

近年來，研究人員對(duì)局部切片圖像的修復(fù)重構(gòu)進(jìn)行了深入研究，其中主要有以下3 種方法：逐像素法[13]、逐區(qū)塊法[14]和基于深度學(xué)習(xí)的方法。逐像素法是基于fast marching method (FMM)的修復(fù)算法，逐區(qū)塊法則是基于exemplar-based inpainting的算法，這兩種方法都是通過統(tǒng)計(jì)單張局部圖像的像素特征來推測(cè)全局信息，其重構(gòu)修復(fù)得到的完整圖像與真實(shí)切片圖像不具有很好的相似性。而基于深度學(xué)習(xí)的方法可以利用局部圖像的局部特征，生成與真實(shí)切片圖像特征相似的完整圖像，有著更好的重構(gòu)效果，如文獻(xiàn)[7]應(yīng)用了條件生成對(duì)抗網(wǎng)絡(luò)實(shí)現(xiàn)了對(duì)多種局部多孔圖像的重構(gòu)，并運(yùn)用統(tǒng)計(jì)相關(guān)函數(shù)證明了生成圖像的合理性。

本文在條件生成對(duì)抗網(wǎng)絡(luò)(pix2pix)[15]的基礎(chǔ)上，對(duì)生成器進(jìn)行改進(jìn)，加入嵌套殘差密集塊(residual-in-residual dense block, RRDB)[16]，同 時(shí)，在判別器中加入極化自注意力模塊(polarized selfattention, PSA)[17]，提出改進(jìn)條件生成對(duì)抗網(wǎng)絡(luò)(RRDB-PSA-pix2pix)對(duì)局部切片圖像進(jìn)行完整重構(gòu)，并通過性能相似性與優(yōu)異性分析，驗(yàn)證了重構(gòu)生成的圖像與真實(shí)骨多孔圖像的相似程度。

1 網(wǎng)絡(luò)原理與改進(jìn)方法

1.1 條件生成對(duì)抗網(wǎng)絡(luò)概述

條件生成對(duì)抗網(wǎng)絡(luò)(conditional generative adversarial networks, cGAN)是生成對(duì)抗網(wǎng)絡(luò)[18](generative adversarial network, GAN)的擴(kuò)展。

cGAN 主要由兩部分構(gòu)成：一個(gè)是生成器(generator network)，用來生成數(shù)據(jù)；另一個(gè)是判別器(discriminator network)，用來分辨原始數(shù)據(jù)和生成數(shù)據(jù)的真?zhèn)巍Ｉ善骱团袆e器相互訓(xùn)練對(duì)抗，從而使生成的數(shù)據(jù)越來越接近真實(shí)數(shù)據(jù)。另外，由于cGAN 的生成器輸入的是隨機(jī)噪聲和約束條件，因此，判別器需要加上約束條件一起辨別。

近年來，適用于cGAN 的約束條件越來越多，因此誕生出許多cGAN 的變體[19-21]，這使得cGAN的應(yīng)用越來越廣泛。

1.2 改進(jìn)條件生成對(duì)抗網(wǎng)絡(luò)

pix2pix 是cGAN 的變體之一，它可以進(jìn)行圖像對(duì)圖像的風(fēng)格變換，由于該網(wǎng)絡(luò)有強(qiáng)大的變現(xiàn)能力，使其在圖像修復(fù)中得到了廣泛應(yīng)用。pix2pix的生成器是跳躍連接形成的U 型結(jié)構(gòu)，維度從高到低又從低到高，且左右兩邊對(duì)應(yīng)維度通道相拼接。跳躍連接可以有效保留輸入數(shù)據(jù)的潛在信息，因此，對(duì)于局部切片圖像的局部特征信息，可以得到很好的保留和傳遞，這有利于局部切片圖像的完整重構(gòu)。

為了更好地獲得圖像中的細(xì)節(jié)，網(wǎng)絡(luò)需要搭建到一定的深度才能重構(gòu)出優(yōu)質(zhì)的圖像。然而，一味地增加網(wǎng)絡(luò)的深度會(huì)造成深度模型退化，這樣的網(wǎng)絡(luò)是很難訓(xùn)練的。本文在pix2pix 網(wǎng)絡(luò)的生成器上加入嵌套殘差密集塊RRDB 來減弱網(wǎng)絡(luò)模型的退化，RRDB 的殘差塊里有殘差組合的密集塊，因此其擁有殘差塊[22]和密集塊[23]各自的優(yōu)勢(shì)。

另外，極化自注意力(PSA)模塊是能耗盡通道注意力和空間注意力表現(xiàn)能力的自注意力模塊，將其應(yīng)用到網(wǎng)絡(luò)中，可以幫助提高網(wǎng)絡(luò)的辨別能力，因此，本文把PSA 模塊加入到pix2pix 網(wǎng)絡(luò)的判別器中。

最后，將網(wǎng)絡(luò)中的各個(gè)模塊進(jìn)行搭建和調(diào)試，得到了改進(jìn)條件生成對(duì)抗網(wǎng)絡(luò)(RRDB-PSA-pix2pix)。圖1 為該網(wǎng)絡(luò)的生成器和判別器的架構(gòu)，其中一個(gè)方塊代表一個(gè)網(wǎng)絡(luò)塊，塊中傳遞順序從左往右，一個(gè)字母表示一個(gè)網(wǎng)絡(luò)層，如C 為卷積層，L 為L(zhǎng)eakyReLU 激活函數(shù)層，R 為RRDB 模塊層，B為批歸一化層，U 為上采樣層，D 為Dropout 層，Re 為ReLU 激活函數(shù)層，T 為Tanh 激活函數(shù)層，P 為PSA 注意力層，F(xiàn) 為Flatten 層。

如圖1a 所示，生成器的左半邊(編碼作用)的主線輸出為out1，跳躍連接上輸出的是out2。生成器的右半邊(解碼作用)整合局部圖像信息和左半邊傳遞過來的out2 信息，其RRDB 不輸出out2，只輸出out1。另外，RRDB 中的殘差尺度參數(shù)在密集模塊的支線上，這會(huì)阻礙密集模塊數(shù)據(jù)的傳遞，并使局部切片圖像中的噪聲信息傳遞到下層，進(jìn)而影響生成圖像的質(zhì)量。因此，本文將RRDB 架構(gòu)進(jìn)行了改進(jìn)，即將殘差尺度參數(shù)放置在“短路連接”上，這樣局部切片圖像的信息可以得到保留，并且密集模塊可以學(xué)習(xí)局部切片圖像中剩余的細(xì)節(jié)。同時(shí)，考慮到資源分配的問題，改進(jìn)后的架構(gòu)使用2 個(gè)密集殘差和2 層密集模塊。如圖2 所示，改進(jìn)的RRDB 架構(gòu)有兩個(gè)輸出，分別為out1 和out2，out1 比out2 包含更多的局部圖像信息。

圖1 RRDB-PSA-pix2pix 網(wǎng)絡(luò)架構(gòu)

圖2 本文改進(jìn)的RRDB 架構(gòu)

對(duì)于條件生成對(duì)抗網(wǎng)絡(luò)，判別器損失采用Wasserstein 距離[24]和梯度懲罰[25]進(jìn)行目標(biāo)優(yōu)化，同時(shí)需要加上其條件概率，可以表示為：

同時(shí)，為了讓生成器生成更高質(zhì)量的圖像，在生成器損失中加入感知損失，計(jì)算過程表示為：

式中， Φ (G(x))表示利用VGG19 網(wǎng)絡(luò)[26]從生成圖像中提取的特征映射(不通過最后激活層)； Φ (x)表示從真實(shí)圖像中提取的特征。

孔隙率是多孔結(jié)構(gòu)最重要的結(jié)構(gòu)參數(shù)。為了得到與真實(shí)圖像孔隙率相似的生成圖像，將生成圖像孔隙率和真實(shí)圖像孔隙率的曼哈頓距離作為孔隙率損失，如下所示：

式中， φporosity為 孔隙率；Avoid為孔隙(即空相)的面積；Avoid+Asolid為總面積。

另外，為了更準(zhǔn)確地描述生成圖像和真實(shí)圖像的相似情況，本文在生成器損失中加入L1損失。根據(jù)式(1)～式(3)，可以得到網(wǎng)絡(luò)的生成器損失表達(dá)：

式中，y是約束條件；L1損失表示真實(shí)圖像數(shù)據(jù)分布x與 生成器生成圖像分布G(x)的L1范數(shù)，可寫為：

2 局部切片圖像重構(gòu)

2.1 局部切片圖像的獲取

利用工業(yè)級(jí)CT 掃描機(jī)(型號(hào)：Y.CT modular,3.5 mA, 200 kV)對(duì)生物骨(豬肋骨)進(jìn)行掃描，得到40 張切片圖像。對(duì)切片圖像進(jìn)行二值化處理和形態(tài)處理，得到更高質(zhì)量的圖像。然后對(duì)處理過的切片圖像進(jìn)行裁剪，每張切片圖像可裁剪出4～5 張骨多孔圖像?？紤]到圖像質(zhì)量和計(jì)算成本，本文對(duì)裁剪后的骨多孔圖像進(jìn)行縮放，縮放至512×512 大小，則得到190 張圖像，如圖3 所示。實(shí)際上，190 張圖像作為深度學(xué)習(xí)的數(shù)據(jù)集是遠(yuǎn)遠(yuǎn)不夠的，本文通過對(duì)骨多孔圖像進(jìn)行旋轉(zhuǎn)和翻轉(zhuǎn)來實(shí)現(xiàn)增強(qiáng)數(shù)據(jù)的目的，最終得到1 140 張圖像作為數(shù)據(jù)集。

圖3 數(shù)據(jù)集獲取過程

將上述獲得的骨多孔圖像與掩膜矩陣相乘得到局部切片圖像，如下所示：

式中，f(x,y)為 骨多孔圖像；f(x′,y′)為局部切片圖像；掩膜矩陣M中缺失區(qū)域信息為1，背景區(qū)域?yàn)?，如圖4 所示。

圖4 局部骨切片圖像的生成

2.2 網(wǎng)絡(luò)的輸入

本文采用在掩膜矩陣的背景圖像上加隨機(jī)噪聲作為網(wǎng)絡(luò)生成器的輸入，如圖5 所示。

圖5 生成器輸入形式

將局部切片圖像設(shè)置為網(wǎng)絡(luò)的約束條件，同時(shí)，對(duì)局部切片圖像進(jìn)行重構(gòu)，需要構(gòu)建局部切片圖像與重構(gòu)生成圖像之間的關(guān)系，否則當(dāng)生成器生成出與輸入不對(duì)應(yīng)但又清晰的圖像時(shí)，判別器仍然會(huì)給高分，影響重構(gòu)結(jié)果。將不帶噪聲的局部切片圖像與重構(gòu)生成圖像拼接在一起，輸入到判別器中進(jìn)行判斷，如圖6 所示。

圖6 判別器輸入形式

2.3 網(wǎng)絡(luò)模型訓(xùn)練

本文提出的網(wǎng)絡(luò)模型在搭載單張RTX 2070super GPU 及Intel-i5 10600KF 主機(jī)的前提下，使用pytorch框架進(jìn)行訓(xùn)練和推理，訓(xùn)練過程采用如下超參數(shù)。

對(duì)于LeakyReLU 激活函數(shù)，參數(shù)a設(shè)置為0.2，式(5)中 λG設(shè) 為0.1， λperceptual設(shè) 為1， λL1設(shè)為0.02， λporosity取0.3；對(duì) 于 所 有 模 型，均 使 用Adam 優(yōu)化器(β1=0， β2=0.9)進(jìn)行訓(xùn)練；使用雙時(shí)間尺度更新規(guī)則，即判別器的學(xué)習(xí)速率設(shè)為0.000 4，生成器的學(xué)習(xí)速率設(shè)為0.000 1；式(1)中判別器梯度懲罰的參數(shù)λ 設(shè)為10。

骨多孔圖像數(shù)據(jù)集被按照3:1 比例劃分為訓(xùn)練集和測(cè)試集。將訓(xùn)練集中的圖像數(shù)據(jù)進(jìn)行處理，得到局部切片圖像數(shù)據(jù)集。為了完整地描述骨多孔結(jié)構(gòu)的特征，本文以閔可夫斯基泛函(minkowski functionals, MF)[27]作為測(cè)試集測(cè)試網(wǎng)絡(luò)的判據(jù)。對(duì)于二維骨切片圖像，它們的 M F 由 總表面積S、總周長(zhǎng)C和2D 歐拉示性數(shù)組成，表示為 MF=[S,C,χ]。同時(shí)，設(shè)置測(cè)試集T={t1,t2,···,tN}， 生成圖像G(T)={G(t1),G(t2),···,G(tN)}， 判據(jù) ?MF可以表示為式(8)，其值越小說明網(wǎng)絡(luò)越好：

網(wǎng)絡(luò)訓(xùn)練的周期為500 輪，每一輪訓(xùn)練后的生成器模型對(duì)測(cè)試集圖像進(jìn)行重構(gòu)，用閔可夫斯基泛函判據(jù)來確定最好的模型，整個(gè)訓(xùn)練過程如圖7 所示。

圖7 RRDB-PSA-pix2pix 網(wǎng)絡(luò)訓(xùn)練過程

2.4 不同網(wǎng)絡(luò)的對(duì)比分析

近年來，文獻(xiàn)[7-8]提出U-NetGAN 網(wǎng)絡(luò)和Dense-W-pix2pix 網(wǎng)絡(luò)對(duì)多孔結(jié)構(gòu)圖像進(jìn)行重構(gòu)。本文分別對(duì)U-NetGAN 網(wǎng)絡(luò)和Dense-W-pix2pix 網(wǎng)絡(luò)進(jìn)行修改使其能夠訓(xùn)練本文數(shù)據(jù)集，并使用MF 判據(jù)，將二者與本文提出的RRDB-PSA-pix2pix網(wǎng)絡(luò)進(jìn)行對(duì)比。骨多孔圖像(圖8a)與掩膜矩陣相乘得到了如圖8b 所示的局部骨切片圖像。3 種網(wǎng)絡(luò)分別對(duì)局部骨切片圖像進(jìn)行重構(gòu)，得到了圖8c～8e 所示的生成圖像。

圖8 3 種網(wǎng)絡(luò)的生成圖像比較

如圖8 所示，U-NetGAN 網(wǎng)絡(luò)生成的圖像質(zhì)量最差，圖中噪聲多，孔隙邊緣輪廓不清晰，出現(xiàn)偽影，且生成的圖像輪廓與真實(shí)圖像相同。而另外兩種網(wǎng)絡(luò)生成的圖像表現(xiàn)較好。同時(shí)，由表1 可以看出， RRDB-PSA-pix2pix 網(wǎng)絡(luò) ?MF判據(jù)值比Dense-W-pix2pix 網(wǎng)絡(luò)更小，由此可知，本文提出的RRDB-PSA-pix2pix 網(wǎng)絡(luò)重構(gòu)出的圖像質(zhì)量最好。

表1 3 種網(wǎng)絡(luò)MF 判據(jù)比較

2.5 重構(gòu)圖像與真實(shí)骨多孔圖像相似性分析

本文通過3 種形態(tài)學(xué)函數(shù)分析和局部孔隙率分布研究[28]，評(píng)價(jià)了RRDB-PSA-pix2pix 網(wǎng)絡(luò)重構(gòu)生成的圖像與真實(shí)圖像的相似程度。這3 種形態(tài)學(xué)函數(shù)分別為兩點(diǎn)相關(guān)函數(shù)[29]、線性路徑函數(shù)[30]和弦長(zhǎng)分布函數(shù)[31]。兩點(diǎn)統(tǒng)計(jì)函數(shù)S2(r)用來表示空間中兩點(diǎn)的關(guān)系。線性路徑函數(shù)L(r) 統(tǒng) 計(jì)給定長(zhǎng)度為r的線段完全位于同一相的概率，包含了多孔結(jié)構(gòu)的連通性信息。弦長(zhǎng)分布函數(shù)C LD(r)計(jì)算位于同一相位相同弦長(zhǎng)的概率，可以較為準(zhǔn)確地捕獲多孔圖像內(nèi)部的空間信息，用來描述孔洞大小、形狀和空間排布等特征。局部孔隙率分布通常用于表征數(shù)字化模型中不同長(zhǎng)度下的孔隙率和連通性的波動(dòng)。與孔隙率不同，局部孔隙率? (r,L) 是在位于r為中心且長(zhǎng)度為L(zhǎng)的二維測(cè)量單元K(r,L)內(nèi)測(cè)得。圖9 為3 種形態(tài)學(xué)函數(shù)和局部孔隙率分布的定義示意圖。

圖9 3 種形態(tài)學(xué)函數(shù)和局部孔隙率分布的定義示意圖

本文以Dense-W-pix2pix 網(wǎng)絡(luò)重構(gòu)生成的圖像作為對(duì)比，來評(píng)價(jià)分析RRDB-PSA-pix2pix 網(wǎng)絡(luò)生成的圖像與真實(shí)圖像的相似程度。圖10 為兩種網(wǎng)絡(luò)生成圖像的形態(tài)學(xué)函數(shù)和局部孔隙率分布統(tǒng)計(jì)圖，其中局部孔隙率分布中的測(cè)量單元取大小為320×320。

圖中的統(tǒng)計(jì)分布均進(jìn)行了兩組獨(dú)立樣本的Kolmogorov-Smirnov 檢驗(yàn)(K-S 檢驗(yàn))，且p值皆大于0.05。由圖10 可以看出， RRDB-PSA-pix2pix網(wǎng)絡(luò)生成的圖像分布曲線更接近于真實(shí)圖像的分布曲線，因此可以得出，本文提出的網(wǎng)絡(luò)能夠重構(gòu)出高質(zhì)量的完整骨多孔切片圖像，并且其兩相位置分布關(guān)系和連通性更符合真實(shí)情況。

圖10 形態(tài)學(xué)函數(shù)與局部孔隙率分布統(tǒng)計(jì)圖

3 結(jié) 束 語

孔隙率是骨多孔結(jié)構(gòu)最重要的結(jié)構(gòu)參數(shù)之一，孔隙率的不同，直接影響著骨多孔結(jié)構(gòu)的質(zhì)量、滲透性能、機(jī)械性能以及相關(guān)生物性能，實(shí)際上，除孔隙率外，還有一些參數(shù)可用來評(píng)價(jià)骨多孔結(jié)構(gòu)，但是，不同參數(shù)的評(píng)價(jià)方法，所采用的網(wǎng)絡(luò)算法與數(shù)學(xué)模型也各不相同，因此，考慮到篇幅問題，本文重點(diǎn)考慮通過孔隙率來評(píng)價(jià)骨多孔結(jié)構(gòu)。

1) 本文在pix2pix 基礎(chǔ)上，對(duì)生成器進(jìn)行改進(jìn)，加入嵌套殘差密集塊(RRDB)，在判別器中加入PSA 模塊，提出對(duì)局部骨切片圖像進(jìn)行完整修復(fù)重構(gòu)的RRDB-PSA-pix2pix 網(wǎng)絡(luò)。

2) 研究并分析了閔可夫斯基泛函對(duì)骨多孔結(jié)構(gòu)唯一性的表征，建立了閔可夫斯基泛函判據(jù)，為快速篩選網(wǎng)絡(luò)提供了方法。

3) 利用3 種形態(tài)學(xué)函數(shù)與局部孔隙率分布對(duì)本文提出的RRDB-PSA-pix2pix 網(wǎng)絡(luò)進(jìn)行分析評(píng)價(jià)。研究表明，該網(wǎng)絡(luò)重構(gòu)生成的骨切片圖像中的多孔結(jié)構(gòu)特征，更加接近于真實(shí)圖像，并且生成的圖像質(zhì)量更高。

4) 本文所提出的網(wǎng)絡(luò)框架可以很容易地?cái)U(kuò)展到其他應(yīng)用，如2D-to-3D 或3D-to-3D 的重構(gòu)。理論上，它可以將任意數(shù)量的任何類型的對(duì)象函數(shù)以及任何用戶定義的條件數(shù)據(jù)合并到重構(gòu)中。當(dāng)通知某一區(qū)域存在特定結(jié)構(gòu)時(shí)，這在實(shí)踐中可能特別有用。此外，它還可以通過與其他方法相結(jié)合來降低計(jì)算成本。