STIR散射校正使用说明

STIR(Software for Tomographic Image Reconstruction)是一款受业界认可的开源医学图像重建软件，使用C++编写，在linux平台下可实现PET图像重建、数据校正等功能。（源代码：https://github.com/tokkot/STIR)

为验证自己的散射校正算法正确性，就想使用STIR来测试一下，但发现STIR使用的正弦图存储格式，图像存储方式与我现在用的不同，而且修改脚本后总是报错且结果总是不正确。为此我阅读了STIR散射校正部分的代码，并在此记录一下STIR散射校正算法的使用方法。

1、  STIR安装

(1)    在STIR官网上(http://stir.sourceforge.net/)下载.zip压缩包，解压后发送到linux虚拟机中

(2)    /STIR/src/CMakeLists.txt中设置"Compile with OpenMP" 为ON，使能openmp功能

(3)    STIR使用boost库，所以在编译之前要先安装boost。安装指令：”sudo apt-get install libboost-dev”

(4)    在STIR/下新建build文件夹，在build/中依次执行”sudo cmake ..”, “sudo make -j5”, “sudo make install”. 其中j5表示使用5个线程进行make，可以大大加快速度。

(5)    安装完成后会在/usr/local/bin中生成大量可执行文件，若”echo $PATH”后发现该目录在PATH下，则可以直接调用STIR中命令，即安装完成。STIR使用示例：”simulate_scatter xxx.par”、”estimate_scatter xxx.par”

2、  STIR数据输入

(1)    投影数据

a)     有关投影数据的专业术语(详见STIR-glossary)

• Michelogram：即我们常说的4D正弦图数据

• Sinogram: 即我们常说的2D正弦图数据，在STIR中slice表示图像横截面

• Segment: 具有相同环差的一组sinogram为一个segment，在下图中倾斜的黑色直线表示环差相同的一组sinogram

• Span: 轴向压缩(Axial compression)，将具有相同平均环差((ring1+ring2)/2)的sinogram相加，从而减少sinogram的数目，一般对奇数个segment进行span，如下图，7个segment融合成一个segment。融合后的每个segment包含的sinogram是不同的，具体计算公式现在不清楚。若span=1，即不做轴向压缩，michelogram中包含环差从-ringNum+1 : ringNum -1共2*ringNum-1个segment。

• 下图中rdmin和rdmax分别为最小和最大环差，表示哪些环差范围的sinogram进行了融合；average_delta表示平均环差，即((ring1+ring2)/2)。

b)    STIR的投影数据格式

• STIR使用修改过的INTERFILE数据格式(STIR专用-_-)读入投影数据，包含2个文件，头文件(.hs)和数据文件(.s)

• 可以使用”create_projdata_template”命令生成头文件

• 投影数据头文件中包含3部分数据：实验数据(设备名称、患者姿态、扫描时间、能窗信息)，投影数据信息(数据格式、segment, view, bin的数目)和探测器结构(环数，每环晶体数，能量分辨率等)

• STIR以segment为单位管理投影数据。读取数据时，segment默认存储顺序由头文件中的”minimum ring difference per segment”和”maximum ring difference per segment”决定

• 一个segment中包含多个sinogram，sinogram个数由头文件的”axial coordinate”决定，STIR默认sinogram根据所在环的位置从小到大排列，所以存储数据时要注意顺序

• Sinogram中的view数目一般为CrystalNumPerRing/2；bin的数目最大为CrystalNumPerRing-1，但一般根据FOV大小在两边减去部分(最好减去，可以减少计算量；不减的话可能会报错)

• STIR支持的segment的实际物理存储顺序有两种：Segment_View_AxialPos_TangPos和Segment_AxialPos_View_TangPos。例如，一个segment中包含10个sinogram，1个sinogram包含20view和39bin。前一种方式先存储这10个sinogram中view=0的所有bin，再存储view1…view10；后一种方式先存储第1个sinogram，再存储第2个sinogram，依次类推。在InterfilePDFSHeader::find_storage_order()中说明了，投影数据的存储顺序是由interfile的头文件确定的。matrix axis label [3] := view; matrix axis label [2] := axial coordinate对应Segment_View_AxialPos_TangPos。matrix axis label [3] := axial coordinate; matrix axis label [2] := view对应Segment_AxialPos_View_TangPos。

• 在Scanner.h中有crystal, block等探测器模块的介绍。crystal: 最小探测单元；block: 多个晶体组成block；layer: DOI探测器中有多个layer；bucket: 多个block探测到数据会由同一个bucket发送；singles_unit: 绝大多数探测器可探测到单光子，有的探测器每个晶体可以输出单光子，有的是每bucket输出，singles_units就是可以输出单光子计数率的晶体模块

• 注意：头文件以”;”作为注释符号，必须在单独行中进行注释，不能在实际参数后进行注释，否则读取参数错误

c)    投影数据头文件部分参数说明

• 示例头文件见附件，下面只介绍部分需要设置的参数

• name of data file := AttnCoff_STIR.s     ;数据文件名

• !number format := float         ;数据类型

• !number of bytes per pixel := 4     ;每个数据所占字节数

• number of dimensions := 4   ;数据维数，一般是4维，segment, segment, view, axial coordinate

• matrix axis label [4] := segment  ;第4维的名称，不能改变，顺序检索投影数据时，该维度变化最慢

• !matrix size [4] := 95     ;第4维的数目，若不进行span，对48环的脑PET，三维投影数据包含2*48-1=95个环差，也即有95个segment

• matrix axis label [3] := view ;第3D为view，说明存储顺序是Segment_View_AxialPos_TangPos

• !matrix size [3] := 132   ;共有132个view

• matrix axis label [2] := axial coordinate      ;第2D为轴向坐标

• !matrix size [2] := {1, 2, 3, … , 45, 46, 47, 48, 47,46, 45, … ,3, 2, 1}     ;每个segment包含的sinogram数目，如有95个segment，则该数组的长度为95

• matrix axis label [1] := tangential coordinate     ;第1维为bin，顺序检索投影数据时，该维度变化最快

• !matrix size [1] := 191   ;每个view包含191个bin

• minimum ring difference per segment := { -47,-46,-45, … , 45,46,47 }

• maximum ring difference per segment := { -47,-46,-45, … , 45,46,47 }    ;与上式共同表示每个segment中包含的最大和最小环差，若没有span，则这两项是相同的。如有95个segment，则该数组的长度为95

• number of energy windows:=1

• energy window lower level[1]:=410

• energy window upper level[1]:=700  ;设置能窗

• Energy resolution := 0.19

• Reference energy (in keV) := 511         ;设置能量分辨率

• Number of rings := 48   ;探测环数目

• Number of detectors per ring := 264   ;每环晶体数

• Inner ring diameter (cm) := 37.7          ;探测环内径

• Average depth of interaction (cm) := 1.0     ;平均探测深度

• Distance between rings (cm) := 0.42   ;每环的间隔

• Maximum number of non-arc-corrected bins := 191          ;随便设？

• Default number of arc-corrected bins := 191      ;随便设？

• Number of blocks per bucket in transaxial direction := 1 ;下面这些参数在散射校正中没有使用

• Number of blocks per bucket in axial direction := 2

• Number of crystals per block in axial direction := 6

• Number of crystals per block in transaxial direction := 6

• Number of detector layers := 1

• Number of crystals per singles unit in axial direction := 1

• Number of crystals per singles unit in transaxial direction := 1

(2)    图像数据

a)    STIR的图像数据格式

• STIR同样使用修改过的INTERFILE数据格式(STIR专用-_-)读入图像数据，包含2个文件，头文件(.hv)和数据文件(.v)

• STIR只考虑环形探测器，坐标系与探测器相关，以探测器轴向为z轴，纵向y轴，横向x轴。探测器坐标轴原点位于探测器中心。图像坐标原点默认为第一个平面中心体素的中心，以远离床的那一边的xy平面作为第一个平面

• STIR中默认图像X和Y方向的体素数目为奇数；Z方向的体素数目为2*ringNum-1，Z方向体素尺寸为晶体尺寸的一半。若XY方向体素数目为偶数，图像原点位于XY轴正方向第一个体素的中心，不建议设置为偶数。Z方向体素数目最好设置为默认值，这样可以使用代码中的对称，且不会报错。

• 对图像进行下采样时(例如将尺寸为(127,127,95)的衰减图像下采样为(17,17,13)的图像以用于获取散射点)，需要保证边缘体素中心到图像中心距离一致，而不是图像边缘到图像中心一致，即：

• (image_size1-1)*voxel_size1/2 = (image_size2-1)*voxel_size2/2

• 可以使用AMIDE直接读入.hv文件

• 注意：使用SimpleITK重采样图像时，可能会有Y轴镜像翻转的情况，目前还不清楚具体原因，转换后应检查图像！

b)    图像数据头文件部分参数说明

• 示例头文件见附件，下面只介绍部分需要设置的参数

• name of data file := Zubal_Umap_17_17_13.v     ;图像数据文件名

• !number format := float         ;数据存储格式

• !number of bytes per pixel := 4     ;每个数据所占字节数

• number of dimensions := 3   ;数据维度

• matrix axis label [1] := x        ;第1维为x

• !matrix size [1] := 17     ;x轴有17个体素

• scaling factor (mm/pixel) [1] := 16.5375    ;x轴体素尺寸

• matrix axis label [2] := y        ;第2D为y

• !matrix size [2] := 17     ;y轴有17个体素

• scaling factor (mm/pixel) [2] := 16.5375    ;y轴体素尺寸

• matrix axis label [3] := z        ;第3D为z轴

• !matrix size [3] := 13     ; z轴有13个体素

• scaling factor (mm/pixel) [3] := 16.45         ; z轴体素尺寸

• first pixel offset (mm) [1] := -132.3     ;索引值为0的体素中心距图像中心的距离

• first pixel offset (mm) [2] := -132.3

• first pixel offset (mm) [3] := 0

• number of energy windows := 1

• energy window lower level[1] := 410

• energy window upper level[1] :=  700      ;设置能窗

3、  STIR散射校正

(1)    STIR散射校正整体说明及注意事项

• STIR散射校正分为2部分，simulate scatter和estimate scatter

• simulate scatter是实际使用SSS估计散射事件的代码，输入活度图像、衰减图像和模板投影数据头文件，输出散射估计Michelogram，输出格式由模板投影数据头文件决定。

• simulate scatter使用2种方法减少计算量。首先是下采样探测器和使用放大的衰减图像设置散射点

• 可以自定义下采样探测器结构，默认每环64个晶体，轴向20mm设置一环晶体，从而减少LOR数目。该方法需与estimate scatter中的upsample_and_fit_single_scatter结合使用，由于上采样函数只能对segment = 0的投影数据（直层投影数据）上采样，所以为了下采样探测器，3D投影数据要先SSRB为2D数据后才能进行散射估计

• STIR是在衰减图像的每个体素中采集一个散射点，所以使用放大的衰减图像可以减少散射点。若不设置放大的衰减图像，STIR根据下采样探测器（若设置了）的结构，放大衰减图像（XY方向的体素尺寸约为bin的间距，Z方向体素尺寸为晶体尺寸的一半）；还可以通过设置放大参数（参数文件中的4个zoom参数）放大衰减图像（建议不使用）；还可以直接输入放大后的衰减图像文件名

• 只使用simulate scatter是可以进行3D散射估计的，但不能下采样探测器，而且由于使用3D投影数据，计算时间大大增加；另外simulate scatter不包含尾部拟合

• estimate scatter在simulate scatter的基础上规定了一套估计散射事件的流程，包含SSRB，上采样等操作

• 对输入的3D投影数据进行SSRB，转化为2D投影数据

• 设置图像重建输入，设置2D投影数据作为输入参数，将衰减校正因子（>1），归一化校正因子，随机校正因子SSRB后加入到图像重建中

• 下采样探测器，如果设置了

• 设置尾部拟合掩膜。依据衰减图像计算mask_image（通过阈值设置0-1图像）；对mask_image进行滤波后，正投影得到mask_projdata，并加1（create_tail_mask_from_ACFs是处理衰减校正因子的，加1模拟衰减校正因子）；对mask_projdata中每个view两端小于阈值（默认1.1）的bin置1，中间部分置0，用于尾部拟合

• 计算初始活度图像（无法在参数文件中设置），使用SSRB后的2D投影数据作为输入，全1图像作为初始值

• 对2D投影数据进行simulate scatter

• 进行上采样，将下采样探测器后得到的投影数据上采样为SSRB后的2D投影数据格式，使用二次B样条采样方法，将投影数据看作图像，确定新投影数据中的bin在原数据中的位置，然后进行插值

• 对2D投影数据进行尾部拟合，每个直层采用实际数据与估计数据直接相除的方式计算一个缩放因子，再使用一个算术平均滤波器进行平滑，缩放因子的范围由minimum /maximum scatter scaling factor决定

• 将2D投影数据inverseSSRB为3D投影数据，直接将平均轴向位置相同的直层数据赋值给斜层

• 3D投影数据除以归一化因子，即散射校正考虑探测器效率因素

• estimate scatter会覆盖simulate scatter参数文件中定义的衰减图像和模板投影数据，其中模板投影数据由SSRB后的2D投影数据决定；而且simulate scatter参数文件中放大衰减图像也无法使用，只能通过下采样探测器减少LOR和放大衰减图像尺寸

• 使用estimate scatter时，simulate scatter参数文件中的” downsample scanner”不能同时置1，不然程序会报错

• 只能通过参数文件来定义重建方式和散射估计方式，而不能像用户指南里说的只输入一个重建类型名或散射估计类型名

• 参数文件中不输入归一化校正因子，只输入衰减校正因子的写法：

  Bin Normalisation type := from projdata

  Bin Normalisation From ProjData :=

  normalisation projdata filename:= AttnCoff_STIR.hs

  End Bin Normalisation From ProjData:=

• 若没有归一化因子，可以自己创建一个全1归一化校正因子，需要在边缘处设置一个不为1的数。因为代码会判断全1归一化因子，会报错

• 若第一次散射迭代的散射估计很大，那做散射校正时，由于prmp – scat < 0时，会直接置0，相当于没做散射，每次迭代的变化很小。这时，尾部拟合最大缩放限制可以设小一点，不让散射估计变得太大，从而让散射估计的变化加快。现在尝试的是第3次迭代效果最好

(2)    simulate scatter参数文件部分参数说明

• 示例头文件见附件，下面只介绍部分需要设置的参数

• attenuation_threshold :=.01         ;选取散射点的阈值，单位为cm-1，输入衰减图像单位也是cm-1

• randomly_place_scatter_points := 1     ;在衰减图像内随机取点

• use_cache := 1         ;缓存散射点到晶体的积分

• activity_image_filename := BrainZubal_Emap_STIR.hv      ;活度图像文件名

• attenuation_image_filename := BrainZubal_Umap_STIR.hv      ;衰减图像文件名

• template_projdata_filename :=  ZubalPrmp_STIR.hs        ;模板投影数据头文件

• output_filename_prefix := Scat_Coarse_STIR      ;输出文件名

• downsample scanner := 0     ;下采样探测器

• downsampled scanner number of detectors per ring := 64       ;下采样探测器每环晶体数

• downsampled scanner number of rings := 11    ;下采样探测器环数

• attenuation_image_for_scatter_points_filename := Zubal_Umap_17_17_13.hv        ;计算散射点的衰减图像

• zoom XY for attenuation image for scatter points := -1      ;原图像体素尺寸/新图像体素尺寸，应该小于1

• zoom Z for attenuation image for scatter points := -1        ;

• XY size of downsampled image for scatter points := -1      ;新图像体素数目

• Z size of downsampled image for scatter points := -1        ;

• attenuation image for scatter points output filename := zoomed_atten_image       ;放大后图像输出文件名

(3)    estimate scatter参数文件部分参数说明

• 示例头文件见附件，下面只介绍部分需要设置的参数

• run in debug mode := 1 ;输出中间变量到extras文件夹中，输出初始活度图像，未上采样的2D散射估计结果，上采样后的2D散射估计结果，重建结果和生成的中间投影数据

• input file := ZubalPrmp_STIR.hs  ;输入3D投影数据

• attenuation image filename := BrainZubal_Umap_STIR.hv        ;输入衰减图像

• recompute mask image := 1 ;重计算掩膜图像

• mask image filename := mask_image.hv      ;掩膜图像文件名

• mask attenuation image filter filename := postfilter_Gaussian_for_mask.par ;掩膜图像正投影前的滤波器

• mask attenuation image min threshold := 0.003        ;计算掩膜图像时，大于阈值置1，否则置0

• recompute mask projdata := 1      ;重新正投影掩膜图像

• mask projdata filename := mask_projdata.hs     ;掩膜正投影文件名

• tail fitting parameter filename := tail_fitting.par        ;尾部拟合参数文件名

• background projdata filename := ${randoms3d}        ;随机校正因子

• Bin Normalisation type := Chained      ;第一个文件默认归一化因子，第二个默认衰减校正因子

• Chained Bin Normalisation Parameters:=  ;只输入一个文件则默认是衰减校正因子

•     Bin Normalisation to apply first:= from ProjData

•         Bin Normalisation From ProjData :=

•         normalisation projdata filename:= norm.hs

•     End Bin Normalisation From ProjData:=

• Bin Normalisation to apply second:= From ProjData

•     Bin Normalisation From ProjData :=

•         normalisation projdata filename:= AttnCoff_STIR.hs

•     End Bin Normalisation From ProjData:=

• END Chained Bin Normalisation Parameters:=

  
  

• reconstruction parameter filename := run_reconstruction.par        ;重建参数文件

• number of scatter iterations := 2 ;散射估计次数，每次改变的是活度图像

• scatter simulation parameter filename := scatter_simulation.par   ;散射估计参数文件

• use scanner downsampling in scatter simulation := 0       ;下采样探测器

• export scatter estimates of each iteration := 1    ;每次迭代输出结果

• output scatter estimate name prefix := scatter_temp         ;输出文件名

• output additive estimate name prefix:= total_additive_temp   ;散射加随机

• do average at 2 := 1        ;第2次迭代的活度图像是初始活度图像和第一次散射校正后重建图像的平均

• maximum scatter scaling factor := 2    ;尾部拟合时，每个直层的缩放因子范围

• minimum scatter scaling factor := 0.4 ;

• upsampling half filter width := 3 ;尾部拟合缩放因子平均滤波器长度，这里是2*3+1=7，每个参数为1/7

• remove interleaving before upsampling := 0      ;上采样插值时，解耦view中bin的Z字型排列，view翻倍

• run in 2d projdata := 1  ;固定，现在只能对2D投影数据进行散射估计

4、  附件 使用STIR进行散射估计的流程

(1)    实验室投影数据转化为STIR类型

a)    数据文件处理如下，更改类型为”.s”

b)    设置投影文件头文件”.hs”

(2)    设置图像，让XY方向有奇数个体素，Z方向有2*ringNum-1个体素

a)    使用simpleITK重采样图像，更改类型为”.v”

b)    设置图像头文件”.v”

(3)    设置simulate scatter参数文件scatter_simulation.par

(4)    设置estimate scatter参数文件scatter_estimate.par

(5)    设置用于mask_image的滤波器参数文件postfilter_Gaussian_for_mask.par

(6)    设置尾部拟合参数文件tail_fitting.par

(7)    设置图像重建参数文件run_reconstruction.par

(8)    执行散射估计指令