基于可控卷積曲面的三維神經元建模

朱曉強，陳 琦

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

0 概述

自現代神經科學誕生以來，深入認識神經元形態是神經科學領域重要的研究方向之一。大腦的復雜結構使得密集神經網絡特征和行為過程的分析變得極為困難。三維可視化已被證明是評估復雜神經系統的有效工具。

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

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

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

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

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

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

1 神經元數據預處理

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

表1 標準神經元數據結構Table 1 Structure of standard neuron data

原始神經元數據存在非常短的神經元骨架節段以及不規則的頂點采樣密度。此外，數據還可能存在潛在異常采樣點和分支，如果直接構建網格模型會導致網格不連續、拓撲異常和分支重疊，因此，需要對神經元數據進行預處理。

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

圖1 神經元節段線骨架示意圖Fig.1 Schematic diagram of neuron segmental line skeleton

由于短邊現象幾乎存在于所有的神經元形態數據中，因此針對采樣點屬于不同類型的拓撲角色，提出以下預處理算法流程：

輸入SWC 格式神經元數據

輸出去除異常采樣點和短邊后的神經元數據

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

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

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

步驟4如果PA是內部點，PB是分支結束點，處理方法同步驟3。

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

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

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

步驟8根據集合U，刪除單個分離的神經元分支。

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

圖2 三維神經元預處理結果Fig.2 Preprocessing result of 3D neuron

2 三維神經元建模

2.1 卷積曲面

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

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

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

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

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

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

2.2 卷積曲面可控混合

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

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

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

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

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

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

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

2.3 稀疏體素結構

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

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

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

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

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

2.4 網格后處理

2.4.1 網格重構

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

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

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

2.4.2 網格細分

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

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

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

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

每條邊的中點ec 為：

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

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

之后，更新原始點，設v為原始點，則更新后的點v'為：

2.4.3 網格體素化

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

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

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

3 實驗結果與分析

3.1 建模結果分析

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

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

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

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

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

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

圖9 基于稀疏體素結構的卷積曲面Fig.9 Convolution surface based on sparse voxels structure

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

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

3.2 網格后處理分析

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

圖10 三維神經元網格模型優化結果Fig.10 Optimization results of 3D neuron grid model

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

圖11 三維神經元網格體素化Fig.11 Grid voxelization of 3D neuron

4 結束語

針對生成表面網格的光滑度和網格分支處的平滑度問題，結合卷積曲面技術與線性插值思想，本文提出可控卷積曲面混合方法，用于構造三維神經元模型。采用稀疏體素技術加快網格生成速度，解決模型生成效率的問題。與現有神經元建模方法相比，本文方法具有較優的魯棒性和高效性。下一步將研究復雜神經元骨架數據產生的環狀結構問題，實現高質量神經元數據的建模。