CN109389651B - 磁共振化学位移编码成像方法与装置 - Google Patents

磁共振化学位移编码成像方法与装置 Download PDF

Info

Publication number
CN109389651B
CN109389651B CN201710674824.0A CN201710674824A CN109389651B CN 109389651 B CN109389651 B CN 109389651B CN 201710674824 A CN201710674824 A CN 201710674824A CN 109389651 B CN109389651 B CN 109389651B
Authority
CN
China
Prior art keywords
field
field map
value
image
map
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201710674824.0A
Other languages
English (en)
Other versions
CN109389651A (zh
Inventor
郑海荣
刘新
邹超
程传力
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Shenzhen Institute of Advanced Technology of CAS
Original Assignee
Shenzhen Institute of Advanced Technology of CAS
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Shenzhen Institute of Advanced Technology of CAS filed Critical Shenzhen Institute of Advanced Technology of CAS
Priority to CN201710674824.0A priority Critical patent/CN109389651B/zh
Publication of CN109389651A publication Critical patent/CN109389651A/zh
Application granted granted Critical
Publication of CN109389651B publication Critical patent/CN109389651B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T11/002D [Two Dimensional] image generation
    • G06T11/003Reconstruction from projections, e.g. tomography
    • G06T11/008Specific post-processing after tomographic reconstruction, e.g. voxelisation, metal artifact correction
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01RMEASURING ELECTRIC VARIABLES; MEASURING MAGNETIC VARIABLES
    • G01R33/00Arrangements or instruments for measuring magnetic variables
    • G01R33/20Arrangements or instruments for measuring magnetic variables involving magnetic resonance
    • G01R33/44Arrangements or instruments for measuring magnetic variables involving magnetic resonance using nuclear magnetic resonance [NMR]
    • G01R33/48NMR imaging systems
    • G01R33/54Signal processing systems, e.g. using pulse sequences ; Generation or control of pulse sequences; Operator console
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/10Image acquisition modality
    • G06T2207/10072Tomographic images
    • G06T2207/10088Magnetic resonance imaging [MRI]

Landscapes

  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Signal Processing (AREA)
  • High Energy & Nuclear Physics (AREA)
  • Condensed Matter Physics & Semiconductors (AREA)
  • Magnetic Resonance Imaging Apparatus (AREA)

Abstract

本发明涉及磁共振化学位移编码成像方法与装置以及执行所述方法的各步骤的计算机可读介质。所述方法包括:(1)求得二维原始磁共振图像的两组场图候选解;(2)在不同的尺度下,将前述两组场图候选解分别分割为不交叠的子块;(3)在每一个尺度下分别采用局部迭代场图提取算法,得到相应的原始图像空间下的场图;(4)将各个尺度下得到的场图合并,保留具有相同场图值的像素点的场图值,其余像素点的场图值置0;(5)利用局部迭代场图提取算法对置0像素点的场图值进行选取,得到最终场图;(6)对最终场图进行平滑处理,分离原始磁共振图像,得到准确的各待分离成分的图。本发明所述方法能实现正确的场图估计和水脂分离。

Description

磁共振化学位移编码成像方法与装置
技术领域
本发明涉及磁共振成像技术领域,尤其涉及基于多尺度局部迭代场图提取算法的磁共振化学位移编码成像方法与装置。
背景技术
磁共振化学位移编码成像是一种基于组织中各成分之间化学位移差异的成像方法。化学位移编码成像在不同的回波时间采集信号,使得每一组信号中各成分之间存在相位差。临床上最常用的化学位移编码成像是水脂分离成像,它主要用于脂肪抑制和脂肪定量等应用中。
多回波化学位移编码成像假设组织信号由多种不同成分(本发明中以水和脂肪分离成像为例)被同时激发得到,并已知各成分相对于水的化学位移,然后拟合信号与各成分、B0场图和回波时间(echo time,TE)之间的数学模型。实际应用最多的化学位移编码成像是水脂分离成像,即通过多回波信号分离水和脂肪两种成份,以下均以水脂分离成像为例说明两点化学位移编码成像方法。
1984年,Dixon[1]最早提出了两点水脂分离技术,该技术假设场图分布均匀,选取两个特定的TE分别采集水和脂肪质子相位同向(In-phase,IP)和反向(opposed-phase,OP)的图像,通过简单的相加/相减操作得到水图和脂肪图。为了解决B0场不均匀造成的影响,2006年,Xiang[2]提出了一种新的两点水脂分离技术,该技术不再限制第二幅采集图像中水和脂肪相位反向,而是采集水脂相位部分反向(partially-opposed-phase,POP)图像,利用余弦定理得到两组水和脂肪的幅值图像,在此基础上解到两个场图候选值然后将原始图像分割为多个不交叠的子块,根据子块内场图分布的光滑特性,选取出每一个子块内具有最大场图相似性的场图集合作为初始场图,利用区域迭代场图提取(Regional IterativePhasor Extraction,RIPE)算法从场图候选解中将整个原始图像空间内的场图估计出来,最后利用最小二乘法得到分离的水和脂肪图像。2011年,Eggers等[3]在Xiang方法的基础上提出了新的两点分离技术,该技术也不再限制第一幅采集图像中水和脂肪相位同向,只是限制第一幅采集图像中的水脂相位差θ1和第二幅图像中的水脂相位差θ2之间的差或和不能等于2π或其整数倍。该技术借用Xiang的方法对场图进行估计。
多回波化学位移编码成像面临的最大挑战是由主磁场B0不均匀性(本发明称之为场图,B0均匀时,场图连续缓慢变化)造成的局部分反,因此该技术的分离效果依赖于B0场图的估计。
前述现有算法存在其缺点,文献[2-3]在不同的子块得到各自的场图集合,将其直接拼接后作为初始场图进行区域迭代场图提取,这样做的缺陷是在场图提取过程中容易出现选取矛盾,即从不同方向进行场图选取时得到不同的选取结果,造成最后的场图中存在跳变。
发明内容
针对上述缺点,本发明提出一种基于多尺度局部迭代场图提取(multi-scaleRIPE)算法的磁共振化学位移编码成像方法,以实现正确的场图估计和待分离成分的分离(例如,水脂分离)。
本发明所述的磁共振化学位移编码成像方法,所述方法包括:
(1)求得二维原始磁共振图像的两组场图候选解(又称场图候选值);
(2)在不同的分块尺度下,将前述两组场图候选解分别分割为不交叠的子块;
(3)在每一个分块尺度下分别采用局部迭代场图提取算法,得到相应的原始图像空间下的场图;
(4)将各个分块尺度下得到的场图合并,保留具有相同场图值的像素点的场图值,其余像素点的场图值置0;
(5)利用局部迭代场图提取算法对置0像素点的场图进行选取,得到最终场图;
(6)对最终场图进行平滑处理,分离原始磁共振图像(即,利用平滑处理后的场图对原始磁共振图像进行分离),得到准确的各待分离成分(例如水和脂肪,或者水与另一待分离成分)的图。
在一实施方式中,所述待分离成分为水和脂肪。
在一实施方式中,所述待分离成分为水和另一成分,所述成分可以是硅胶或其它任意一种和水具有不同化学位移的成分;也可以分离任意两种具有不同化学位移的成分。
在一实施方式中,所述的步骤(1)根据以下执行:
描述二维原始磁共振图像中每个像素点的场图值拟合误差曲线,并根据该场图值拟合误差曲线寻找使之达到局部极小值的所有场图值,构成该像素点的两组场图候选解。
在一子实施方式中,所述的步骤(1)可根据以下执行:根据式(5)描述二维原始磁共振图像中每个像素点的场图值拟合误差曲线(phasor-error谱),
p=argminerr(p)=argmin||(I-AA+)S||2 (5)
其中p是相位矢量(或相位因子),也称为场图值;I是4×4单位矩阵;A+=(ATA)-1AT,上标“+”表示取矩阵伪逆,上标“T”表示矩阵的转置;err(p)描述了相位矢量p和拟合误差之间的关系,表示为phasor-error谱,其定义为不同相位矢量值下对应的模型拟合误差;S如下面式(4)中所定义;
根据该phasor-error谱寻找使式(5)表征的二范数达到局部极小值的所有场图值p,构成该像素点的场图候选解,个数为1个或2个,以符号p1和p2表示,场图值数目为1时,p1=p2
在一实施方式中,所述式(5)通过以下获得:
采用以下式(1)的含有水和另一成分(例如脂肪)的两种成分的MR信号模型:
Figure BDA0001373998240000031
其中TE为回波时间,Sn是在回波时间TEn下的信号强度;ρw和ρf是水和另一成分(例如脂肪)的信号强度值;水的化学位移设定为0,fF是另一成分相对于水的化学位移,为已知量,例如脂肪,为-3.5ppm;fB是静磁场B0的场图值,i表示虚数;
利用式(2)定义以下参数c1、c2、q和p:
Figure BDA0001373998240000041
其中ΔTE=TE1-TE2,q代表TE1时间内由于fB造成的水的相位,p代表ΔTE时间内由于fB造成的水的相位累积,即如上所定义的相位矢量或场图值,其角度与fB成正比,(角度)范围是[-ππ];
其中利用式(2)的p,由式(1)和(2)解得以下式(3)表示的q:
Figure BDA0001373998240000042
通过利用式(2)、(3),将式(1)改写成以下式(4)的矩阵形式:
Figure BDA0001373998240000043
其中
Figure BDA0001373998240000044
是4×1矩阵,
Figure BDA0001373998240000045
Figure BDA0001373998240000046
S1、S2分别独立地为在回波时间TE1、TE2下一个像素点的信号强度,
Figure BDA0001373998240000047
Re代表取复数实部,Im代表取复数虚部;
通过式(4),根据变量投影方法,使用最小二乘法解得式(5)表示的p:
p=argminerr(p)=argmin||(I-AA+)S||2 (5)。
在一实施方式中,所述的步骤(2)根据以下执行:
根据预先定义的分块尺度值将待分割的二维原始磁共振图像进行分块,得到的每一组图像中均包含数目不等的子块,每一个子块的大小对应其分块尺度。
在一实施方式中,所述分块尺度值为整数,选自1至16中的一个或多个数值。例如,所述分块尺度值可以是任选的4个数值:6、8、10和12。
在一实施方式中,所述步骤(3)根据以下执行:
(3a)场图的初始化
对每一个分块尺度下的每一个子块,求得子块场图候选解,利用每一个子块场图候选解为初始值,对子块内所有像素点的两组场图候选解进行选取,选取规则为像素点场图候选解与初始值差值的幅值距离最小,符合所述规则的像素点场图候选解组成一个场图值集合,未被选取的像素点场图候选解组成另一个场图值集合;计算每一个场图候选解对应场图值集合的矢量和;比较所述矢量和的幅值,选取最大幅值对应的场图值集合;
对每一个分块尺度下每一个子块均进行上述场图初始化,将每一个分块尺度下所有子块的场图值集合合并在一起,即得到该分块尺度下的初始场图值集合;
(3b)场图的更新
对上述初始场图值集合进行低通滤波,且对于下面每一步迭代后得到的场图值集合均进行低通滤波处理;然后以每一个像素点的场图值为初始值,对图像内所有像素点的两组场图候选解再次进行选取,选取规则为像素点场图候选解与初始值差值的幅值距离最小;
并记录场图值发生改变的像素点的个数N;
(3c)迭代更新场图
重复步骤(3b),直至N小于一定的数值或在数次迭代中均固定在某一数值则迭代停止,得到一组场图;
对应于多组分块尺度,经过步骤(3a)、(3b)和(3c)得到多组的场图。
在一子实施方式中,所述步骤(3)可根据以下执行:
(3a)场图的初始化
每一个分块尺度下的每一子块内所有像素点作为整体根据式(6)和(7)描述场图值拟合误差曲线,根据该场图值拟合误差曲线寻找使式(6)达到局部极小值的场图值,所述场图值为该子块的场图候选解,数目最多为2个,该子块的场图候选解以符号pu和pv表示,当场图值数目为1时,pu=pv
Figure BDA0001373998240000066
式(6)由式(5)改写得到,其中,I是4×4单位矩阵;A+=(ATA)-1AT,上标“+”表示取矩阵伪逆,上标“T”代表矩阵的转置,
Figure BDA0001373998240000067
表示点乘;err(p)描述了场图值p和拟合误差之间的关系,表示为场图值拟合误差曲线,其定义为不同场图值p下对应的模型拟合误差,
Figure BDA0001373998240000061
式(7)中,
Figure BDA0001373998240000062
Figure BDA0001373998240000063
分别代表一组像素点的信号强度,定义如式(4)中表述,
Figure BDA0001373998240000064
Figure BDA0001373998240000065
都是行向量,行向量中每一个元素代表一个像素点的信号强度;
以上述子块的一个场图候选解为该子块的初始值,对该子块内所有像素点的两个场图候选解(p1,p2)进行选取,选取的规则是:计算该子块内每个像素点的两个场图候选解与子块初始值之间的差值,并对差值取绝对值,较小绝对值对应的场图候选解即为该像素点的场图值,以此方法对子块内所有像素点的场图值进行选取,并计算这些场图值的矢量和,同时计算所有像素点未被选取的场图候选解的矢量和,比较上述两组矢量和的幅值,保留较大值对应的幅值,舍弃另一个幅值;
同样地,以该子块另一个场图候选解为初始值,对该子块内所有像素点的场图值进行选取,也得到两组场图值矢量和的幅值,比较这两组幅值,保留较大值对应的幅值,舍弃另一个幅值;
比较该子块2个场图候选解所得到的保留幅值,并保留较大值的保留幅值对应的子块场图候选解,舍弃另一组较小值对应的场图候选解;根据保留的子块场图候选解对应的两组场图值的矢量和得到加权因子,进而得到该子块内所有像素点的场图值集合P;
或者,根据下文描述的具体方法对该子块内所有像素点的两个场图候选解(p1,p2)进行选取;
对每一个分块尺度下每一个子块均进行上述场图初始化,将每一个分块尺度下所有子块的场图值集合P合并在一起,即得到该分块尺度下的初始场图值集合P0
(3b)场图的更新
对上述初始场图值集合P0进行低通滤波,且对于下面每一步迭代后得到的场图值集合均进行低通滤波处理;
假定迭代后得到的场图值集合为Pn,则滤波后得到Pn s,然后以每一个像素点的场图值pn s为初始值,对图像内所有像素点的两个场图候选解(p1,p2)再次进行选取,选取方法类似于步骤(3a)场图的初始化中对子块内所有像素点的两个场图候选解(p1,p2)的选取方法,
并记录场图值发生改变的像素点的个数N;
(3c)迭代更新场图
重复步骤(3b),直至N小于一定的数值或在数次迭代中均固定在某一数值则迭代停止,得到一组场图;
对应于多组分块尺度,经过步骤(3a)、(3b)和(3c)得到多组的场图。
在一实施方式中,在前述步骤(3a)中,在根据式(6)和(7)描述每一个分块尺度下的每一子块内所有像素点整体的场图值拟合误差曲线时,根据该场图值拟合误差曲线寻找使式(6)达到局部极小值的场图值时,设置以下阈值排除一定数目的场图值,使之数目为1个或2个,阈值设定为:
Amin+α*(Amax-Amin),
其中Amax和Amin分别表示场图值拟合误差曲线的最大值和最小值,α的选取要排除拟合误差大的场图值,α取0.3-0.8。在一实施方式中,α取0.4。
在一实施方式中,在前述步骤(3a)中,对每一子块内所有像素点的两个场图候选解(p1,p2)进行选取的方法具体为:
针对一个子块,以该子块的场图值pu和pv为初始值,分别对该子块内所有像素点的两个场图候选解(p1,p2)进行选取,以pu为例,选取方法如下:
Figure BDA0001373998240000071
这样对应pu,子块内所有像素点都选取了相应的场图值,组成一个集合Pu,余下的场图值也组成一个集合Puc
同样地,对应pv,存在集合Pv和Pvc
分别计算上述四个集合的矢量和的幅值,记为Mu、Muc、Mv和Mvc
取{Mu,Muc}最大值,记为Umax,取{Mv,Mvc}最大值,记为Vmax,将Umax和Vmax进行比较,舍弃较小数值对应的数据集合;假设保留的集合为Pu和Puc,且Mu大于Muc;根据保留的子块场图候选解对应的两组矢量和幅值计算得到一个加权因子C,定义如下:
Figure BDA0001373998240000081
则该子块所有像素点对应的场图值集合为P=Pu*C;如果Muc大于Mu,则P=Puc*C。
在一实施方式中,在步骤(3b)中,以每一个像素点的场图值pn s为初始值,对图像内所有像素点的两个场图候选解(p1,p2)再次进行选取的方法类似于步骤(3a)中的式(8),即采用如下式(10)进行选取:
Figure BDA0001373998240000082
在一实施方式中,在步骤(3c)中,N为像素点个数的8%-10%。
在一实施方式中,所述步骤(4)根据以下执行:
(4a)场图合并
将上述多组场图进行合并得到一幅场图,具体实现方式是:如果像素点在所有的场图中具有相同的数值,则该保留该像素点的场图值;如果像素点在所有的场图中具有不同的数值,则该像素点的场图值重新置0。
所述重新置0是指重新设置为未知或未选取状态。
在一实施方式中,所述步骤(5)根据以下执行:
对合并的场图中置0像素点的场图值重复步骤(3b)和(3c),得到最终的场图。
在一实施方式中,所述步骤(6)根据以下执行:
利用得到的最终场图,进行平滑处理,然后计算得到水图和另一成分(例如脂肪)的图。
在一子实施方式中,所述步骤(6)可根据以下执行:利用得到的最终场图,进行平滑处理,然后利用下式计算得到水图和另一成分(例如脂肪)的图:
Figure BDA0001373998240000091
其中ρw、ρf、A+和S如上所定义,根据ρw和ρf得到的水图和另一成分(例如脂肪)的图。
术语“场图值拟合误差曲线”或“phasor-error谱”,是指由两点水脂分离信号模型推导出了单个像素点和子块内所有像素点为一个整体的phasor-error谱。
术语“场图的初始化”,是指根据子块phasor-error谱得到的场图值对该子块内所有像素的场图值进行了初始化。
在一实施方式中,当待分离成分为水和其他成分时,只需将式(1)中的脂肪相对于水的化学位移fF替换该其他成分对应的化学位移值,即可分离水图和该其他成分图。
此外,本发明还提供一种计算机可读介质,该计算机可读介质具有存储在其中的程序,该程序是计算机可执行的,以使计算机执行上述各种步骤的处理。
本发明还提供一种利用上述方法的磁共振化学位移编码成像装置。该装置包含了不同模块用于实现上述方法中提及的各个步骤。
本发明的振化学位移编码成像装置可以包括:场图候选解求解模块,其从二维原始磁共振图像求得该图像的两组场图候选解;场图候选解分割模块,其在预先定义的不同的分块尺度下,将在场图候选解求解模块中求得的前述两组场图候选解分别分割为不交叠的子块;局部迭代场图提取模块,其在场图候选解分割模块定义的每一个分块尺度下分别采用局部迭代场图提取算法处理在场图候选解分割模块中获得的场图候选解子块,得到相应的原始图像空间下的场图;场图合并模块,其将在局部迭代场图提取模块中各个分块尺度下得到的场图合并,保留具有相同场图值的像素点的场图值,其余像素点的场图值置0;最终场图获取模块,其利用局部迭代场图提取算法对在场图合并模块中置0像素点的场图值进行选取,得到最终场图;图像分离模块,其对在最终场图获取模块中得到的最终场图进行平滑处理,分离原始磁共振图像,得到准确的各待分离成分的图。在一子实施方式中,所述待分离成分为水和另一成分。
在一实施方式中,所述场图候选解求解模块还包括:描述各像素点场图值拟合误差曲线模块,其描述二维原始磁共振图像中每个像素点的场图值拟合误差曲线;以及获取各像素点场图候选解模块,其根据所述场图值拟合误差曲线寻找使之达到局部极小值的所有场图值。
在一子实施方式中,所述场图候选解求解模块可包括:描述各像素点场图值拟合误差曲线模块,其根据权利要求3中所采用的式(5)描述二维原始磁共振图像中每个像素点的场图值拟合误差曲线;以及获取各像素点场图候选解模块,其根据在描述各像素点场图值拟合误差曲线模块中获得的所述场图值拟合误差曲线寻找使式(5)达到局部极小值的所有场图值p,构成该像素点的场图候选解,个数为1个或2个,以符号p1和p2表示,其中场图值数目为1时,p1=p2
在一实施方式中,所述场图候选解分割模块根据预先定义的分块尺度值将二维原始磁共振图像进行分块,得到的每一组图像中均包含数目不等的子块,每一个子块的大小对应其分块尺度。在一实施方式中,所述分块尺度值为整数,选自1至16中的一个或多个数值。在另一实施方式中,所述分块尺度值为6、8、10和12。
在一实施方式中,所述局部迭代场图提取模块包括:场图初始化模块,其实现前面提及的本发明所述的磁共振化学位移编码成像方法中所述的(3a)场图的初始化步骤;场图更新模块,其实现前面提及的(3b)场图的更新步骤;以及迭代更新场图模块,其实现前面提及的(3c)迭代更新场图步骤。
在一实施方式中,所述场图合并模块通过以下实现方式将在局部迭代场图提取模块中获得的多组场图进行合并得到一幅场图:如果像素点在所有的场图中具有相同的数值,则该保留该像素点的场图值;如果像素点在所有的场图中具有不同的数值,则该像素点的场图值重新置0。
在一实施方式中,所述最终场图获取模块对在场图合并模块中合并的场图中置0像素点的场图值重复步骤(3b)和(3c),得到最终的场图。
在一子实施方式中,所述图像分离模块对在最终场图获取模块中得到的最终场图进行平滑处理,然后(可选地根据前面提及的式(11))计算得到水图和另一成分的图。
前述任一实施方式或子实施方式中,所述另一成分为脂肪。
本发明与现有技术相比具有以下优点:
本发明提出的两点磁共振化学位移编码成像技术不对TE(回波时间)施加任何的限制,达到更快的图像采集速度,使其适用于一些快速成像的医学应用中,如动态成像和腹部屏气成像等;应用多尺度的局部迭代场图提取算法可独立得到多组场图,然后利用自检验机制将其进行合并,保证了场图估计的正确性。
附图说明
图1为本发明提出的基于多尺度局部迭代场图提取算法的磁共振化学位移编码成像方法应用于两点水脂分离算法的流程图。
图2显示根据本发明所述方法对二维原始磁共振图像进行处理,最终分离得到水图和脂肪图;其中图中K=6、8、10和12分别表示一个尺度,在每一个尺度下独立地完成区域迭代场图提取算法(RIPE)。
图3显示了本发明所述方法在体组织实验的颈部的分离结果图;(a)颈部冠状位分离结果;(b)颈部矢状位分离结果,图像从左向右依次为分离脂肪图和分离水图。
具体实施方式
为使本发明解决的技术问题、采用的技术方案和达到的技术效果更加清楚,下面结合附图和实施例对本发明作进一步的详细说明。可以理解的是,此处所描述的具体实施例仅仅用于解释本发明,而非对本发明的限定。另外还需要说明的是,为了便于描述,附图中仅示出了与本发明相关的部分而非全部内容。
通过以下获得本发明所述方法的式(5):
含有水和脂肪的两种成分的MR信号模型为:
Figure BDA0001373998240000121
其中TE为回波时间,Sn是在回波时间TEn下的信号强度;ρw和ρf是水和脂肪的信号强度值;水的化学位移设定为0,fF是脂肪相对于水的化学位移,为-3.5ppm;fB是静磁场B0的场图值,i表示虚数。
利用式(2)定义参数c1、c2、q和p如下:
Figure BDA0001373998240000122
Figure BDA0001373998240000123
Figure BDA0001373998240000124
Figure BDA0001373998240000125
其中ΔTE=TE1-TE2,q代表TE1时间内由于fB造成的水的相位,p代表ΔTE时间内由于fB造成的水的相位累积,本发明将p称之为相位矢量(phasor),其角度与fB成正比,范围是[-ππ]。在本发明中,在实际应用中,只要估计出p,无需得到实际的fB,就可得到分离的图像,因此本发明将p也称之为场图。
利用式(2)的p,由式(1)和(2)可解得q:
Figure BDA0001373998240000126
通过利用式(2)、(3),将式(1)写作矩阵的形式:
Figure BDA0001373998240000127
其中
Figure BDA0001373998240000128
是4×1矩阵,
Figure BDA0001373998240000129
Figure BDA00013739982400001210
其中S1、S2分别独立地为在回波时间TE1、TE2下一个像素点的信号强度,
Figure BDA0001373998240000131
Re代表取复数实部,Im代表取复数虚部。
通过式(4),根据变量投影方法(variable projection,VARPRO),使用最小二乘法可解得p:
p=argminerr(p)=argmin||(I-AA+)S||2 (5)。
上式(5)可以改写为以下形式:
Figure BDA0001373998240000132
I是4×4单位矩阵;A+=(ATA)-1AT,上标“+”表示取矩阵伪逆,上标“T”代表矩阵的转置,
Figure BDA0001373998240000133
表示点乘。err(p)描述了场图p和拟合误差之间的关系,表示为场图值拟合误差曲线(phasor-error谱),其定义为不同场图(phasor)值p下对应的模型拟合误差。
公式(6)用于描述感兴趣区域内所有像素点作为一个整体的phasor-error谱。公式(6)达到局部极小值时的p值,数目为1个或2个。
对于描述感兴趣区域(volume of interest,VOI)内所有像素点作为一个整体的phasor-error谱,式(6)依然成立,此时S为:
Figure BDA0001373998240000134
其中
Figure BDA0001373998240000135
Figure BDA0001373998240000136
分别代表一组像素点的信号强度,
Figure BDA0001373998240000137
Figure BDA0001373998240000138
都是行向量,定义如式(4)中表述,行向量中每一个元素代表一个像素点的信号强度。
本发明的基于多尺度局部迭代场图提取算法的磁共振化学位移编码成像方法,实现了对场图值的正确选取和待分离成分的分离。流程图见附图1。
参考图1和图2,本发明的所述方法具体实施如下:
(1)求得二维原始磁共振图像的两组场图候选解,根据以下执行:
1)场图候选解的求解:
根据式(5)描述二维原始磁共振图像中每个像素点的phasor-error谱,根据该phasor-error谱寻找使式(5)达到局部极小值的所有场图值,构成该像素点的场图候选解,个数为1个或2个,以符号p1和p2表示,在场图值数目为1时,p1=p2
(2)在不同的尺度下,将前述两组场图候选解分别分割为不交叠的子块:
1)图像的分块:
本发明将二维原始图像根据预先定义的分块尺度值将图像进行分块,比如选取分块尺度6、8、10和12可得到四组分块的图像,每一组图像中均包含数目不等的子块,每一个子块的大小对应其分块尺度,如当分块尺度是8时,图像中每一个子块的大小均为8×8。
所述分块尺度值可以选自1至16中的一个或多个数值。例如,所述分块尺度值是6、8、10和12。
(3)在每一个尺度下分别采用局部迭代场图提取算法,得到相应的原始图像空间下的场图:
(3a)场图的初始化:
每一个尺度下的每一子块内所有像素点作为整体根据式(6)和(7)描述phasor-error谱,根据该phasor-error谱寻找使式(6)达到局部极小值的场图值,所述场图值为该子块的场图候选解,并设置以下阈值排除一定数目的场图值,使之数目为1个或2个,阈值设定为Amin+α*(Amax-Amin),其中Amax和Amin分别表示phasor-error谱的最大值和最小值,α的选取要排除拟合误差大的场图值,α可以取0.4。该子块的场图候选解以符号pu和pv表示,当场图值数目为1时,pu=pv
以上述子块的场图值pu和pv为初始值,分别对该子块内所有像素点的两个场图候选解(p1,p2)进行选取,以pu为例,选取方法如下:
Figure BDA0001373998240000141
这样对应pu,子块内所有像素点都选取了相应的场图值,组成一个集合Pu,余下的场图值也组成一个集合Puc;同样地,对应pv,存在集合Pv和Pvc
分别计算上述四个集合的矢量和的幅值,记为Mu、Muc、Mv和Mvc。取{Mu,Muc}最大值,记为Umax,取{Mv,Mvc}最大值,记为Vmax,将Umax和Vmax进行比较,舍弃较小数值对应的数据,比如Umax大于Vmax,则舍弃集合Pv和Pvc,保留Pu和Puc,反之,保留Pv和Pvc。为方便表述,假设保留的集合为Pu和Puc,且Mu大于Muc。根据保留的子块场图候选解对应的两组矢量和幅值计算得到一个加权因子C,定义如下:
Figure BDA0001373998240000151
则该子块所有像素点对应的场图值集合为P=Pu*C;如果Muc大于Mu,则P=Puc*C。
对每一个尺度下每一个子块均进行上述场图初始化,将每一个尺度下所有子块的场图值集合P合并在一起,即得到该尺度下的初始场图值集合P0
(3b)场图的更新:
对上述初始场图值集合P0进行低通滤波,且对于下面每一步迭代后得到的场图值集合均进行低通滤波处理。假定迭代后得到的场图值集合为Pn,则滤波后得到Pn s,然后以每一个像素点的场图值pn s为初始值,对图像内所有像素点的两个场图候选解(p1,p2)再次进行选取,选取方法类似于式(8),如下式(10)所示:
Figure BDA0001373998240000152
并记录场图值发生改变的像素点的个数N。
(3c)迭代更新场图:
重复步骤(3b),直至N小于一定的数值或在数次迭代中均固定在某一数值则迭代停止,得到一组场图。N可以为像素点个数的8%-10%。
对应于多组分块尺度,经过步骤4)和5)可得到多组的场图。
(4)将各个尺度下得到的场图合并,保留具有相同场图值的像素点的场图值,其余像素点的场图值置0:
(4a)场图合并:
将上述多组场图进行合并得到一幅场图,具体实现方式是:如果像素点在所有的场图中具有相同的数值,则该保留该像素点的场图值;如果像素点在所有的场图中具有不同的数值,则该像素点的场图值重新置0。
(5)利用局部迭代场图提取算法对置0像素点的场图值进行选取,得到最终场图:
对合并的场图中置0像素点的场图值重复步骤(3b)和(3c),得到最终的场图。
(6)对最终场图进行平滑处理,分离原始磁共振图像,得到准确的各待分离成分的图:
水、脂肪图重建:
利用得到的最终场图,进行平滑处理,然后利用下式计算得到水图和脂肪图:
Figure BDA0001373998240000161
其中ρw、ρf、A+和S如上所定义,根据ρw和ρf得到的水图和脂肪图。
本发明的多尺度的局部迭代场图提取算法(multi-scale RIPE),是在多个尺度下分别进行场图的估计,各尺度之间是相互独立的,提供了一种自校验机制,保证场图估计的正确性。
体组织实验
采用本发明方法在体组织实验进行了测试,其中磁共振扫描操作采用本领域已知的磁共振扫描标准流程,对象为人体颈部,磁共振扫描采用仪器为磁共振成像系统(MR)(联影,上海)。采集序列为TSE,B0=1.5T,成像参数为:矩阵大小=304×320,TE=[0 1.7]ms。实验数据的处理环境是具有Intel E5-2650v2 CPU和64GB RAM的工作站,采用的数据处理软件是MATLAB。
参考图3。图3为颈部的水脂分离结果。(a)颈部冠状位分离结果;(b)颈部矢状位分离结果,图像从左向右依次为分离脂肪图和分离水图。可以看出在体组织实验中,得到的水图和脂肪图中无明显分离错误的组织。
本发明实施例所述的方法不仅可以用于水脂分离成像,也可以用于分离其他化学位移成份的应用中,区别仅仅在于在式(1)中将脂肪相对于水的化学位移值fF替换为待分离的其他化学位移成份对应的化学位移值。
本发明不受上面所显示和描述的实施例的限制,而是在权利要求的范围内可以改变、置换以及等同替代。
参考文献
[1]Dixon WT.Simple proton spectroscopic imaging.Radiology 1984;153(1):189-194.
[2]Xiang QS.Two-point water-fat imaging with partially-opposed phase(POP)acquisition:An asymmetric Dixon method.Magnetic resonance in medicine2006;56(3):572-584.
[3]Eggers H,Brendel B,Duijndam A,Herigault G.Dual-echo Dixon imagingwith flexible choice of echo times.Magnetic resonance in medicine 2011;65(1):96-107.

Claims (20)

1.磁共振化学位移编码成像方法,所述方法包括:
(1)求得二维原始磁共振图像的两组场图候选解;
(2)在不同的分块尺度下,将前述两组场图候选解分别分割为不交叠的子块;
(3)在每一个分块尺度下分别采用局部迭代场图提取算法处理场图候选解子块,得到相应的原始图像空间下的场图;
(4)将各个分块尺度下得到的场图合并,保留具有相同场图值的像素点的场图值,其余像素点的场图值置0;
(5)利用局部迭代场图提取算法对置0像素点的场图值进行选取,得到最终场图;
(6)对最终场图进行平滑处理,分离原始磁共振图像,得到准确的各待分离成分的图;
其中,步骤(2)根据以下执行:
根据预先定义的分块尺度值将所述二维原始磁共振图像进行分块,得到的每一组图像中均包含数目不等的子块,每一个子块的大小对应其分块尺度;
其中,步骤(3)根据以下执行:
(3a)场图的初始化
对所述每一个分块尺度下的每一个子块,求得子块场图候选解,利用每一个子块场图候选解为初始值,对子块内所有像素点的所述两组场图候选解进行选取,选取规则为像素点场图候选解与初始值差值的幅值距离最小,符合所述规则的像素点场图候选解组成一个场图值集合,未被选取的像素点场图候选解组成另一个场图值集合;计算每一个场图候选解对应场图值集合的矢量和;比较所述矢量和的幅值,选取最大幅值对应的场图值集合;
将所述每一个分块尺度下所有子块的场图值集合合并在一起,即得到该分块尺度下的初始场图值集合;
(3b)场图的更新
对上述初始场图值集合进行低通滤波,且对于下面每一步迭代后得到的场图值集合均进行低通滤波处理;然后以每一个像素点的场图值为初始值,对图像内所有像素点的两个场图候选解再次进行选取,选取规则为像素点场图候选解与初始值差值的幅值距离最小;
并记录场图值发生改变的像素点的个数N;
(3c)迭代更新场图
重复步骤(3b),直至N小于一定的数值或在数次迭代中均固定在某一数值则迭代停止,得到一组场图;
对应于多组分块尺度,经过步骤(3a)、(3b)和(3c)得到多组的场图;
其中,步骤(4)根据以下执行:
(4a)场图合并
将所述多组场图进行合并得到一幅场图,具体实现方式是:如果像素点在所有的场图中具有相同的数值,则该保留该像素点的场图值;如果像素点在所有的场图中具有不同的数值,则该像素点的场图值重新置0;
其中,步骤(5)根据以下执行:
对合并的场图中置0像素点的场图值重复步骤(3b)和(3c),得到最终的场图;
其中,所述待分离成分为水和另一成分。
2.根据权利要求1所述的磁共振化学位移编码成像方法,其特征在于,步骤(1)根据以下执行:
描述二维原始磁共振图像中每个像素点的场图值拟合误差曲线,并根据该场图值拟合误差曲线寻找使之达到局部极小值的所有场图值,构成该像素点的两组场图候选解。
3.根据权利要求1所述的磁共振化学位移编码成像方法,其特征在于,所述分块尺度值为整数,选自1至16中的一个或多个数值。
4.根据权利要求3所述的磁共振化学位移编码成像方法,其特征在于,所述分块尺度值为6、8、10和12。
5.根据权利要求1所述的磁共振化学位移编码成像方法,其特征在于,在步骤(3c)中,N为像素点个数的8%-10%。
6.根据权利要求1所述的磁共振化学位移编码成像方法,其特征在于,步骤(6)根据以下执行:
利用得到的最终场图,进行平滑处理,然后计算得到水图和另一成分的图。
7.根据权利要求1-6中任一项所述的磁共振化学位移编码成像方法,其特征在于,所述另一成分为脂肪。
8.一种磁共振化学位移编码成像装置,其特征在于,其包括:
场图候选解求解模块,其从二维原始磁共振图像求得该图像的两组场图候选解;
场图候选解分割模块,其在预先定义的不同的分块尺度下,将在场图候选解求解模块中求得的前述两组场图候选解分别分割为不交叠的子块;
局部迭代场图提取模块,其在场图候选解分割模块定义的每一个分块尺度下分别采用局部迭代场图提取算法处理在场图候选解分割模块中获得的场图候选解子块,得到相应的原始图像空间下的场图;
场图合并模块,其将在局部迭代场图提取模块中各个分块尺度下得到的场图合并,保留具有相同场图值的像素点的场图值,其余像素点的场图值置0;
最终场图获取模块,其利用局部迭代场图提取算法对在场图合并模块中置0像素点的场图值进行选取,得到最终场图;
图像分离模块,其对在最终场图获取模块中得到的最终场图进行平滑处理,分离原始磁共振图像,得到准确的各待分离成分的图;
其中,场图候选解分割模块根据所述预先定义的分块尺度值将二维原始磁共振图像进行分块,得到的每一组图像中均包含数目不等的子块,每一个子块的大小对应其分块尺度;
其中,局部迭代场图提取模块包括:
场图初始化模块,其实现权利要求1中所述的(3a)场图的初始化步骤;
场图更新模块,其实现权利要求1中所述的(3b)场图的更新步骤;以及
迭代更新场图模块,其实现权利要求1中所述的(3c)迭代更新场图步骤;
其中,场图合并模块通过以下实现方式将在局部迭代场图提取模块中获得的多组场图进行合并得到一幅场图:如果像素点在所有的场图中具有相同的数值,则该保留该像素点的场图值;如果像素点在所有的场图中具有不同的数值,则该像素点的场图值重新置0;
其中,最终场图获取模块对在场图合并模块中合并的场图中置0像素点的场图值重复步骤(3b)和(3c),得到最终的场图;
其中,所述待分离成分为水和另一成分。
9.根据权利要求8所述的磁共振化学位移编码成像装置,其特征在于,所述场图候选解求解模块包括:
描述各像素点场图值拟合误差曲线模块,其描述二维原始磁共振图像中每个像素点的场图值拟合误差曲线;以及
获取各像素点场图候选解模块,其根据所述场图值拟合误差曲线寻找使之达到局部极小值的所有场图值。
10.根据权利要求8所述的磁共振化学位移编码成像装置,其特征在于,所述分块尺度值为整数,选自1至16中的一个或多个数值。
11.根据权利要求10所述的磁共振化学位移编码成像装置,其特征在于,所述分块尺度值为6、8、10和12。
12.根据权利要求11所述的磁共振化学位移编码成像装置,其特征在于,所述图像分离模块对在最终场图获取模块中得到的最终场图进行平滑处理,然后计算得到水图和另一成分的图。
13.根据权利要求8-12中任一项所述的磁共振化学位移编码成像装置,其特征在于,所述另一成分为脂肪。
14.一种计算机可读介质,该计算机可读介质具有存储在其中的程序,该程序是计算机可执行的以使计算机执行包括以下步骤的处理:
(1)求得二维原始磁共振图像的两组场图候选解;
(2)在不同的分块尺度下,将前述两组场图候选解分别分割为不交叠的子块;
(3)在每一个分块尺度下分别采用局部迭代场图提取算法处理场图候选解子块,得到相应的原始图像空间下的场图;
(4)将各个分块尺度下得到的场图合并,保留具有相同场图值的像素点的场图值,其余像素点的场图值置0;
(5)利用局部迭代场图提取算法对置0像素点的场图值进行选取,得到最终场图;
(6)对最终场图进行平滑处理,分离原始磁共振图像,得到准确的各待分离成分的图;
其中,步骤(2)根据以下执行:
根据预先定义的分块尺度值将所述二维原始磁共振图像进行分块,得到的每一组图像中均包含数目不等的子块,每一个子块的大小对应其分块尺度;
其中,步骤(3)根据以下执行:
(3a)场图的初始化
对所述每一个分块尺度下的每一个子块,求得子块场图候选解,利用每一个子块场图候选解为初始值,对子块内所有像素点的所述两组场图候选解进行选取,选取规则为像素点场图候选解与初始值差值的幅值距离最小,符合所述规则的像素点场图候选解组成一个场图值集合,未被选取的像素点场图候选解组成另一个场图值集合;计算每一个场图候选解对应场图值集合的矢量和;比较所述矢量和的幅值,选取最大幅值对应的场图值集合;
将所述每一个分块尺度下所有子块的场图值集合合并在一起,即得到该分块尺度下的初始场图值集合;
(3b)场图的更新
对上述初始场图值集合进行低通滤波,且对于下面每一步迭代后得到的场图值集合均进行低通滤波处理;然后以每一个像素点的场图值为初始值,对图像内所有像素点的两组场图候选解再次进行选取,选取规则为像素点场图候选解与初始值差值的幅值距离最小;
并记录场图值发生改变的像素点的个数N;
(3c)迭代更新场图
重复步骤(3b),直至N小于一定的数值或在数次迭代中均固定在某一数值则迭代停止,得到一组场图;
对应于多组分块尺度,经过步骤(3a)、(3b)和(3c)得到多组的场图;
其中,步骤(4)根据以下执行:
(4a)场图合并
将所述多组场图进行合并得到一幅场图,具体实现方式是:如果像素点在所有的场图中具有相同的数值,则该保留该像素点的场图值;如果像素点在所有的场图中具有不同的数值,则该像素点的场图值重新置0;
其中,步骤(5)根据以下执行:
对合并的场图中置0像素点的场图值重复步骤(3b)和(3c),得到最终的场图;
其中,所述待分离成分为水和另一成分。
15.根据权利要求14所述的计算机可读介质,其特征在于,步骤(1)根据以下执行:
描述二维原始磁共振图像中每个像素点的场图值拟合误差曲线,并根据该场图值拟合误差曲线寻找使之达到局部极小值的所有场图值,构成该像素点的两组场图候选解。
16.根据权利要求14所述的计算机可读介质,其特征在于,所述分块尺度值为整数,选自1至16中的一个或多个数值。
17.根据权利要求16所述的计算机可读介质,其特征在于,所述分块尺度值为6、8、10和12。
18.根据权利要求14所述的计算机可读介质,其特征在于,在步骤(3c)中,N为像素点个数的8%-10%。
19.根据权利要求14所述的计算机可读介质,其特征在于,步骤(6)根据以下执行:
利用得到的最终场图,进行平滑处理,然后计算得到水图和另一成分的图。
20.根据权利要求14-19中任一项所述的计算机可读介质,其特征在于,所述另一成分为脂肪。
CN201710674824.0A 2017-08-09 2017-08-09 磁共振化学位移编码成像方法与装置 Active CN109389651B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201710674824.0A CN109389651B (zh) 2017-08-09 2017-08-09 磁共振化学位移编码成像方法与装置

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201710674824.0A CN109389651B (zh) 2017-08-09 2017-08-09 磁共振化学位移编码成像方法与装置

Publications (2)

Publication Number Publication Date
CN109389651A CN109389651A (zh) 2019-02-26
CN109389651B true CN109389651B (zh) 2022-12-13

Family

ID=65414744

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201710674824.0A Active CN109389651B (zh) 2017-08-09 2017-08-09 磁共振化学位移编码成像方法与装置

Country Status (1)

Country Link
CN (1) CN109389651B (zh)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN118067158A (zh) * 2022-11-23 2024-05-24 中国科学院深圳先进技术研究院 基于转换区域和局部场图迭代的化学位移编码成像方法
CN118091514A (zh) * 2022-11-28 2024-05-28 中国科学院深圳先进技术研究院 基于相位解卷绕的化学位移编码成像方法、装置及设备
CN117058265A (zh) * 2023-10-10 2023-11-14 安徽福晴医疗装备有限公司 基于区域迭代向量提取和新的加权图的水脂分离方法

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104382597A (zh) * 2014-11-11 2015-03-04 奥泰医疗系统有限责任公司 一种磁共振成像中的Dixon水脂分离及辨析方法及系统
CN105809661A (zh) * 2014-12-30 2016-07-27 中国科学院深圳先进技术研究院 磁共振成像的图像水脂分离方法和系统
WO2017067130A1 (zh) * 2015-10-21 2017-04-27 华中科技大学 一种气动光学热辐射噪声校正方法与系统
CN106618571A (zh) * 2016-11-16 2017-05-10 深圳先进技术研究院 一种磁共振成像方法和系统

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104382597A (zh) * 2014-11-11 2015-03-04 奥泰医疗系统有限责任公司 一种磁共振成像中的Dixon水脂分离及辨析方法及系统
CN105809661A (zh) * 2014-12-30 2016-07-27 中国科学院深圳先进技术研究院 磁共振成像的图像水脂分离方法和系统
WO2017067130A1 (zh) * 2015-10-21 2017-04-27 华中科技大学 一种气动光学热辐射噪声校正方法与系统
CN106618571A (zh) * 2016-11-16 2017-05-10 深圳先进技术研究院 一种磁共振成像方法和系统

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
Fat-water separation using a region-growing algorithm with self-feeding phasor estimation;chuanli cheng,et al;《Magnetic Resonance in Medicine》;20160614;全文 *
two-point water-fat imaging with partially-opposed phase acquisition:an asymmetric Dixon method;Qing-San Xiang;《Magnetic Resonance in Medicine》;20061231;全文 *

Also Published As

Publication number Publication date
CN109389651A (zh) 2019-02-26

Similar Documents

Publication Publication Date Title
Novikov et al. Rotationally-invariant mapping of scalar and orientational metrics of neuronal microstructure with diffusion MRI
Scherrer et al. Super-resolution reconstruction to increase the spatial resolution of diffusion weighted images from orthogonal anisotropic acquisitions
US9448289B2 (en) Background field removal method for MRI using projection onto dipole fields
CN108375746B (zh) 一种相位反卷折方法与设备
Greve et al. Accurate and robust brain image alignment using boundary-based registration
Rensonnet et al. Towards microstructure fingerprinting: estimation of tissue properties from a dictionary of Monte Carlo diffusion MRI simulations
US9797974B2 (en) Nonrigid motion correction in 3D using autofocusing with localized linear translations
Malavé et al. Reconstruction of undersampled 3D non‐Cartesian image‐based navigators for coronary MRA using an unrolled deep learning model
JP5719968B2 (ja) Mrデータを収集する方法及び装置
CN105556326B (zh) 具有dixon类型的水/脂肪分离的mr成像
CN109389651B (zh) 磁共振化学位移编码成像方法与装置
CN111505553B (zh) 磁共振成像系统和方法
Auría et al. Structured sparsity for spatially coherent fibre orientation estimation in diffusion MRI
Bhushan et al. Improved B0‐distortion correction in diffusion MRI using interlaced q‐space sampling and constrained reconstruction
US11948676B2 (en) Qualitative and quantitative MRI using deep learning
Tabelow et al. Local estimation of the noise level in MRI using structural adaptation
CN113924503B (zh) 时域磁共振的参数图确定
Ma et al. Shearlet-based compressed sensing for fast 3D cardiac MR imaging using iterative reweighting
CN104504657A (zh) 磁共振弥散张量去噪方法和装置
US20170299682A1 (en) Field-mapping and artifact correction in multispectral imaging
Çetingül et al. Group action induced averaging for HARDI processing
US10310044B2 (en) Method of characterizing molecular diffusion within a body from a set of diffusion-weighted magnetic resonance signals and apparatus for carrying out such a method
US11047944B2 (en) Multi-point magnetic resonance imaging
US20160296126A1 (en) Off-resonance correction for pseudo-continuous arterial spin labeling
Close et al. Fourier Tract Sampling (FouTS): A framework for improved inference of white matter tracts from diffusion MRI by explicitly modelling tract volume

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant