基于装置成像坐标的MR图像配准方法、分析方法及其装置与流程

基于装置成像坐标的mr图像配准方法、分析方法及其装置
技术领域
1.本发明属于磁共振成像领域，具体涉及一种mr图像的配准、分析方法及装置。


背景技术：

2.磁共振化学交换饱和转移(chemical exchange saturation transfer,cest)成像是在磁化转移技术及化学交换理论基础上发展起来的一种新的磁共振成像方法，可产生与常规磁共振成像序列不一样的对比度。其成像原理是利用特定的偏共振饱和脉冲，对特定物质中的可交换质子或分子进行充分的预饱和，这种饱和通过化学交换，进一步影响自由水的信号强度，因此通过检测水的信号，可间接反映这种物质的信息。cest mri可以检测内源性和外源性物质，例如蛋白质和多肽中的酰胺质子(
‑
nh)，谷氨酸和肌酸中的胺质子(
‑
nh2)，糖胺聚糖、糖原和葡萄糖中的羟基质子(
‑
oh)等，具有广泛的研究和应用价值。酰胺质子转移(amideproton transfer,apt)成像是cest mri的一种细分技术，通过检测游离蛋白中的酰胺质子和水的交换，间接反映蛋白质水平。有研究表明，肿瘤恶性程度通常与细胞增殖和密度呈正相关，导致整体蛋白质水平升高，使得这种酰胺质子在肿瘤中的浓度较正常组织中的更高。apt mri已应用于检测人类脑肿瘤，预测胶质瘤级别、评估缺血性中风和评估肿瘤对化疗的反应等。糖胺聚糖化学交换饱和转移技术、谷氨酸化学交换饱和转移技术和葡萄糖化学交换饱和转移技术也是cest领域的重要技术，已分别应用于诊断骨关节炎、癫痫症和脑肿瘤。
3.cest定量分析通常需进行一系列后处理，步骤主要为：(1)对cest图像进行配准以消除运动伪影。由于cest mri通常需要采集不同饱和偏移频率下的多个帧，在量化cest参数时需要采集的帧数尤其多。因此该技术的采集时间很长，容易受呼吸等运动影响；(2)主磁场b0的不均匀部分会导致空间编码错误,扭曲cest成像,产生伪影并严重影响定量分析的准确性,因此需要对cest图像进行b0场不均匀性校正；(3)计算cest量化图；(4)在cest图像上选取感兴趣区域(region of interest,roi)进行量化分析，感兴趣区域一般为脑区或病灶区域。但cest成像的空间分辨率较低，组织对比度较差，难以精准显示解剖结构，尚无法实现cest图像的脑区自动分割，而人工选取roi的要求较高，准确度也不能保障。
4.cest作为磁共振新兴成像方法的重要分支，在分子影像和组织代谢功能测定方面具有巨大的研究和应用潜力。然而在将其应用于临床诊断的过程中，仍然有许多问题亟待解决，比如cest后处理流程较为繁琐，感兴趣区域的准确提取存在困难。这些问题会对后续定量分析的准确度产生影响，不利于临床医生和医学影像研究工作者对cest mri进行研究。


技术实现要素：

5.本发明的目的在于解决现有技术中cest后处理流程较为繁琐，感兴趣区域的准确提取存在困难的缺陷，并提供一种基于装置成像坐标的mr图像配准方法及装置。本发明可将其他高分辨率的mr结构像与cest图像进行配准，借助解剖结构图像的解剖信息提取感兴
趣区域，以便于后续在配准后的cest图像上进行定量分析，提高定量分析的准确性。
6.本发明所采用的具体技术方案如下：
7.第一方面，本发明提供了一种基于装置成像坐标的mr图像配准方法，用于实现mr结构像与cest图像之间的跨模态配准，其中mr结构像和cest图像两种mr图像中一个为参考图像，另一个为待配准图像，该mr图像配准方法包括：
8.s1、基于分别从参考图像和待配准图像的dicom信息中提取的旋转变换矩阵、平移变换矩阵和缩放变换矩阵，计算将参考图像由图像坐标系变换到解剖学空间坐标系所需的第一仿射变换矩阵a1，以及将待配准图像由图像坐标系变换到解剖学空间坐标系所需的第二仿射变换矩阵a2；
9.s2、根据第一仿射变换矩阵a1，将参考图像中每个体素在图像坐标系下的坐标[i1,j1,k1]转换为解剖学空间坐标系下的坐标[x,y,z]：
[0010][0011]
s3、基于待配准图像和参考图像的解剖学空间坐标系的统一性，利用第二仿射变换矩阵a2的逆变换，将参考图像中每个体素在解剖学空间坐标系下的坐标[x,y,z]转换为该体素在待配准图像的图像坐标系下的坐标[i2,j2,k2]：
[0012][0013]
s4、针对参考图像中每个体素经过转换后的坐标[i2,j2,k2]，通过插值方法利用待配准图像中该坐标周边的体素灰度值计算该坐标处的灰度值，从而得到与参考图像中每个体素一一对应的待配准图像，实现待配准图像在参考图像上的跨模态配准。
[0014]
作为第一方面的优选，所述s1中，待配准图像与参考图像由同一患者在同一次mr扫描时采集得到，且两次序列采集之间无明显头动。
[0015]
作为第一方面的优选，所述s1中，参考图像和待配准图像的仿射变换矩阵计算方法相同，步骤如下：
[0016]
s11、从目标图像的dicom信息中提取六元double数组[r
ix
,r
iy
,r
iz
,r
jx
,r
jy
,r
jz
]形式的(0020,0037)image orientation(patient)字段，再由[r
ix
,r
iy
,r
iz
]和[r
jx
,r
jy
,r
jz
]的叉积计算得到[r
kx
,r
ky
,r
kz
]，从而得到旋转变换矩阵rotation：
[0017][0018]
s12、从目标图像的dicom信息中提取三元double数组[t
x
,t
y
,t
z
]形式的(0020,0032)image position(patient)字段，从而得到平移变换矩阵translation：
[0019]
[0020]
s13、从目标图像的dicom信息中提取二元double数组[s
x
,s
y
]形式的(0028,0030)pixelspacing字段，同时提取(0018,0050)slicethickness或(0018,0088)spacingbetweenslices字段作为一元double变量s
z
，从而得到缩放变换矩阵scaling：
[0021][0022]
s14、根据目标图像对应的旋转变换矩阵、平移变换矩阵和缩放变换矩阵，得到该目标图像的仿射变换矩阵a为：
[0023][0024]
其中：若目标图像为参考图像，则a＝a1，若目标图像为待配准图像，则a＝a2。
[0025]
作为第一方面的优选，所述s4中，插值方法包括线性插值或三次样条插值法。
[0026]
作为第一方面的优选，若所述参考图像的层厚slicethickness大于1mm，则采用多轮配准并取均值的方法得到最终配准结果，步骤如下：
[0027]
首先，针对从参考图像的dicom信息中提取的原始平移变换矩阵将其沿图像法向量方向移动n次，n为以mm为单位的参考图像层厚的四舍五入值；n次移动后得到n个新的平移变换矩阵为：
[0028][0029]
然后，以所述新的平移变换矩阵作为从参考图像的dicom信息中提取的平移变换矩阵，重复s1～s4得到一张配准至参考图像上的待配准图像，将得到的n张待配准图像进行平均，作为最终配准参考图像上的待配准图像。
[0030]
作为第一方面的优选，所述s4中，需利用flirt或spm工具将配准得到的待配准图像与参考图像再次进行配准，以校正参考图像和待配准图像采集期间的轻微头动。
[0031]
第二方面，本发明提供了一种基于图谱的cest图像定量分析方法，其具体做法为：首先对采集的cest数据进行配准和校正处理得到cest图像，然后再按照第一方面任一方案所述的mr图像配准方法，将待配准图像和cest图像进行配准；再基于配准后的mr结构像确定感兴趣区域并将感兴趣区域映射至cest图像上，从而计算该感兴趣区域内的cest量化参数。
[0032]
作为第二方面的优选，所述采集的cest数据需首先将cest图像的所有帧配准到目标代谢产物集团所对应的帧上以消除运动伪影，然后进行b0场不均匀性校正，得到cest图像。
[0033]
第三方面，本发明提供了一种基于装置成像坐标的mr图像配准装置，其包括存储器和处理器；
[0034]
所述存储器，用于存储计算机程序；
[0035]
所述处理器，用于当执行所述计算机程序时，实现如上述第一方面任一方案所述的mr图像配准方法。
[0036]
第四方面，本发明提供了一种基于图谱的cest图像定量分析装置，其包括存储器和处理器；
[0037]
所述存储器，用于存储计算机程序；
[0038]
所述处理器，用于当执行所述计算机程序时，实现如上述第二方面任一方案所述的cest图像定量分析方法。
[0039]
本发明相对于现有技术而言，具有以下有益效果：
[0040]
mr结构像以较高的分辨率提供解剖形态信息，但无法反映组织的代谢功能情况；cest图像的分辨率较差，但它提供的组织代谢信息是解剖图像不能替代的。本发明结合了两者的优点，提出一种基于装置成像坐标系的跨模态mr图像配准的方法，将cest图像和mr结构像配准，在高分辨率的mr结构像上提取感兴趣区域，提高感兴趣区域勾画的精确度。然后计算cest图像上对应区域的cest量化指标，本发明可计算多个常用的cest参数：cestr、cestr
nr
、mtr
rex
以及arex，为后续研究提供多种选择。基于本发明可进一步将cest数据后处理、图像配准以及定量分析的各个步骤整合，实现自动处理，为临床医生和医学影像研究工作者研究cest mri提供便利，推进cest mri技术的临床应用的发展步伐。
[0041]
另外，现有的cest图像定量分析大多数针对2d cest图像，本发明既可用于2d cest图像，也可用于3d cest图像。3d cest图像较2d cest图像在定量分析优势明显，对全脑的代谢分析具有极其重要的意义。因此，本发明提出的基于装置成像坐标的跨模态mr图像配准方法，即能自动快速实现3d
‑
3d配准，也能实现2d
‑
3d配准，有效解决了不同序列图像fov、分辨率不一致导致难以相互匹配、融合的问题。
附图说明
[0042]
图1为本发明提供的一种基于图谱的cest图像定量分析的流程图；
[0043]
图2为健康志愿者a的t1wi
‑
mprage图(a)，t1wi
‑
mprage图自动分割后的label图(b)，头动校正后以t1wi
‑
mprage图为参考图像进行配准得到的cest不饱和帧源图像(c)，以t1wi
‑
mprage图为参考图像进行配准得到的t1mapping图(d)，以及用配准后的cest图像和t1 mapping图计算的aptw图(e)，cestr
nr
图(f)，mtr
rex
图(g)和arex图(h)。列为不同层面。
[0044]
图3为脑肿瘤患者的t1wi
‑
mprage图(a)，头动校正后以t1wi
‑
mprage图为参考图像进行配准得到的cest不饱和帧源图像(b)，用配准后的cest图像计算的aptw图(c)，cestr
nr
图(d)，mtr
rex
图(e)。列为不同层面。
[0045]
图4为健康志愿者b头动校正后的2d cest不饱和帧源图像(a)，以2d cest图像为参考图像进行配准后的flair图(b)，t1wi
‑
mprage图(c)和t1wi
‑
mprage
‑
label图(d)，(e)为aptw图。
具体实施方式
[0046]
下面结合附图和具体实施方式对本发明做进一步阐述和说明。本发明中各个实施方式的技术特征在没有相互冲突的前提下，均可进行相应组合。
[0047]
为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图对本发明的具体实施方式做详细的说明。在下面的描述中阐述了很多具体细节以便于充分理解本发明。但是本发明能够以很多不同于在此描述的其它方式来实施，本领域技术人员可以在不违背本发明内涵的情况下做类似改进，因此本发明不受下面公开的具体实施例的限制。本发明各个实施例中的技术特征在没有相互冲突的前提下，均可进行相应组合。
[0048]
在本发明的描述中，需要理解的是，术语“第一”、“第二”仅用于区分描述目的，而不能理解为指示或暗示相对重要性或者隐含指明所指示的技术特征的数量。
[0049]
在通过核磁共振设备对患者做mr检查时，除了cest图像通常还会扫描mr结构像作为参考，理想情况下患者扫描两个序列之间无头部运动，则cest图像和mr结构像基于同一装置成像坐标系。这样，可以将由同一成像设备获取的mr结构像与cest图像这两幅mr图像的配准称为跨模态配准，本发明期望利用两种序列图像的坐标信息进行图像配准：一种方式是将mr结构像配准到cest图像上，即参考图像为cest图像，待配准图像为mr结构像；另一种方式是将cest图像配准到mr结构像上，即参考图像为mri结构像，待配准图像为cest图像。例如，后续实施例中的图2和3是将cest图像配准到mr结构像上，而图4是将mr结构像配准到cest图像上。
[0050]
装置成像坐标系即解剖学空间坐标体系，用来描述标准的人体在解剖学上的位置。这个坐标体系由三个解剖轴组成，分别为人体的前
‑
后(anterior
‑
posterior)，头
‑
脚(superior
‑
inferior)，左
‑
右(left
‑
right)。本发明以lps为坐标轴，即从右到左(x轴)，从前到后(y轴)，从脚到头(z轴)，以下称解剖学坐标系为lps坐标系。图像坐标体系描述一幅与解剖相关的图像是怎样被获取的。医学扫描仪器创建了规则的点和网格的矩形数组，它的原点在左上角，图像的三条坐标轴中：i坐标轴向右递增，j坐标轴向下递增，k坐标轴向后递增，以下称图像坐标系为ijk坐标系。
[0051]
医学数字成像和通信标准(digital imaging and communication in medicine,dicom)，是医学图像和相关信息的国际标准(iso 12052)，定义了质量能满足临床需要的可用于数据交换的医学图像格式。dicom标准中通过(0020,0032)image position(patient)和(0020,0037)image orientation(patient)字段来规定患者的空间定位，(0028,0030)pixel spacing、(0018,0050)slice thickness等字段描述图像体素之间的间距。本发明利用这些字段信息构建仿射矩阵，实现lps坐标系和ijk坐标系之间的转换，进而利用插值法实现图像配准。
[0052]
在本发明的一个较佳实施例中，先以cest图像为参考图像，mr结构像为待配准图像为例，提供了一种基于装置成像坐标的mr图像配准方法，该方法用于实现mr结构像与cest图像之间的跨模态配准，具体是将作为待配准图像的mr结构像配准到作为参考图像的cest图像上，利用mr结构像的解剖信息来提取感兴趣区域，以便于在配准后的cest图像上进行后续的定量分析。进行配准的mr结构像与cest图像，需要是同一患者在同一次mr扫描时采集得到的，且两次序列采集之间无明显头动，即二者基于同一lps坐标系，满足跨模态配准的要求，由此方能采用本发明的方法进行配准。
[0053]
下面对该实施例中mr图像配准方法的其具体步骤进行详细描述，其具体过程分为s1～s4共四个步骤。
[0054]
s1、基于从cest图像的dicom信息中提取的旋转变换矩阵、平移变换矩阵和缩放变
换矩阵，计算将cest图像由图像坐标系变换到解剖学空间坐标系所需的第一仿射变换矩阵a1；同时，基于从mr结构像的dicom信息中提取的旋转变换矩阵、平移变换矩阵和缩放变换矩阵，计算将mr结构像由图像坐标系变换到解剖学空间坐标系所需的第二仿射变换矩阵a2。
[0055]
需注意的是，用于计算第一仿射变换矩阵a1和第二仿射变换矩阵a2所需的旋转变换矩阵、平移变换矩阵和缩放变换矩阵是不同的，需各自从cest图像和mr结构像的dicom信息中提取，但两者的提取和计算过程是一致的。下面将cest图像以及mr结构像分别记为目标图像，然后描述从目标图像中获取仿射变换矩阵的基本流程，仿射变换包括旋转、平移和缩放三个变换，具体步骤如下：
[0056]
s11、旋转变换矩阵计算：
[0057]
旋转变换可由dicom头文件中的(0020,0037)image orientation(patient)字段计算得到。该字段为一六元double数组[r
ix
,r
iy
,r
iz
,r
jx
,r
jy
,r
jz
],r
ix
,r
iy
,r
iz
分别代表行向量(i轴)与x,y,z轴夹角的余弦，r
jx
,r
jy
,r
jz
分别为列向量(j轴)与x,y,z轴夹角的余弦，另外k轴与x,y,z轴夹角的余弦[r
kx
,r
ky
,r
kz
]可由[r
ix
,r
iy
,r
iz
]和[r
jx
,r
jy
,r
jz
]的叉积计算得到。旋转变换矩阵表示如下：
[0058][0059]
s12、平移变换矩阵计算：
[0060]
平移变换可由dicom头文件中的(0020,0032)image position(patient)字段得到。该字段为三位double数组[t
x
,t
y
,t
z
]，用于表示当前图像的首个像素即ijk坐标的原点(左上角)在lps坐标体系下的坐标，单位为毫米。平移变换矩阵表示如下：
[0061][0062]
s13、缩放变换矩阵计算：
[0063]
缩放变换则由dicom头文件中的(0028,0030)pixelspacing配合(0018,0050)slicethickness或(0018,0088)spacingbetweenslices字段表示。其中，(0028,0030)pixelspacing为二元double数组[s
x
,s
y
]，用来表示相邻像素中心点之间的距离，即i,j轴的单位长度在lps坐标系中所占据的长度。(0018,0050)slicethickness或(0018,0088)spacingbetweenslices为一元double变量s
z
，用于表示相邻slice之间的距离，即z轴的单位长度在lps坐标系中所占据的长度。(0018,0050)slicethickness或(0018,0088)spacingbetweenslices两者任取其一即可，若存在(0018,0088)spacingbetweenslices则优先采用(0018,0088)spacingbetweenslices，若不存在(0018,0088)spacingbetweenslices则可以用(0018,0050)slicethickness代替。缩放变换矩阵表示如下：
[0064]
[0065]
s14、仿射变换矩阵计算：
[0066]
根据目标图像对应的旋转变换矩阵、平移变换矩阵和缩放变换矩阵，得到该目标图像的仿射变换矩阵a为：
[0067][0068]
需注意的是，上述目标图像可选自cest图像和mr结构像，若目标图像为cest图像，则a＝a1，即得到的仿射变换矩阵a为第一仿射变换矩阵a1，若目标图像为mr结构像，则a＝a2，即得到的仿射变换矩阵a为第二仿射变换矩阵a2。为了便于描述，可将上述两种仿射变换矩阵分别记为：
[0069][0070][0071]
当获取到两个仿射变换矩阵后，即可进行cest图像和mr结构像之间的后续坐标转换和配准。
[0072]
s2、根据第一仿射变换矩阵a1，将cest图像中每个体素在图像坐标系下的坐标[i1,j1,k1]转换为解剖学空间坐标系下的坐标[x,y,z]：
[0073][0074]
即
[0075][0076]
其中i1＝0:rows
‑
1，j1＝0:columns
‑
1，k1＝0:slices
‑
1，row、columns、slices分别为cest图像的行数、列数、层数。[x,y,z]为cest图像每个体素对应的lps坐标。
[0077]
s3、由于mr结构像和cest图像基于同一lps坐标系，因此基于mr结构像和cest图像的解剖学空间坐标系的统一性，利用第二仿射变换矩阵a2的逆变换，将cest图像中每个体素在解剖学空间坐标系下的坐标[x,y,z]转换为该体素在mr结构像的图像坐标系下的坐标[i2,j2,k2]：
[0078]
[0079]
即
[0080][0081]
其中i1＝0:rows
‑
1，j1＝0:columns
‑
1，k1＝0:slices
‑
1，[i2,j2,k2]为cest图像每个体素在mr结构像中对应的ijk坐标。
[0082]
s4、针对cest图像中每个体素经过转换后的坐标[i2,j2,k2]，通过插值方法利用mr结构像中该坐标周边的体素灰度值计算该坐标处的灰度值，从而得到与cest图像中每个体素一一对应的待配准图像，实现mr结构像在cest图像上的跨模态配准。
[0083]
上述插值的目的是为了获取mr结构像中坐标[i2,j2,k2]处的灰度值，因为ijk坐标[i2,j2,k2]是经过了两次转换得到的，而mr结构像中的体素坐标在三个ijk维度上并非完全连续的，因此[i2,j2,k2]位置大概率无法完全映射至mr结构像中的某一体素，需要通过插值方法计算其所在位置的体素灰度值。当cest图像中每一个体素处均通过插值方法得到其在mr结构像中的体素灰度值后，相当于将mr结构像配准到了cest图像上。具体的插值方法可根据实际需要调整，例如可采用线性插值或三次样条插值法等，理论上任意适用于mr结构像中插值的方法均可。
[0084]
在本实施例中，作为待配准图像的mr结构像是三维图像，而作为参考图像为cest图像既可以是三维图像也可以是二维图像。当cest图像是二维图像时，mr结构像的配准过程实际上是形成了一层与cest图像中体素一一对应的二维mr结构像，该mr结构像反应了cest图像中每一个体素在对应的三维mr结构像中的灰度值。当cest图像时三维图像时，本质上是对每一层二维的cest形成一层对应的二维mr结构像。
[0085]
需要注意的是，在上述配准过程中，(0020,0032)image position(patient)指向了图像首个像素的中心坐标，因此对于作为参考图像的cest图像层厚slicethickness大于1mm的情况，由于其层厚较大，单独用一层基于(0020,0032)image position(patient)得到的mr结构像作为cest图像对应的解剖结构图像的解剖信息，可能存在不准确的问题。因此，为提高配准结果的准确性，若cest图像的层厚slicethickness大于1mm，则可以采用多轮配准并取均值的方法得到最终配准结果，具体步骤如下：
[0086]
首先，针对从cest图像的dicom信息中提取的原始平移变换矩阵将其沿图像法向量方向移动n次，n为以mm为单位的cest图像层厚的四舍五入值，即n＝round(slicethickness)，round()为四舍五入函数，每次移动步长为(slicethickness)/n mm。将原始平移变换矩阵移动后得到的新的平移变换矩阵记为由此n次移动得到n个新的平移变换矩阵为：
[0087]
[0088]
然后，依次以每一个新的平移变换矩阵作为从cest图像的dicom信息中提取的平移变换矩阵，重复s1～s4，得到一张配准至cest图像上的mr结构像，一共得到的n张mr结构像，对这n张mr结构像进行平均，平均结果即可作为最终配准cest图像上的mr结构像。
[0089]
需说明的是，在上述每一次重复s1～s4时，从cest图像的dicom信息中提取的原始平移变换矩阵需要用移动后得到的新的平移变换矩阵代替，进而用于构建第一仿射变换矩阵a1，而从mr结构像的dicom信息中提取的原始平移变换矩阵则不需要移动，直接用原始平移变换矩阵构建第二仿射变换矩阵a2即可。上述对于平移变换矩阵的移动，使得能够将多层不同法向位置的mr结构像均配准至cest图像上，其平均信息能更准确反映厚度较大的cest图像对应的解剖结构图像的解剖信息。
[0090]
另外，虽然上述参与配准的mr结构像与cest图像是同一患者在同一次mr扫描时采集得到的，原则上保证两次序列采集之间无明显头动，但在实际情况中，依然不可避免地存在轻微的头动现象。因此，如果cest图像和mr结构像采集期间存在轻微头动，本实施例中需在上述s4步骤之后，利用flirt或spm等工具将配准得到的mr结构像与参考图像再次进行配准，以校正cest图像和mr结构像采集期间的轻微头动，提高配准的准确性。当然，假如cest图像和mr结构像采集期间的头动可以忽略，则该步骤亦可不做。
[0091]
上述配准方法的意义在于：mr结构像的灰度值体现了解剖结构图像的解剖信息，可以利用配准到cest图像上的mr结构像来勾画感兴趣区域(roi)，并基于两种图像经过配准后的空间一致性在cest图像上也得到对应的roi，便于在cest图像上进行定量分析，解决cest图像上感兴趣区域难以准确提取的缺陷。
[0092]
因此，在本发明的另一较佳实施例中，还基于上述实施例中的mr图像配准方法，进一步提供了一种基于图谱的cest图像定量分析方法，该方法如下：
[0093]
首先对采集的cest数据进行配准和校正处理得到cest图像，然后再按照上述s1～s4所述的mr图像配准方法，将mr结构像配准至cest图像上；再基于配准后的mr结构像确定感兴趣区域roi并将感兴趣区域映射至cest图像上，从而计算该感兴趣区域内的cest量化参数。
[0094]
需说明的是，cest量化参数可以在勾画roi后针对roi内的体素进行计算，也可以在cest数据的配准和校正处理后针对整幅cest图像中每个体素预先进行计算然后再根据roi进行统计，对此不作限定。
[0095]
还需要特别说明的是，本发明的cest图像定量分析方法中，目的仅仅是根据实际需要得到指定的roi内的cest量化参数，但该参数与相应的疾病之间并不存在必然关联。在部分应用中，这些cest量化参数可以作为科学研究等非诊断治疗之用。
[0096]
为获取cest图像，对cest数据的具体处理方式可参照现有技术进行。一般而言，采集的cest数据需首先将cest图像的所有帧配准到目标代谢产物集团所对应的帧上，然后进行b0场不均匀性校正，方可得到cest图像。本实施例中，下面提供了一种对cest数据进行后处理，得到已计算cest指标的cest量化图的具体方法，该方法按照s11～s13的流程进行：
[0097]
s11：为校正潜在的运动伪影，利用flirt或spm等工具将cest图像的所有帧配准到目标代谢产物集团所对应的帧上，配准包括刚体变换、相关比例代价函数和三线性插值法
等。
[0098]
s12：使用wassr(water saturation shift referencing)技术或梯度双回波序列采集b0图校正b0场不均匀性。
[0099]
s13：计算cest量化图，本发明可计算以下常见的cest指标：
[0100]
(1)化学交换饱和转移比率(cest ratio,cestr),利用磁化转移率不对称性(magnetization transfer ratio asymmetry,mtr
asym
)计算得到。
[0101][0102]
s0是没有射频饱和时测量的图像信号强度，s
ref
和s
lab
分别为参考帧图像信号强度和目标帧图像信号强度，z
ref
和z
lab
是对应的z谱信号强度。
[0103]
(2)利用参考值归一化cestr(cestr normalized with the reference value，cestr
nr
)
[0104][0105]
(3)基于倒z谱的磁化转移率(inverse z
‑
spectrum
‑
based magnetizationtransfer ratio,mtr
rex
),校正了直接水饱和(direct water saturation,ds)和半固态磁化转移对cest信号的影响。
[0106][0107]
(4)表观交换弛豫(apparent exchange
‑
related relaxation，arex),可消除纵向弛豫(t1)对cest信号的影响。
[0108][0109]
上述具体的cest指标仅作为参考，可根据实际需要进行选择一种或多种。后续将mr结构像配准至cest图像上后，即可在mr结构像上提取感兴趣区域并根据cest量化图计算该区域的cest量化参数。例如，如需计算表观交换弛豫arex，需要t1 mapping图像，该t1 mapping图像需与cest图像对齐，因此先按照s1～s4的方法将t1 mapping图像配准到cest图像上，再计算roi内的arex map。在mr结构像上提取感兴趣区域可以采用手动勾画，也可使用自动分割工具自动分割脑区，计算各个脑区内cest量化参数的均值。
[0110]
另外，在上述实施例中的mr图像配准方法以及基于图谱的cest图像定量分析方法中，虽然都是以cest图像为参考图像，mr结构像为待配准图像，但是两者也可以相互调换，即以mr结构像为参考图像，以cest图像为待配准图像，其具体做法是类似的，只需要对换两种图像即可。
[0111]
另外，在本发明的另一实施例中，还可以进一步提供一种基于装置成像坐标的mr图像配准装置，其包括存储器和处理器；
[0112]
所述存储器，用于存储计算机程序；
[0113]
所述处理器，用于当执行所述计算机程序时，实现如上述s1～s4所述的mr图像配准方法。
[0114]
同样的，在本发明的另一实施例中，还可以进一步提供一种基于图谱的cest图像
定量分析装置，其包括存储器和处理器；
[0115]
所述存储器，用于存储计算机程序；
[0116]
所述处理器，用于当执行所述计算机程序时，实现上述基于图谱的cest图像定量分析方法。
[0117]
需要注意的是，存储器可以包括随机存取存储器(random access memory，ram)，也可以包括非易失性存储器(non
‑
volatile memory，nvm)，例如至少一个磁盘存储器。上述的处理器可以是通用处理器，包括中央处理器(centralprocessing unit，cpu)、网络处理器(network processor，np)等；还可以是数字信号处理器(digital signal processing，dsp)、专用集成电路(application specificintegrated circuit，asic)、现场可编程门阵列(field
‑
programmable gate array，fpga)或者其他可编程逻辑器件、分立门或者晶体管逻辑器件、分立硬件组件。当然，还装置中还应当具有实现程序运行的必要组件，例如电源、通信总线等等。
[0118]
下面将上述实施例中的基于图谱的cest图像定量分析方法应用于一个具体实例中，以进一步展示本发明所能够实现的技术效果，以便本领域技术人员更好地理解本发明的实质。
[0119]
实施例
[0120]
1、mri实验
[0121]
用64通道头部线圈的3特斯拉siemens prisma扫描仪(erlangen,germany)扫描了2位健康志愿者和一名脑肿瘤患者的大脑。
[0122]
所有志愿者均采用3d space
‑
cest序列，视野(fov)＝212
×
212
×
201mm3,分辨率＝2.8
×
2.8
×
2.8mm3,重复时间(tr)＝3seconds,回波时间(te)＝17ms,加速因子(nominal turbo factor)＝140,信号平均次数(number of signalaverages，nsa)＝1.2,grappa加速因子＝2
×
2。施加7个cest饱和频率用于酰胺质子转移加权(aptw)成像，包括不饱和(s0)和
±
3ppm、
±
3.5ppm和
±
4ppm的饱和频率。b0图采用与space
‑
cest序列相同视野、方向和分辨率的3d gre序列，tr＝30ms，光谱预饱和，反演恢复脂肪抑制。gre序列在双回波模式下运行，te分别为4.92ms和9.84ms，使得水和脂肪旋转的相位一致。结构像包括flair(水抑制反转恢复)序列(t2_tirm_tra_dark
‑
fluid，tr＝8s；te＝81ms；fov＝220
×
220
×
130mm3；分辨率＝0.7
×
0.7
×
5mm3；反转时间ti＝2.4s)。t1wi
‑
3d
‑
mprage(t1加权三维预磁化快速梯度回波)序列(tr＝2.3s；te＝2.9ms；fov＝256
×
256
×
176mm3；分辨率＝1
×1×
1mm3；ti＝0.9s)。
[0123]
健康志愿者a采用t1wi
‑
3d
‑
mp2rage(t1加权三维双预磁化快速梯度回波)序列,tr＝5s,te＝3.59ms，fov＝230
×
216
×
132mm3；分辨率＝0.7
×
0.7
×
3mm3，获得t1 mapping图。另一名健康志愿者b在轴向上采用2d快速自旋回波(tse)cest序列，采集参数为tr＝3s，te＝7.2ms，fov＝212
×
186
×
5mm3，分辨率＝1.1
×
1.1
×
5mm3，使用“频谱预饱和反转恢复”(spir)脂肪抑制。施加一个不饱和频率和62个饱和频率，范围为
‑
6ppm到80ppm,频率间隔0.25ppm，在
±
3.5ppm做了多次平均用于提高apt加权(aptw)图像的snr(信噪比)。用于b0场校正的2d gre b0序列与2d cest序列使用相同的视野、方向和分辨率，tr＝30m，光谱预饱和，反演恢复脂肪抑制。gre序列在双回波模式下运行，te分别为4.92ms和9.84ms。
[0124]
2、mr结构像和cest图像的配准以及cest图像定量分析
[0125]
在对3d cest数据进行定量分析时，首先按照s11步骤用spm进行运动校正，再按上述步骤s1～s4将cest图像、b0图像、t1 mapping配准到mr结构像上，然后根据步骤s12～s13对配准后的cest图像进行处理，得到cest量化参数图。步骤s1～s4的mr图像配准过程如前所述，不再赘述，而本实例中3dcest数据的处理流程如下：
[0126]
s11：为校正研究中潜在的运动伪影，利用flirt或spm等工具将cest图像的所有帧配准到目标代谢产物集团所对应的帧上，实施例选用酰胺质子为研究对象，故将cest图像所有帧配准到 3.5ppm帧上，具体的配准参考帧根据研究对象的不同而不同。配准包括刚体变换、相关比例代价函数和三线性插值法等。
[0127]
s12：使用wassr(water saturation shift referencing)技术或梯度双回波序列采集b0图校正b0场不均匀性。
[0128]
首先，对cest和b0源图像进行空间插值，使其分辨率提高一倍。
[0129]
wassr技术通过施加足够小的功率以及持续时间足够短的rf脉冲，使磁化转移(magnetization transfer)和cest效应被忽略，而直接水饱和(directwater saturation,ds)起主要作用，通过该方法可以得到自由水信号的频率偏移信息，将自由水的频率校正到0ppm，从而提供b0场映射图。
[0130]
用双回波梯度序列采集了b0场映射图，将两个不同回波时间获得的b0相位图像之间的相位差与回波时间差4.92ms相除获得b0映射图，量化了每个体素b0场偏移映射，基于该偏移映射对b0场进行校正。该方法的采集时间短于wassr方法，而wassr方法提供了通常更优的b0场映射图，实施例选用双回波梯度序列采集b0映射图进行b0校正。
[0131]
s13：cest源图像首先用不饱和帧进行归一化，逐体素将z谱的 2.5～ 4.5ppm和
‑
2.5～
‑
4.5ppm的范围内线性插值至1hz的分辨率。然后计算cest量化图，本发明可计算以下常见的cest指标：
[0132]
(1)化学交换饱和转移比率(cest ratio,cestr),利用磁化转移率不对称性(magnetization transfer ratio asymmetry,mtr
asym
)计算得到。
[0133][0134]
s0是没有射频饱和时测量的图像信号强度，s
ref
和s
lab
分别为参考帧图像信号强度和目标帧图像信号强度，z
ref
和z
lab
是对应的z谱信号强度。实施例选用酰胺质子为研究对象，cestr( 3.5ppm)也称aptw(酰胺质子加权)，δω
lab
＝ 3.5ppm，δω
ref
＝
‑
3.5ppm。
[0135]
(2)利用参考值归一化cestr(cestr normalized with the reference value，cestr
nr
)
[0136][0137]
(3)基于倒z谱的磁化转移率(inverse z
‑
spectrum
‑
based magnetizationtransfer ratio,mtr
rex
),校正了直接水饱和(direct water saturation,ds)和半固态磁化转移对cest信号的影响。
[0138]
[0139]
(4)表观交换弛豫(apparent exchange
‑
related relaxation，arex),可消除纵向弛豫(t1)对cest信号的影响。
[0140][0141]
图2为健康志愿者a的t1wi
‑
mprage图(a)，t1wi
‑
mprage图自动分割后的label图(b)，头动校正后以t1wi
‑
mprage图为参考图像进行配准得到的cest不饱和帧源图像(c)，以t1wi
‑
mprage图为参考图像进行配准得到的t1mapping图(d)，以及用配准后的cest图像和t1 mapping图计算的aptw图(e)，cestr
nr
图(f)，mtr
rex
图(g)和arex图(h)。列为不同层面。计算各脑区内cest量化参数的均值，健康志愿者a的各脑区cest量化参数均值如下表所示：
[0142]
[0143][0144]
图3为脑肿瘤患者的t1wi
‑
mprage图(a)，头动校正后以t1wi
‑
mprage图为参考图像进行配准得到的cest不饱和帧源图像(b)，用配准后的cest图像计算的aptw图(c)，cestr
nr
图(d)，mtr
rex
图(e)。列为不同层面。aptw图、cestr
nr
图，mtr
rex
图上的肿瘤区域的信号都高于正常脑组织区域。
[0145]
另外，在对2d cest数据进行定量分析的流程与前述一致，仅需将参考图像改为2d cest图像。图4为健康志愿者b头动校正后的2d cest不饱和帧源图像(a)，以2d cest图像为参考图像进行配准后的flair图(b)，t1wi
‑
mprage图(c)和t1wi
‑
mprage
‑
label图(d)，(e)为aptw图。
[0146]
以上所述的实施例只是本发明的一种较佳的方案，然其并非用以限制本发明。有关技术领域的普通技术人员，在不脱离本发明的精神和范围的情况下，还可以做出各种变化和变型。因此凡采取等同替换或等效变换的方式所获得的技术方案，均落在本发明的保护范围内。