基于深度学习方法实现SPECT放射性核素定量测量（二）

基于深度学习方法实现SPECT放射性核素定量测量（二）—GATE仿真

（1）GATE介绍

GATE（GEANT4 Application for Emission Tomography）由OpenGATE组织开发，是一个基于GEANT4的软件包。它封装了GEANT4库，是一个专用于核医学领域的免费开源、通用性、脚本化、模块化的仿真平台。GATE继承了GEANT4仿真工具包的优势，其中包括经过验证的物理模型、复杂几何的描述、强大的可视化功能和3D渲染工具等，并且具有发射断层成像特有的功能。GATE软件的架构可由图1所示的层级结构表示，以GEANT4为内核（GEANT4 kernel），由核心层（Core layer）、应用层（Application layer）以及用户层（User level）组成。

图1：GATE层级结构

GATE仿真步骤由图2所示，包括设计探测器几何结构、设置体模参数、设置物理过程、初始化、设置电子学参数、设置放射源、设置随机数以及定义数据输出格式、开始仿真。在仿真中针对各步骤编写GATE宏语言，通过可视化检查探测器结构，最后通过主文件调用各模块代码执行程序。其中设置探测器结构、物理过程、电子学参数和放射源设置是GATE仿真重点。本课题仿真过程基于GATE9.0版本（目前最新版本貌似为9.1），在linux系统下实现，并且采用多线程并行计算加快仿真速度，缩短仿真时长。

图2：GATE仿真过程

GATE具体使用说明可以参考官方说明：http://opengatecollaboration.org/

（2）设置探测器

仿真探测器设置以西门子公司Symbia T2型双探头SPECT/CT系统为原型，构建SPECT系统结构。使用SPECTHead系统作为探测器几何结构的模板，自上而下定义的几何结构分别为：World、SPECTHead、 Shielding、Collimator、Scintillator、Back Compartment，然后需要将Scintillator通过attach命令将其连接到CrystalSD层次，使Scintillator具备记录粒子与物质相互作用信息的功能。定义的SPECTHead系统结构如图3所示。设置探测器选择360度，为了加快仿真速度，可以使用四个探头，每个探头只需选择90度即可，如图4所示。各结构材料如表1所示。

图3：SPECT探测器简图1
图4：SPECT探测器简图2

表1

探测器设置的代码如下所示：

# V I S U A L I S A T I O N
#####
/vis/disable

# M A N D A T O R Y
#####
/gate/geometry/setMaterialDatabase GateMaterials.db
# G E O M E T R Y
#####

# World
# Define the world dimensions
##
/gate/world/geometry/setXLength 100 cm
/gate/world/geometry/setYLength 100 cm
/gate/world/geometry/setZLength 100 cm

# Scanner Head
# Create a new box representing the main head-volume
# SPECThead is the name of the predefined SPECT system
# Create the SPECT system, which will yield an Interfile output of projection data
##
/gate/world/daughters/name SPECThead
/gate/world/daughters/insert box
/gate/SPECThead/geometry/setXLength  7. cm
/gate/SPECThead/geometry/setYLength 21. cm
/gate/SPECThead/geometry/setZLength 30. cm
/gate/SPECThead/placement/setTranslation  0. 25. 0. cm
/gate/SPECThead/placement/setRotationAxis 0 0 1
/gate/SPECThead/placement/setRotationAngle 90 deg
/gate/SPECThead/repeaters/insert ring
/gate/SPECThead/ring/setRepeatNumber 4
/gate/SPECThead/setMaterial Air
/gate/SPECThead/moves/insert orbiting
#/gate/SPECThead/orbiting/setSpeed 2.8125 deg/s
/gate/SPECThead/orbiting/setSpeed 1.40625 deg/s
#/gate/SPECThead/orbiting/setSpeed 0.703125 deg/s
/gate/SPECThead/orbiting/setPoint1 0 0 0 cm
/gate/SPECThead/orbiting/setPoint2 0 0 1 cm
/gate/SPECThead/vis/forceWireframe

# Shielding
# Create the shielding volume
##
/gate/SPECThead/daughters/name shielding
/gate/SPECThead/daughters/insert box
/gate/shielding/geometry/setXLength  7. cm
/gate/shielding/geometry/setYLength 21. cm
/gate/shielding/geometry/setZLength 30. cm
/gate/shielding/placement/setTranslation  0. 0. 0. cm
/gate/shielding/setMaterial Lead
/gate/shielding/vis/setColor red
/gate/shielding/vis/forceWireframe


# Collimator
# Create a full volume defining the shape of the collimator
##
/gate/SPECThead/daughters/name collimator
/gate/SPECThead/daughters/insert box
/gate/collimator/geometry/setXLength 3. cm
/gate/collimator/geometry/setYLength 19. cm
/gate/collimator/geometry/setZLength 28. cm
/gate/collimator/placement/setTranslation  -2. 0. 0. cm
/gate/collimator/setMaterial Lead
/gate/collimator/vis/setColor red
/gate/collimator/vis/forceWireframe

# Insert the first hole of air in the collimator
##
/gate/collimator/daughters/name hole
/gate/collimator/daughters/insert hexagone
/gate/hole/geometry/setHeight 3. cm
/gate/hole/geometry/setRadius .15 cm
/gate/hole/placement/setRotationAxis 0 1 0
/gate/hole/placement/setRotationAngle 90 deg
/gate/hole/setMaterial Air
/gate/hole/vis/forceSolid

#
# Repeat the hole in an array
##
/gate/hole/repeaters/insert cubicArray
/gate/hole/cubicArray/setRepeatNumberX 1
/gate/hole/cubicArray/setRepeatNumberY 52
/gate/hole/cubicArray/setRepeatNumberZ 44
/gate/hole/cubicArray/setRepeatVector 0. 0.36  0.624 cm

# Repeat these holes in a linear
##
/gate/hole/repeaters/insert linear
/gate/hole/linear/setRepeatNumber 2
/gate/hole/linear/setRepeatVector 0. 0.18 0.312 cm

# CRYSTAL
# Create the crystal volume
##
/gate/SPECThead/daughters/name crystal
/gate/SPECThead/daughters/insert box
/gate/crystal/geometry/setXLength 1. cm
/gate/crystal/geometry/setYLength  19. cm
/gate/crystal/geometry/setZLength  28. cm
/gate/crystal/placement/setTranslation  0. 0. 0. cm
/gate/crystal/setMaterial NaI
/gate/crystal/vis/setColor yellow
/gate/crystal/vis/forceSolid

# BACK-COMPARTMENT
# Create the back-compartment volume
##
/gate/SPECThead/daughters/name compartment
/gate/SPECThead/daughters/insert box
/gate/compartment/geometry/setXLength 2.5 cm
/gate/compartment/geometry/setYLength 19. cm
/gate/compartment/geometry/setZLength 28. cm
/gate/compartment/placement/setTranslation   1.75 0. 0. cm
/gate/compartment/setMaterial Glass
/gate/compartment/vis/setColor grey
/gate/compartment/vis/forceSolid

#########################################
# S Y S T E M
######
/gate/systems/SPECThead/crystal/attach crystal
/gate/systems/SPECThead/describe

# Crystal SD
/gate/crystal/attachCrystalSD
/gate/compartment/attachPhantomSD
/gate/shielding/attachPhantomSD
/gate/SPECThead/attachPhantomSD
/gate/collimator/attachPhantomSD

（3）体模及放射源设置

使用两种体素化模型，一种是材质均匀且放射源均匀分布的Shepp-Logan头模型，另一种是真实的患者数据。两种模型分辨率均为128×128×128。Shepp-Logan头模型如图5所示。通过调整模型内部椭球的位置、大小、形状以及放射源活度来获取多组不同的源数据。在仿真中设置源发射光子能量为140 keV，模拟放射性核素${}^{99}T{c}^m$。
真实患者数据在医学图像公开数据集：https://wiki.cancerimagingarchive.net/pages/viewpage.action?pageId=70224216下载，部分图像如图6、7所示。

图5：Shepp-Logan模型

图6：真实患者CT图像

图7：真实放射源分布图像

GATE中提供了三种体素化模型的导入方法，使用使用Regular parameterizatigon method导入体素化模型，使用imageReader方法导入放射源分布数据。

GATE支持多种数据格式的导入，包括：
ASCII、Interfile格式、Analyze格式、DICOM格式、MetaImage格式。在本次仿真中使用Interfile格式的数据导入。

导入体素化模型的代码如下：

# 1-CT16-128.h33为CT图像的header文件，图像以二进制存储-1-CT16-128.i33
# V O X E L I Z E D  M A T R I X   H O F F M A N   B R A I N   P H A N T O M
# Generate materials from Hounsfield units
#/gate/HounsfieldMaterialGenerator/SetMaterialTable                  data/SimpleMaterialsTable.txt
/gate/HounsfieldMaterialGenerator/SetMaterialTable                  data/Schneider2000MaterialsTable.txt
/gate/HounsfieldMaterialGenerator/SetDensityTable                   data/Schneider2000DensitiesTable.txt
/gate/HounsfieldMaterialGenerator/SetDensityTolerance               0.1 g/cm3
/gate/HounsfieldMaterialGenerator/SetOutputMaterialDatabaseFilename data/patient-HUmaterials.db
/gate/HounsfieldMaterialGenerator/SetOutputHUMaterialFilename       data/patient-HU2mat.txt
/gate/HounsfieldMaterialGenerator/Generate

/gate/world/daughters/name VoxPhantom
/gate/world/daughters/insert ImageRegularParametrisedVolume
/gate/VoxPhantom/geometry/setImage data/1-CT16-128.h33
/gate/geometry/setMaterialDatabase data/patient-HUmaterials.db
/gate/VoxPhantom/geometry/setHUToMaterialFile data/patient-HU2mat.txt

/gate/VoxPhantom/placement/setTranslation  0. 0. 0. mm
/gate/VoxPhantom/placement/setRotationAxis 1 0 0
/gate/VoxPhantom/placement/setRotationAngle 0 deg
/gate/VoxPhantom/attachPhantomSD

/gate/VoxPhantom/setSkipEqualMaterials 1

导入体素化放射源的分布代码如下：

#
#1-PET128.h33为放射源分布图像的header文件
# V O X E L   S O U R C E   B A S E D  O N   T H E   H O F F M A N   B R A I N   P H A N T O M

/gate/source/addSource VoxelizedSource voxel

/gate/source/VoxelizedSource/reader/insert image
/gate/source/VoxelizedSource/imageReader/translator/insert linear
/gate/source/VoxelizedSource/imageReader/linearTranslator/setScale 2.6 Bq

/gate/source/VoxelizedSource/imageReader/readFile data/1-PET128.h33
/gate/source/VoxelizedSource/setPosition   -128 -128 -128 mm

/gate/source/VoxelizedSource/setType backtoback
/gate/source/VoxelizedSource/gps/particle gamma
/gate/source/VoxelizedSource/gps/energytype Mono
/gate/source/VoxelizedSource/gps/monoenergy 140. keV
/gate/source/VoxelizedSource/gps/angtype iso
/gate/source/VoxelizedSource/gps/mintheta 0. deg
/gate/source/VoxelizedSource/gps/maxtheta 90. deg
/gate/source/VoxelizedSource/gps/minphi 0. deg
/gate/source/VoxelizedSource/gps/maxphi 360. deg
/gate/source/VoxelizedSource/gps/confine NULL

/gate/source/list

（3）设置物理过程以及电子学参数

在SPECT仿真中定义系统输运参数即物理过程包括：电子对效应、康普顿散射、瑞利散射、电离作用、轫致辐射以及多重散射，并且设置阈值，只需要记录射程在一定范围内的粒子的信息。在GATE中数字化模块又称为数字化链，这一部分模拟了前端电子学电路的功能，可进行一系列的信号处理，可将探测到的粒子信息转换为可观测的物理量，如能量、时间、位置等。

表2：数字化模块参数值

物理过程设置代码如下：

#  P H Y S I C S
#####
/gate/physics/addProcess PhotoElectric
/gate/physics/processes/PhotoElectric/setModel StandardModel

/gate/physics/addProcess Compton
/gate/physics/processes/Compton/setModel PenelopeModel

/gate/physics/addProcess RayleighScattering
/gate/physics/processes/RayleighScattering/setModel PenelopeModel

/gate/physics/addProcess ElectronIonisation
/gate/physics/processes/ElectronIonisation/setModel StandardModel e-

/gate/physics/addProcess Bremsstrahlung
/gate/physics/processes/Bremsstrahlung/setModel StandardModel e-

/gate/physics/addProcess MultipleScattering e-

/gate/physics/processList Enabled
/gate/physics/processList Initialized

#  C U T S
#####
# Cuts for particle in WORLD
##
/gate/physics/Gamma/SetCutInRegion      SPECThead 0.1 cm
/gate/physics/Electron/SetCutInRegion   SPECThead 1.0 cm

#/gate/physics/Gamma/SetCutInRegion      VoxPhantom 0.1 mm
#/gate/physics/Electron/SetCutInRegion   VoxPhantom 0.1 mm

# I N I T I A L I Z A T I O N
#####
# 仿真初始化
/gate/run/initialize

数字化模块设置代码如下：

# D I G I T I Z E R
#####

/gate/digitizer/Singles/insert adder
/gate/digitizer/Singles/insert blurring
/gate/digitizer/Singles/blurring/setResolution 0.10
/gate/digitizer/Singles/blurring/setEnergyOfReference 140. keV
/gate/digitizer/Singles/insert spblurring
/gate/digitizer/Singles/spblurring/setSpresolution 2.0 mm
/gate/digitizer/Singles/spblurring/verbose 0
/gate/digitizer/Singles/insert thresholder
/gate/digitizer/Singles/thresholder/setThreshold 20. keV
/gate/digitizer/Singles/insert upholder
/gate/digitizer/Singles/upholder/setUphold 190. keV

（4）设置仿真输出以及随机数

在仿真过程中可以使用并行计算加快速度，使用并行计算只需要修改随机数种子即setEngineSeed。每次使用不同的setEngineSeed，将模拟得到的计数结果相加即可。

# O U T P U T
#####
#   R A N D O M
# JamesRandom Ranlux64 MersenneTwister
/gate/random/setEngineName Ranlux64
#/gate/random/setEngineSeed default
#/gate/random/setEngineSeed auto
/gate/random/setEngineSeed 123456869
#/gate/random/resetEngineFrom fileName
/gate/random/verbose 1

# P R O J E C T I O N
#####
/gate/output/projection/enable
/gate/output/projection/setFileName out/YourProjection1
/gate/output/projection/pixelSizeX 0.148 cm
/gate/output/projection/pixelSizeY 0.218 cm
/gate/output/projection/pixelNumberX 128
/gate/output/projection/pixelNumberY 128

# Specify the projection plane (XY, YZ or ZX)
##
/gate/output/projection/projectionPlane YZ

# E X P E R I M E N T 
#####
#  (Proj = 64)=180 deg (Rotation speed = 2.8125 deg/sec)
/gate/application/setTimeSlice    4.0 s
/gate/application/setTimeStart     0. s
/gate/application/setTimeStop    64 s
# V E R B O S I T Y
####
/control/execute Verbose.mac
# L E T' S   R U N   T H E   S I M U L A T I O N  !
#####
/gate/application/startDAQ