基于可控卷積曲面的三維神經(jīng)元建模

朱曉強，陳 琦

（上海大學 通信與信息工程學院，上海 200444）

0 概述

自現(xiàn)代神經(jīng)科學誕生以來，深入認識神經(jīng)元形態(tài)是神經(jīng)科學領(lǐng)域重要的研究方向之一。大腦的復(fù)雜結(jié)構(gòu)使得密集神經(jīng)網(wǎng)絡(luò)特征和行為過程的分析變得極為困難。三維可視化已被證明是評估復(fù)雜神經(jīng)系統(tǒng)的有效工具。

GLASER 等［1］提出Neurolucida 軟件應(yīng)用，使用錐形圓柱體表示神經(jīng)元節(jié)段，在分支點處通過球體將神經(jīng)元節(jié)段相連接，生成的網(wǎng)格模型在段之間是不連續(xù)的，不滿足水密性的特點。NeuGen［2］、GENESIS［3］等軟件應(yīng)用可以生成神經(jīng)元模型表面，但均使用球體表示胞體，生成的體細胞網(wǎng)格質(zhì)量較差。AGUIAR 等［4］開發(fā)Py3DN 工具箱，通過幾何近似的方法實現(xiàn)更真實的軀體重建，但沒有將生成的細胞體與樹突相連接。LASSERRE 等［5］研究用于電生理仿真的皮層神經(jīng)元的重建，提出從一組采樣神經(jīng)元形態(tài)數(shù)據(jù)中生成一個連續(xù)的表示神經(jīng)元膜表面的算法，但代碼難以復(fù)現(xiàn)，因分支上存在尖銳的邊緣，導(dǎo)致重建網(wǎng)格的質(zhì)量較差。

BRITO 等［6］提出使用Neuronize 工具構(gòu)建神經(jīng)元細胞三維模型，主要使用胡克定律和質(zhì)量彈簧模型在合理的物理基礎(chǔ)上重建神經(jīng)元的胞體模型，但與GLASER［1］相似都采用圓柱體表示神經(jīng)元節(jié)段，改進之處是使用剪切和拼接的方式連接神經(jīng)元段與段、段與胞體之間的縫隙。GARCIA-CANTERO等［7］在Neuronize 的基礎(chǔ)上提出NeuroTessMesh 硬件加速工具，能夠解決復(fù)雜場景下存儲和計算成本較高的問題，以生成具有自適應(yīng)動態(tài)細化功能的神經(jīng)元網(wǎng)格模型，但是該工具生成的網(wǎng)格在分支處不平滑且不滿足水密性要求。

ABDELLAH 等［8］提出從新皮層神經(jīng)元電路的形態(tài)骨架重建規(guī)模較大且高度詳細體積模型的新方法。通過對單個神經(jīng)元形態(tài)數(shù)據(jù)進行預(yù)處理，從預(yù)處理后的神經(jīng)元形態(tài)骨架創(chuàng)建分段平滑和符合水密性要求的網(wǎng)格模型，然后構(gòu)造局部體積模型，最后合并成全局體積模型，其缺點在于網(wǎng)格是分段水密的。

文獻［9-11］提出用于反映擴散模擬的四面體體積模型，其網(wǎng)格滿足水密性和雙流形的特點。由于該體積模型具有較高的鑲嵌性且忽略神經(jīng)元形態(tài)的真實有機外觀，因此不能用于交互式可視化分析。

ABDELLAH 等［12］提出使用Skin Modifiers 從神經(jīng)元的形態(tài)骨架生成高保真的新皮層神經(jīng)元表面網(wǎng)格。KARLSSON 等［13］基于符號距離函數(shù)和光線追蹤［14］技術(shù)，提出用于渲染大規(guī)模神經(jīng)元電路數(shù)據(jù)的高保真可視化方法，但是生成網(wǎng)格的質(zhì)量不高且不平滑。

為生成高質(zhì)量、平滑且滿足水密性要求的三維神經(jīng)元網(wǎng)格模型，本文提出一種三維神經(jīng)元建模方法。通過預(yù)處理NeuroMorpho.Org 數(shù)據(jù)庫提供的神經(jīng)元形態(tài)數(shù)據(jù)，對現(xiàn)有基于骨架的卷積曲面混合建模方法進行改進，設(shè)計一種基于可變支撐半徑控制的卷積曲面混合新方法。利用稀疏體素結(jié)構(gòu)提高分辨率并加快等值面提取過程，以優(yōu)化網(wǎng)格生成的神經(jīng)元三角網(wǎng)格模型和體積模型。

1 神經(jīng)元數(shù)據(jù)預(yù)處理

從NeuroMorpho.Org 數(shù)據(jù)庫中獲取的SWC 格式神經(jīng)元形態(tài)數(shù)據(jù)以ASCII 形式進行編碼，每個非空行代表一個神經(jīng)元樣本點，其中包含七個數(shù)據(jù)項，具體如表1 所示。其中：采樣號和父采樣號表示空間中樣本點的父子連接關(guān)系；結(jié)構(gòu)標識符主要包括0-未定義、1-胞體、2-軸突、3-樹突和4-尖端樹突5 種類型；X、Y、Z是以μm 為單位的空間坐標；半徑是樹突厚度的1/2，單位是μm。

表1 標準神經(jīng)元數(shù)據(jù)結(jié)構(gòu)Table 1 Structure of standard neuron data

原始神經(jīng)元數(shù)據(jù)存在非常短的神經(jīng)元骨架節(jié)段以及不規(guī)則的頂點采樣密度。此外，數(shù)據(jù)還可能存在潛在異常采樣點和分支，如果直接構(gòu)建網(wǎng)格模型會導(dǎo)致網(wǎng)格不連續(xù)、拓撲異常和分支重疊，因此，需要對神經(jīng)元數(shù)據(jù)進行預(yù)處理。

神經(jīng)元節(jié)段線骨架示意圖如圖1 所示。對于線骨架Si，PA、PB分別是左右兩端采樣點，rA、rB分別是采樣點的半徑，lAB是線骨架的長度，引入可變參數(shù)e，如果lAB≤e×Max(rA，rB)，其中e≥2，則定義Si為短邊。

圖1 神經(jīng)元節(jié)段線骨架示意圖Fig.1 Schematic diagram of neuron segmental line skeleton

由于短邊現(xiàn)象幾乎存在于所有的神經(jīng)元形態(tài)數(shù)據(jù)中，因此針對采樣點屬于不同類型的拓撲角色，提出以下預(yù)處理算法流程：

輸入SWC 格式神經(jīng)元數(shù)據(jù)

輸出去除異常采樣點和短邊后的神經(jīng)元數(shù)據(jù)

步驟1定義U為已訪問過的節(jié)點集合，依次遍歷節(jié)點，刪除節(jié)點半徑值為負數(shù)的樣本點，并連接前后節(jié)點。

步驟2移除長度接近0（小于1e-4）的分支。

步驟3處理短邊問題，如果PA、PB都是內(nèi)部點，首先獲取PA的父節(jié)點P0，然后連接P0和PB，最后刪除PA節(jié)點。

步驟4如果PA是內(nèi)部點，PB是分支結(jié)束點，處理方法同步驟3。

步驟5如果PA是分支點，PB是內(nèi)部點，由于刪除分支點會破壞整體拓撲結(jié)構(gòu)，則首先獲取PB的子節(jié)點PC，然后連接PA和PC，最后刪除PB節(jié)點。

步驟6如果PA是分支點，PB是分支結(jié)束點，直接刪除PB。

步驟7如果PA是根節(jié)點，PB是其他節(jié)點，處理方法同步驟5。

步驟8根據(jù)集合U，刪除單個分離的神經(jīng)元分支。

以人類大腦區(qū)域的神經(jīng)元電路數(shù)據(jù)為例，三維神經(jīng)元預(yù)處理結(jié)果如圖2 所示。原始神經(jīng)元數(shù)據(jù)包含1 429 個采樣點，預(yù)處理后的神經(jīng)元數(shù)據(jù)包含1 374 個采樣點。預(yù)處理算法能有效地解決神經(jīng)元形態(tài)數(shù)據(jù)可能存在的異常樣本點、短邊等問題，減少采樣點數(shù)量，為后續(xù)網(wǎng)格模型的生成提供有效的輸入數(shù)據(jù)。

圖2 三維神經(jīng)元預(yù)處理結(jié)果Fig.2 Preprocessing result of 3D neuron

2 三維神經(jīng)元建模

2.1 卷積曲面

卷積曲面最初被定義為由幾何骨架和高斯核函數(shù)的卷積生成的［15］，其數(shù)學定義如式（1）所示：

其中：λi為第i段幾何骨架的權(quán)重值；Fi(x，y，z)為第i段幾何骨架的勢能函數(shù)；參數(shù)T為提取等值面的閾值。

根據(jù)式（1），假設(shè)p為三維空間中的點，S為幾何骨架，骨架的幾何函數(shù)可定義為：

對于一個核函數(shù)f，且三維空間點q∈S，幾何骨架S在點p處的總勢能值可定義為：

核函數(shù)按作用范圍分為無限支撐和有限支撐核函數(shù)。對于無限支撐核函數(shù)，定義域內(nèi)所有的勢能值均大于零，例如高斯函數(shù)［15］和柯西函數(shù)［16］。對于有限支撐核函數(shù)，在指定的定義域內(nèi)有值，在其余范圍內(nèi)的勢能值均為零，例如Metaballs［17］。采用無限支撐核函數(shù)構(gòu)造卷積曲面有兩個缺點：1）在構(gòu)造單個神經(jīng)元骨架模型時，因卷積場計算過遠導(dǎo)致生成曲面效率低；2）當曲面混合時，一旦勢能函數(shù)發(fā)生變化，須重新計算整個骨架的勢能場，無法做到局部控制。鑒于此，本文參考JIN 等［18］提出的方法使用四次多項式作為卷積核函數(shù)，如式（4）所示：

其中：R為核函數(shù)的有效支撐半徑，用于剪切離目標太遠的骨架，可有效控制混合范圍。

2.2 卷積曲面可控混合

基于ZHU 等［19］提出的結(jié)論，對于一條線骨架Si，如果曲面厚度為di，則p點處的權(quán)重值λi滿足如下關(guān)系：

為方便控制局部混合，本文引入比例參數(shù)t=Ri/di，其中，Ri為支撐半徑。對于常規(guī)卷積曲面，di和Ri的默認關(guān)系為di=Ri/2，即t=2。根據(jù)式（5）求得λi為：

通過改變參數(shù)t的值來改變空間每個位置相對骨架的權(quán)重值，從而改變基于骨架的卷積曲面混合程度。單純改變t雖然可調(diào)節(jié)混合程度，但是無法做到局部混合。本文根據(jù)線性插值的思想來局部調(diào)節(jié)線骨架上某一段的支撐半徑值，提出一種基于可變支撐半徑控制的卷積曲面混合新方法。基于可變支撐半徑變化的卷積曲面示意圖如圖3 所示。

圖3 基于可變支撐半徑變化的卷積曲面示意圖Fig.3 Schematic diagram of convolution surface based on variable support radius change

假設(shè)線骨架Si的左端點C1處的有限支撐半徑為o1，右端點C2處的有限支撐半徑為o2，l10為C1C0的長度，l12為C1C2的長度，則在線骨架之間任一點C0處的有限支撐半徑可定義為：

當線骨架混合時，通過為每一段線骨架設(shè)置不同的左右端點支撐半徑，實現(xiàn)局部混合的目的。

因卷積曲面的疊加性，在線骨架段之間自動混合生成光滑連續(xù)的表面，但是在多段線骨架之間會因勢能值的不斷疊加而產(chǎn)生額外的“鼓包”現(xiàn)象。利用上述方法可有效地緩解“鼓包”問題以及解決常規(guī)卷積曲面因支撐半徑范圍過大而產(chǎn)生的過渡混合問題。

2.3 稀疏體素結(jié)構(gòu)

在上節(jié)中，本文將三維神經(jīng)元骨架作為輸入數(shù)據(jù)，并基于卷積曲面生成三維空間離散勢能場，利用經(jīng)典的Marching Cubes 算法［20］提取等值面，但在實際計算中可能會產(chǎn)生問題。基于Marching Cubes 算法的等值面提取示意圖如圖4 所示。假設(shè)在每個體素邊長為1 的多個立方體中計算一條曲面厚度為1，支撐半徑為1.2 的線骨架的卷積曲面，當使用Marching Cubes 算法提取等值面時，等值面和支撐半徑可能會處于同一個體素中，導(dǎo)致最終無法有效提取高質(zhì)量的神經(jīng)元網(wǎng)格模型。

圖4 基于Marching Cubes 算法的等值面提取示意圖Fig.4 Schematic diagram of isosurface extraction based on Marching Cubes algorithm

通過提高體素分辨率的方法，使同一個體素內(nèi)只包含等值面和支撐半徑就可解決上述問題，但密集網(wǎng)格在構(gòu)造多個模型實例時，其內(nèi)存占用隨模型體積的增大而增大。由于三維神經(jīng)元復(fù)雜的分支結(jié)構(gòu)，空間中存在大量的空隙，因此不需要在一個密集的網(wǎng)格中均勻采樣數(shù)據(jù)，可采用稀疏的體積數(shù)據(jù)結(jié)構(gòu)來提高效率。

本文參考MUSETH［21］提出的基于VDB 的稀疏體素結(jié)構(gòu)，該結(jié)構(gòu)具有內(nèi)存效率高、支持時變數(shù)據(jù)模擬、自適應(yīng)采樣等特點。核心思想是在一個類似B+樹的多層次數(shù)據(jù)結(jié)構(gòu)中動態(tài)地安排網(wǎng)格處的體素塊，將空間分為根節(jié)點、兩個內(nèi)部節(jié)點和一個葉子節(jié)點三部分。其中，最小的體元素按活躍值可分為活躍體素和非活躍體素，且每一個活躍體素都有一個區(qū)別于背景值的默認值。

針對分辨率較低的問題，本文通過改變體素大小來解決，設(shè)置默認背景值為0，活躍體素值為0.5，利用上述稀疏體素結(jié)構(gòu)計算三維空間卷積勢能場，以加快等值面提取速度。

2.4 網(wǎng)格后處理

2.4.1 網(wǎng)格重構(gòu)

基于稀疏體素結(jié)構(gòu)生成的神經(jīng)元模型網(wǎng)格頂點分布不均勻，會影響基于網(wǎng)格的數(shù)值模擬實驗等應(yīng)用。因此，本文需要重構(gòu)并均勻化網(wǎng)格且均勻的網(wǎng)格必須滿足三個基本條件：所有邊等長、所有三角形面積相等和所有頂點度數(shù)為60°。

本文參考BOTSCH等［22］提出的Isotropic Remeshing網(wǎng)格優(yōu)化算法，設(shè)定一個目標邊長L，遍歷網(wǎng)格中所有邊長向目標邊長靠近，得到最終優(yōu)化結(jié)果。具體處理步驟如下：1）分割所有長度大于4L/3 的邊；2）坍縮所有長度小于4L/5 的邊；3）翻轉(zhuǎn)所有邊以優(yōu)化度數(shù)；4）將所有新加入的點映射到切平面。

為有效地生成體網(wǎng)格數(shù)據(jù)，本文利用MeshLab［23］工具驗證三維神經(jīng)元網(wǎng)格模型的水密性，使其不存在非流形邊和頂點。

2.4.2 網(wǎng)格細分

為獲得表面更加光滑且包含更多細節(jié)信息的三維神經(jīng)元網(wǎng)格模型，本文采用Loop［24］算法對重構(gòu)后的網(wǎng)格做進一步細分處理。細分過程主要包括計算邊點和更新原始點兩個過程。

在計算邊點的 過程中，假設(shè)四個點V0、V1、V2、V3，其拓撲關(guān)系如圖5 所示。

圖5 計算邊點的示意圖Fig.5 Schematic diagram of calculating edge point

針對由屬于兩個三角面的每條邊所組成的四點面，面點fp 為：

每條邊的中點ec 為：

fp 與ec 點的平均點，即邊點ep 為：

其中：eep表示該邊屬于兩個三角面片，即邊點處于網(wǎng)格內(nèi)部。當邊點位于網(wǎng)格邊界時，邊點直接退化為中點。

之后，更新原始點，設(shè)v為原始點，則更新后的點v'為：

2.4.3 網(wǎng)格體素化

為生成能在腦神經(jīng)元組織中進行光傳播模擬實驗的神經(jīng)元模型，本文利用TetGen［25］軟件生成高質(zhì)量的神經(jīng)元四面體體積模型。TetGen 體素化過程如圖6 所示。

圖6 TetGen 體素化流程Fig.6 Procedure of TetGen voxelization

首先，將細化后的三維神經(jīng)元網(wǎng)格作為輸入數(shù)據(jù)；然后，類比二維Delaunay 三角剖分，在三維空間中對神經(jīng)元網(wǎng)格模型進行四面體剖分，生成保持原有神經(jīng)元模型形態(tài)的約束網(wǎng)格；最后，通過限制插入模型內(nèi)頂點個數(shù)以及四面體體積來生成高質(zhì)量的三維神經(jīng)元體積模型。

3 實驗結(jié)果與分析

3.1 建模結(jié)果分析

為驗證本文提出的基于可控卷積曲面的三維神經(jīng)元建模方法的有效性，將第1 節(jié)中的神經(jīng)元形態(tài)骨架作為輸入數(shù)據(jù)，并與ABDELLAH 等［8］、ABDELLAH 等［12］及KARLSSON 等［13］所提方法的建模結(jié)果進行對比。

不同方法的三維神經(jīng)元建模結(jié)果對比如圖7 所示。文獻［8］所提的方法生成的神經(jīng)元網(wǎng)格是分段水密的，在細胞體與分支之間會出現(xiàn)斷層，即整體網(wǎng)格不滿足水密性要求，且在分支處曲面不平滑，易產(chǎn)生折痕現(xiàn)象。文獻［12］所提的方法相較于文獻［8］生成的神經(jīng)元網(wǎng)格曲面更加平滑，但對于復(fù)雜神經(jīng)元骨架易生成不滿足整體水密性的光滑網(wǎng)格模型，且在多個骨架相交處易生成異常拓撲結(jié)構(gòu)。上述兩種方法的額外貢獻點在于利用網(wǎng)格變形技術(shù)構(gòu)造細胞體的結(jié)構(gòu)。文獻［13］利用平滑函數(shù)解決神經(jīng)元分支處的折痕現(xiàn)象，但產(chǎn)生額外的曲面鼓包現(xiàn)象。本文根據(jù)卷積曲面的疊加性質(zhì)直接混合生成神經(jīng)元結(jié)構(gòu)模型。相較于前三種方法，本文方法生成的網(wǎng)格模型整體呈光滑連續(xù)，能夠更加真實地描繪三維神經(jīng)元形態(tài)。

圖7 不同方法的三維神經(jīng)元建模結(jié)果對比Fig.7 3D neuron modeling results comparison among different methods

由于神經(jīng)元形態(tài)骨架極其復(fù)雜，存在骨架之間距離較近的情況，因此可能會出現(xiàn)曲面相交的問題，影響基于網(wǎng)格的數(shù)值模擬等應(yīng)用。不同方法的近距離骨架建模結(jié)果對比如圖8 所示。文獻［8，12-13］都沒有提出有效的方法解決曲面相交問題。本文利用2.2 節(jié)中提出的基于局部可變支撐半徑控制的卷積曲面混合新方法來控制曲面混合程度，將過渡混合處的骨架支撐半徑和曲面厚度的比例系數(shù)t設(shè)置為1.2～2.0 之間動態(tài)線性插值，有效解決曲面相交問題，結(jié)果如圖8（d）所示。

圖8 不同方法的近距離骨架建模結(jié)果對比Fig.8 Short-range skeleton modeling results comparison among different methods

在線骨架厚度較小的位置，卷積曲面無法提取等值面或生成的網(wǎng)格質(zhì)量較差，本文利用稀疏體素結(jié)構(gòu)改變體素值的大小來提高分辨率，結(jié)果如圖9所示。隨著體素值的減小，生成的曲面越光滑，相應(yīng)的時耗也逐漸增加。

圖9 基于稀疏體素結(jié)構(gòu)的卷積曲面Fig.9 Convolution surface based on sparse voxels structure

為進一步說明本文方法的高效性，在搭載AMD Ryzen7-5800H（3.20 GHz）處理器、16 GB 內(nèi)存、Ryzen 集成顯卡的win11 電腦設(shè)備和Visual Studio 2019 平臺上，不同方法構(gòu)建神經(jīng)元網(wǎng)格模型所需的時間對比結(jié)果如表2 所示。其中，Z1表示分辨率574×1 054×208 像素，Z2表示分辨率1 437×2 635×522 像素。在相同分辨率下，文獻［12］所提的方法解決了文獻［8］所提方法存在的分支處不平滑問題，但平均時耗增加約130%；文獻［13］所提的方法相較于文獻［12］，在解決分支折痕問題的同時將平均時耗減少了約50%，但是在分支處產(chǎn)生額外鼓包問題；本文方法不僅能有效平滑分支處曲面，而且提供參數(shù)來控制曲面混合程度，以避免出現(xiàn)異常拓撲，建模速度相較于文獻［13］方法提升了約1.5 倍。

表2 不同方法構(gòu)建三維神經(jīng)元模型的時耗對比Table 2 Time consumption of constructing 3D neuron model comparison among different methods

3.2 網(wǎng)格后處理分析

為獲得更高質(zhì)量的三維神經(jīng)元網(wǎng)格模型數(shù)據(jù)，本文進行一系列的網(wǎng)格后處理操作包括網(wǎng)格重構(gòu)、網(wǎng)格細化和網(wǎng)格體素化，結(jié)果如圖10 所示。初始網(wǎng)格模型局部如圖10（a）所示，重構(gòu)后的網(wǎng)格局部如圖10（b）所示，細化后的網(wǎng)格局部如圖10（c）所示，優(yōu)化后的網(wǎng)格頂點分布更加均勻且表面更加光滑。

圖10 三維神經(jīng)元網(wǎng)格模型優(yōu)化結(jié)果Fig.10 Optimization results of 3D neuron grid model

為獲得用于光學實驗的高質(zhì)量體積模型，本文使用MeshLab 軟件驗證三角形網(wǎng)格的水密性，將.ply格式的網(wǎng)格數(shù)據(jù)作為輸入，利用TetGen 軟件直接生成，結(jié)果如圖11 所示。本文方法在保持神經(jīng)元形態(tài)的前提下，在網(wǎng)格內(nèi)部生成了等體積密集的四面體網(wǎng)格。

圖11 三維神經(jīng)元網(wǎng)格體素化Fig.11 Grid voxelization of 3D neuron

4 結(jié)束語

針對生成表面網(wǎng)格的光滑度和網(wǎng)格分支處的平滑度問題，結(jié)合卷積曲面技術(shù)與線性插值思想，本文提出可控卷積曲面混合方法，用于構(gòu)造三維神經(jīng)元模型。采用稀疏體素技術(shù)加快網(wǎng)格生成速度，解決模型生成效率的問題。與現(xiàn)有神經(jīng)元建模方法相比，本文方法具有較優(yōu)的魯棒性和高效性。下一步將研究復(fù)雜神經(jīng)元骨架數(shù)據(jù)產(chǎn)生的環(huán)狀結(jié)構(gòu)問題，實現(xiàn)高質(zhì)量神經(jīng)元數(shù)據(jù)的建模。