MCGPULite_v1.5.github.io

MCGPULite_v1.5使⽤⼿册

项目地址：https://github.com/SEU-CT-Recon/MCGPULite_v1.5?tab=readme-ov-file

简介

MCGPU主要⽤途是利⽤GPU多线程计算能⼒，模拟在相关模型中传输X射线，⽣成（⼈体解剖学）真实模型的合成放射图像。

MCGPULite_v1.5是https://github.com/DIDSR/VICTRE_MCGPU/的轻量版本，只允许在Linux系统运行。

其他版本入口:

此使⽤⼿册为作者学习使⽤该程序时整理，对于其中内容的理解还较为粗浅，若有错误的地⽅欢迎批评指正。

代码的一些主要特点

MCGPULite_v1.5 vs MCGPULite

相比于https://github.com/SEU-CT-Recon/MCGPULite，MCGPULite_v1.5新增了一些改进：

MCGPULite_v1.5 vs VICTRE_MCGPU

相比于https://github.com/DIDSR/VICTRE_MCGPU/，MCGPULite_v1.5做了一下删减和修改：

MCGPULite 环境配置

注：以下的环境配置都是基于linux系统。

输⼊⽂件与参数

本程序的输⼊⽂件路径和输⼊参数集中在了.in⽂件中，相关模板为MCGPULite的MCGPULite_v1.3.sample.in或MCGPU-for-MATLAB的MCGPU_MATLAB_sample.in，后者相较于前者多了一小部分内容，本节将会全部提及。

输入文件准备

x-ray能谱⽂件

x 射线源定义为发射 x 射线的点源，其能量从⽤户提供的能谱中随机采样。

关于能谱⽂件的产⽣，可以使⽤程序https://github.com/I-STAR/SPEKTR，在matlab中全部添加路径后运⾏spektr.m⽂件即可。

0

在区域1更改参数可以形成不同的能谱（其他参数可⾃⾏探索），点击按钮Generate Spectrum产⽣能谱，点击按钮Save保存能谱⽂件为txt：

1

根据MCGPULite的输⼊要求，需要做以下更改：

  1. 能量范围需要根据材料的要求进⾏更改，本项⽬中设置为5 ~120keV，且单位改为eV

  2. 去除单位光⼦数（Photons / mm^2 / mAs）为0的⾏，并在最后添加⼀⾏单位光⼦数为负值的项标志能谱⽂件的结束。

经过修改后数据如下（仅展⽰头尾）:

2

3

体素化几何模型文件

体素化几何模型是将每层CT图的各个部位填入相对应的编号和密度，需要按照如下⽰例的格式进⾏创建：

[SECTION VOXELS HEADER v.2008-04-13]
411  190  113      No. OF VOXELS IN X,Y,Z
5.000e-02  5.000e-02  5.000e-02    VOXEL SIZE (cm) ALONG X,Y,Z
1                  COLUMN NUMBER WHERE MATERIAL ID IS LOCATED
2                  COLUMN NUMBER WHERE THE MASS DENSITY IS LOCATED
1                  BLANK LINES AT END OF X,Y-CYCLES (1=YES,0=NO)
[END OF VXH SECTION]
1 0.00120479
1 0.00120479
...

文件头信息中，有五个参数，依次为体素（体积元素）在X、Y、Z轴上的数量，单个体素在X、Y、Z轴上的长度，材料ID所在列，材料质量密度所在列，X、Y循环结束时是否有空行

头信息后为两列数据，每一行分别为材料ID和该材料的密度，材料ID为输入的材料文件在.in中的顺序1,2,3···

注意，ndarray与模体的关系为：[nSlice=VoxZ][height=VoxY][width=VoxX]

下面展示使用python模拟生成water plate，其他示例代码可通过https://github.com/IIIIIsak/vox_vol_interconvert自取。

# Generate water plate phantom .vox

import numpy as np

# The whole volume.

VoxX = 256
VoxY = 256
VoxZ = 256
VoxelSize = 0.125 # cm
WaterMaterialId = 2 # Material index in .in file

for WaterThickness in range(11, 21):  # cm
    # The water plate.
    ModelX = 128  # 4cm
    ModelY = int(WaterThickness / VoxelSize)
    ModelZ = 128
vol = np.ones((VoxX, VoxY, VoxZ), dtype=int)
xl = (VoxX - ModelX) // 2
yl = (VoxY - ModelY) // 2
zl = (VoxZ - ModelZ) // 2

vol[xl:xl + ModelX, yl:yl + ModelY, zl:zl + ModelZ] = WaterMaterialId

voxFile = ""
voxFile += "[SECTION VOXELS phantomER]\n{} {} {} No. OF VOXELS IN X,Y,Z\n{} {} {} VOXEL SIZE (cm) ALONG X,Y,Z\n".format(
    VoxX, VoxY, VoxZ, VoxelSize, VoxelSize, VoxelSize)
voxFile += "1 COLUMN NUMBER WHERE MATERIAL ID IS LOCATED\n"
voxFile += "2 COLUMN NUMBER WHERE THE MASS DENSITY IS LOCATED\n"
voxFile += "0 BLANK LINES AT END OF X,Y-CYCLES (1=YES,0=NO)\n"
voxFile += "[END OF VXH SECTION]\n"

Rhos = [None, '0.00120479', '1'] # Density

vol = vol.flatten()
# Save corresponding .raw volume.
# Use ImageJ Volume Viewer to display. 
# Ray points to +y in ImageJ coordinate when "0.0 1.0 0.0  # SOURCE DIRECTION COSINES: U V W"
with open('WaterPhantom.raw', 'wb') as fp:
    fp.write(vol.astype(np.uint8).tobytes())

# Generate .vox file.
for i in range(vol.shape[0]):
    voxFile += "{} {}\n".format(str(vol[i]), Rhos[vol[i]])

SavePath = '/somewhere/' + str(WaterThickness) + 'cm.vox'
with open(SavePath, 'w') as fp:
    fp.write(voxFile)

print('Done {} cm.'.format(WaterThickness))
材料文件

MC-GPU 使用基于 PENELOPE 数据库的材料属性数据库，已经提供了一组通常用于医学成像模拟的材料的预定义材料文件“material”

目前已有的材料文件：

金属：铝、铯、钨、铅、钢、钛

空气、水、草酸钙、PMMA、硒、聚碳酸酯PC

人体器官：脂肪、血液90_碘10、血液、血_脾、骨头、脑袋、乳房、软骨、眼睛、腺、肾脏、肝、肺、肌肉、红骨髓、皮肤、软组织、骨松质、胃肠

输入参数说明

1. 模拟配置 SIMULATION CONFIG

#[SECTION SIMULATION CONFIG v.2009-05-12]
1.0e10                           # TOTAL NUMBER OF HISTORIES, OR SIMULATION TIME IN SECONDS IF VALUE < 100000
20220216                      # RANDOM SEED (ranecu PRNG)
2                               # GPU NUMBER TO USE WHEN MPI IS NOT USED, OR TO BE AVOIDED IN MPI RUNS
512                             # GPU THREADS PER CUDA BLOCK (multiple of 32)
1000                             # SIMULATED HISTORIES PER GPU THREAD

2. x射线源 SOURCE

#[SECTION SOURCE v.2011-07-12]
./80kVp_bowtie.spec    # X-RAY ENERGY SPECTRUM FILE
8.0   -20.0   8.0              # SOURCE POSITION: X Y Z [cm]
0.0   1.0    0.0                # SOURCE DIRECTION COSINES: U V W
-1 -1                          # TOTAL AZIMUTHAL (WIDTH, X) AND POLAR (HEIGHT, Z) APERTURES OF THE FAN BEAM [degrees] (input negative to automatically cover the whole detector)
0  -1   0             # EULER ANGLES (RxRyRz) TO ROTATE RECTANGULAR BEAM FROM DEFAULT POSITION AT Y=0, NORMAL=(0,-1,0)

3. 图像探测器 IMAGE DETECTOR

./out/somewhat                # OUTPUT IMAGE FILE NAME
512    512                      # NUMBER OF PIXELS IN THE IMAGE: Nx Nz
25     25                     # IMAGE SIZE (width, height): Dx Dz [cm]
50.0                           # SOURCE-TO-DETECTOR DISTANCE (detector set in front of the source, perpendicular to the initial direction)
0.0    0.0                     # IMAGE OFFSET ON DETECTOR PLANE IN WIDTH AND HEIGHT DIRECTIONS (BY DEFAULT BEAM CENTERED AT IMAGE CENTER) [cm]
 0.0200                         # DETECTOR THICKNESS [cm] (DEBUG!!)
 0.004027  # ==> MFP(Se,19.0keV)   # DETECTOR MATERIAL MEAN FREE PATH AT AVERAGE ENERGY [cm] (DEBUG!!)
0.05  3.51795           # PROTECTIVE COVER THICKNESS (detector+grid) [cm], MEAN FREE PATH AT AVERAGE ENERGY [cm] (input 0 or neative to disable)
10.0   78.74   0.0017            # ANTISCATTER GRID RATIO, FREQUENCY, STRIP THICKNESS [X:1, lp/cm, cm] (enter 0 to disable the grid)
0.0157   1.2521   # ==> MFP(lead&polystyrene,19keV)  # ANTISCATTER STRIPS AND INTERSPACE MEAN FREE PATHS AT AVERAGE ENERGY [cm]
1                              # ORIENTATION 1D FOCUSED ANTISCATTER GRID LINES: 0==STRIPS PERPENDICULAR LATERAL DIRECTION (mammo style); 1==STRIPS PARALLEL LATERAL DIRECTION (CBCT style)

4. 扫描轨迹 SCAN TRAJECTORY

3                               # NUMBER OF PROJECTIONS (beam must be perpendicular to Z axis, set to 1 for a single projection)
20.0                            # SOURCE-TO-ROTATION AXIS DISTANCE (rotation radius, axis parallel to Z)
120.0                            # ANGLE BETWEEN PROJECTIONS [degrees] (360/num_projections for full CT)
0                           # ANGULAR ROTATION TO FIRST PROJECTION (USEFUL FOR DBT, INPUT SOURCE DIRECTION CONSIDERED AS 0 DEGREES) [degrees]
 0.0  0.0  1.0                  # AXIS OF ROTATION (Vx,Vy,Vz)
 0.0                            # VERTICAL TRANSLATION BETWEEN PROJECTIONS (HELICAL SCAN)

5. 模型文件 VOXELIZED GEOMETRY FILE

/home/zhuoxu/workspace/scatter-mc/WaterPlate3cm.vox          # VOXEL GEOMETRY FILE (penEasy 2008 format; .gz accepted)
0.0    0.0    0.0              # OFFSET OF THE VOXEL GEOMETRY (DEFAULT ORIGIN AT LOWER BACK CORNER) [cm]

在本程序中，坐标系原点并不是物体的中心，而是物体“包装盒”的后下角，即整体只处于第一象限内，如图所示:

5

在设置其他参数时，需要以此坐标系特征为前提。

假设有一个大小为16*16*16cm的圆柱形水模，源初始方向为(0,1,0)，我们在已知SOURCE-TO-ROTATION AXIS DISTANCE（SOD）的情况下，如何设置射线源的初始位置才能使得旋转轴穿过模型的中心呢？

图片1

如图，射线源的位置为($\frac x2,\frac y2 - SOD,\frac z2$)

间隔120°得到的三张投影结果如下：

7

6. 材料文件列表 MATERIAL FILE LIST

/home/kanshengqi/software/MCGPULite/material/air__5-120keV.mcgpu.gz                         #  1st MATERIAL FILE (.gz accepted)
/home/kanshengqi/software/MCGPULite/material/water__5-120keV.mcgpu.gz                             #  2nd MATERIAL FILE

运行

# .in file. Specify which CUDA GPU to run on here.

3              # GPU NUMBER TO USE WHEN MPI IS NOT USED, OR TO BE AVOIDED IN MPI RUNS

# Command line. Just as it should be.

MCGPULite1.5 ./MCGPULite_v1.5.sample.in

# .in file. Use -1 always here.

-1              # GPU NUMBER TO USE WHEN MPI IS NOT USED, OR TO BE AVOIDED IN MPI RUNS

# Command line. Use CUDA_VISIBLE_DEVICES to specify GPUs, and pass GPU counts to -np.

mpirun -x CUDA_VISIBLE_DEVICES=0,1,2,3 -np 4 MCGPULite1.5 ./MCGPULite_v1.5.sample.in

经过实验，如果 TOTAL NUMBER OF HISTOIES 小于 1e10，因为分布开销，不建议使用多个 GPU。如果光子数达到 1e11 ，时间会缩短了66%。

如果出现类似error while loading shared libraries: libcudart.so.10.0: cannot open shared object file，可以参考此篇博客(23条消息) error while loading shared libraries libcudart.so.8.0: No such file or directory问题解决_坎幽黑尔弥？的博客-CSDN博客

参考修改：

vim .bashrc
在这个文件中添加上述依赖项的位置目录：
export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/usr/local/cuda-10.0/
export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/usr/local/cuda-10.0/lib64/
export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/home/wuxian/app/openmpi-4.1.1/install/lib/
保存并退出之后执行：
source .bashrc

输出文件

MCGPULite_v1.5输出的每个raw文件里包含三张图，依次为：Total(Primary+ Scatter), Primary, Scatter，即含散射的投影，不含散射的投影，散射图， 单个文件导入imageJ中时注意设置Number of images为3。 MCGPU-for-MATLAB 可直接通过代码自行选择保存。

8

除此外，MCGPULite_v1.5还会输出每张投影的模拟报告的二进制文件，只能通过imageJ打开，该报告包含了该张投影图的几何参数、光子到达探测器的平均概率、模拟的光子数、模拟时间等参数。

image-20241018144708511

已知的问题

亮场模拟bug

在模拟亮场时，一般是将模体做成一个全是空气密度的模体，然后其他参数不变就可以了；但是在MCGPU_v1.5b版本中，如果射线源有光子穿过了空体，似乎会发生一些“奇怪”的反应，导致生成的结果出现环形的噪声，目前暂时不知道应该怎么修复，但是可以通过设置射线源的位置，比如“拉高”射线源位置使得光子可以不穿过模体直接打到探测器上，这样模拟出来的亮场是没有问题的。