CN107072592B - 磁共振成像装置以及定量性磁化率匹配方法 - Google Patents

磁共振成像装置以及定量性磁化率匹配方法 Download PDF

Info

Publication number
CN107072592B
CN107072592B CN201580060823.6A CN201580060823A CN107072592B CN 107072592 B CN107072592 B CN 107072592B CN 201580060823 A CN201580060823 A CN 201580060823A CN 107072592 B CN107072592 B CN 107072592B
Authority
CN
China
Prior art keywords
susceptibility distribution
distribution
magnetic
magnetic susceptibility
magnetic field
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
CN201580060823.6A
Other languages
English (en)
Other versions
CN107072592A (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.)
Fujifilm Healthcare Corp
Original Assignee
Hitachi Ltd
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 Hitachi Ltd filed Critical Hitachi Ltd
Publication of CN107072592A publication Critical patent/CN107072592A/zh
Application granted granted Critical
Publication of CN107072592B publication Critical patent/CN107072592B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01R—MEASURING ELECTRIC VARIABLES; MEASURING MAGNETIC VARIABLES
    • G01R33/00—Arrangements or instruments for measuring magnetic variables
    • G01R33/20—Arrangements or instruments for measuring magnetic variables involving magnetic resonance
    • G01R33/44—Arrangements or instruments for measuring magnetic variables involving magnetic resonance using nuclear magnetic resonance [NMR]
    • G01R33/443—Assessment of an electric or a magnetic field, e.g. spatial mapping, determination of a B0 drift or dosimetry
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/05—Detecting, measuring or recording for diagnosis by means of electric currents or magnetic fields; Measuring using microwaves or radio waves 
    • A61B5/055—Detecting, measuring or recording for diagnosis by means of electric currents or magnetic fields; Measuring using microwaves or radio waves  involving electronic [EMR] or nuclear [NMR] magnetic resonance, e.g. magnetic resonance imaging
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01R—MEASURING ELECTRIC VARIABLES; MEASURING MAGNETIC VARIABLES
    • G01R33/00—Arrangements or instruments for measuring magnetic variables
    • G01R33/20—Arrangements or instruments for measuring magnetic variables involving magnetic resonance
    • G01R33/44—Arrangements or instruments for measuring magnetic variables involving magnetic resonance using nuclear magnetic resonance [NMR]
    • G01R33/48—NMR imaging systems
    • G01R33/54—Signal processing systems, e.g. using pulse sequences ; Generation or control of pulse sequences; Operator console
    • G01R33/56—Image enhancement or correction, e.g. subtraction or averaging techniques, e.g. improvement of signal-to-noise ratio and resolution
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01R—MEASURING ELECTRIC VARIABLES; MEASURING MAGNETIC VARIABLES
    • G01R33/00—Arrangements or instruments for measuring magnetic variables
    • G01R33/20—Arrangements or instruments for measuring magnetic variables involving magnetic resonance
    • G01R33/44—Arrangements or instruments for measuring magnetic variables involving magnetic resonance using nuclear magnetic resonance [NMR]
    • G01R33/48—NMR imaging systems
    • G01R33/54—Signal processing systems, e.g. using pulse sequences ; Generation or control of pulse sequences; Operator console
    • G01R33/56—Image enhancement or correction, e.g. subtraction or averaging techniques, e.g. improvement of signal-to-noise ratio and resolution
    • G01R33/5602—Image enhancement or correction, e.g. subtraction or averaging techniques, e.g. improvement of signal-to-noise ratio and resolution by filtering or weighting based on different relaxation times within the sample, e.g. T1 weighting using an inversion pulse
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01R—MEASURING ELECTRIC VARIABLES; MEASURING MAGNETIC VARIABLES
    • G01R33/00—Arrangements or instruments for measuring magnetic variables
    • G01R33/20—Arrangements or instruments for measuring magnetic variables involving magnetic resonance
    • G01R33/44—Arrangements or instruments for measuring magnetic variables involving magnetic resonance using nuclear magnetic resonance [NMR]
    • G01R33/48—NMR imaging systems
    • G01R33/54—Signal processing systems, e.g. using pulse sequences ; Generation or control of pulse sequences; Operator console
    • G01R33/56—Image enhancement or correction, e.g. subtraction or averaging techniques, e.g. improvement of signal-to-noise ratio and resolution
    • G01R33/5608—Data processing and visualization specially adapted for MR, e.g. for feature analysis and pattern recognition on the basis of measured MR data, segmentation of measured MR data, edge contour detection on the basis of measured MR data, for enhancing measured MR data in terms of signal-to-noise ratio by means of noise filtering or apodization, for enhancing measured MR data in terms of resolution by means for deblurring, windowing, zero filling, or generation of gray-scaled images, colour-coded images or images displaying vectors instead of pixels
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T11/00—2D [Two Dimensional] image generation
    • G06T11/003—Reconstruction from projections, e.g. tomography
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01R—MEASURING ELECTRIC VARIABLES; MEASURING MAGNETIC VARIABLES
    • G01R33/00—Arrangements or instruments for measuring magnetic variables
    • G01R33/20—Arrangements or instruments for measuring magnetic variables involving magnetic resonance
    • G01R33/24—Arrangements or instruments for measuring magnetic variables involving magnetic resonance for measuring direction or magnitude of magnetic fields or magnetic flux
    • G01R33/243—Spatial mapping of the polarizing magnetic field

Landscapes

  • Physics & Mathematics (AREA)
  • Health & Medical Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • High Energy & Nuclear Physics (AREA)
  • General Physics & Mathematics (AREA)
  • Radiology & Medical Imaging (AREA)
  • General Health & Medical Sciences (AREA)
  • Condensed Matter Physics & Semiconductors (AREA)
  • Signal Processing (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Artificial Intelligence (AREA)
  • Biophysics (AREA)
  • Animal Behavior & Ethology (AREA)
  • Surgery (AREA)
  • Public Health (AREA)
  • Veterinary Medicine (AREA)
  • Molecular Biology (AREA)
  • Medical Informatics (AREA)
  • Heart & Thoracic Surgery (AREA)
  • Biomedical Technology (AREA)
  • Pathology (AREA)
  • Theoretical Computer Science (AREA)
  • Magnetic Resonance Imaging Apparatus (AREA)

Abstract

本发明提供一种使用MRI计算出磁化率分布时，通过简单处理减少伪影和噪声，提高计算出的磁化率分布的精度的技术。在该技术中，首先，从通过MRI取得的复图像经由相位图像获得磁场分布。接着，在基于磁场和磁化率的关系式的制约条件下，反复对根据磁场分布计算出的磁化率分布进行平滑化处理。具体地，进行通过反复运算使由预先决定的初始磁化率分布和磁场分布定义的误差函数最小化的最小化处理来计算出上述磁化率分布。此时，每算出更新磁化率分布时，对该更新磁化率分布进行平滑化，该更新磁化率分布是在最小化处理中根据误差函数计算出的。

Description

磁共振成像装置以及定量性磁化率匹配方法
技术领域
本发明涉及一种磁共振成像(MRI：Magnetic Resonance Imaging)技术。尤其，涉及一种使用所取得的图像计算生物体内的磁化率分布的图像处理技术。
背景技术
磁共振成像装置是利用设置于静磁场内的氢原子核(质子)与特定频率的高频磁场产生共振的核磁共振现象的非侵袭性医用图像诊断装置。核磁共振信号根据质子密度、衰减时间等各种物性值而变化，因此可以通过MRI得到的图像描绘出生物体组织的构造或组成、细胞性状等各种生物体信息。
近年来，作为能够通过MRI测定的物性值，关注生物体组织间的磁化率差。磁化率是表示静磁场中的物质被磁化的程度的物性值。在生物体内，存在静脉血中的脱氧血红蛋白或铁蛋白质等顺磁性体、以及成为占据生物体组织的大部分的水或钙化的基础的钙等反磁性体。对生物体组织间的磁化率的差即磁化率差定量地进行图像化，从而能够应用于脑缺血性疾病的诊断或癌症放射治疗效果预测，神经变性疾病的诊断。
使用MRI对生物体组织间的磁化率差进行图像化的方法被称为定量磁化率匹配(QSM：Quantitative Susceptibility Mapping)。QSM方法是根据测量出的MR图像的相位信息计算出因生物体组织间的磁化率差产生的空间性的磁场变化，并使用磁场和磁化率的关系式计算出磁化率分布的方法。
然而，磁场分布是对磁化率分布进行空间性卷积积分而得的分布，因此根据磁场分布求出磁化率分布的情况下，成为卷积积分的逆问题，无法得到唯一的解。
一般，通过最小二乘法，根据磁场分布求出磁化率分布。此时，导入误差函数，将这些最小化的值设为解。作为有代表性的方法，例如有通过磁场分布和点偶极子磁场的k空间上的运算计算出磁化率分布的TKD(Truncated-based K-space Division，基于阶段的k空间的划分)(例如，参照非专利文献1)、将通过TKD方法计算出的磁化率分布和通过阈值处理提取出的微细结构而得的磁化率分布通过反复运算进行合成的迭代SWIM(SusceptibilityWeighted Imaging and Mapping，磁敏感加权成像以及匹配)方法(例如，参照非专利文献2、专利文献1)、使用附有正规化的最小二乘法的MEDI方法(Morphology enabled dipoleinversion，形态启用偶极反演)方法(例如，参照非专利文献3、专利文献2)
现有技术文献
专利文献
专利文献1：美国专利申请公开第2011/0262017号说明书
专利文献2：美国专利申请公开第2011/0044524号说明书
非专利文献
非专利文献1：Shmueli K等，“Magnetic Susceptibility Mapping of BrainTissue In Vivo Using MRI Phase Data”、Magnetic Resonance in Medicine、2009年、62卷、1510-1522页
非专利文献2：Tang J等，“Improving Susceptibility Mapping Using aThreshold-Based K-Space/Image Domain Iterative Reconstruction Approach”，Magnetic Resonance in Medicine，2013年，69卷，1396-1407页
非专利文献3：Liu T等，“Morphology enabled dipole inversion(MEDI)from asingle-angle acquisition:comparison with COSMOS in human brain imaging”，Magnetic Resonance in Medicine,2011年，66卷，777-783页
发明内容
发明要解决的课题
TKD方法将不能进行逆运算的区域取舍为0或任意值。因此，TKD方法虽然具有简单且计算时间短的优点，但因对魔角区域近旁进行取舍，因此存在产生伪影，或计算磁化率的准确度下降的问题。
迭代SWIM方法对通过TKD方法得到的磁化率分布进行可维持微细结构的掩码处理，反复进行更新磁化率分布的处理。虽然能够计算出维持微细结构的磁化率分布，但存在不能完全去除伪影的问题。此外，为了提取微细结构，需要进行各种阈值处理，因此存在若不对每个被检体设定恰当的阈值，则无法恰当地提取微细结构的课题。
在MEDI方法中，以磁场和磁化率为制约条件，通过反复运算使被称为正规化项的给出平滑效果的项最小化，来计算出磁化率分布。这样的带有正规化的最小二乘法是对与正规化参数的值对应的大小的信号值进行取舍的处理，因此需要根据混入到磁场分布中的噪声恰当地设定正规化项。然而，通常，无法事先得到混入到磁场分布中的噪声的信息，因此根据使正规化参数变化多次而计算出的残差的变化曲线计算出恰当的正规化参数。因此，在MEDI方法中，为了设定恰当的正规化参数，需要庞大的计算时间。
本发明是鉴于上述问题提出的，其目的是提供一种使用MRI计算出磁化率分布时，通过简单处理减少伪影和噪声，提高计算出的磁化率分布的准确率和精度的技术。
用于解决课题的手段
在本发明中，从通过MRI取得的复图像经由相位图像获得磁场分布。在基于磁场和磁化率的关系式的制约条件下，反复对根据磁场分布计算出的磁化率分布进行平滑化处理。具体地，进行通过反复运算使由预先决定的初始磁化率分布和磁场分布定义的误差函数最小化的最小化处理来计算出磁化率分布。此时，在最小化处理中，每算出更新磁化率分布时，对该更新磁化率分布进行平滑化，该更新磁化率分布是根据误差函数计算出的。
发明效果
使用MRI计算出磁化率分布时，通过简单处理减少伪影和噪声，提高计算出的磁化率分布的精度。
附图说明
图1的(a)是本发明实施方式的垂直磁场方式的磁共振成像装置的外观图，(b)是本发明实施方式的水平磁场方式的磁共振成像装置的外观图，(c)是本发明实施方式的提高了开放感的磁共振成像装置的外观图。
图2是表示本发明实施方式的MRI装置的概要结构的框图。
图3是本发明实施方式的计算机的功能框图。
图4是本发明实施方式的QSM处理的流程图。
图5是用于说明RSSG(RF-spoiled-Steady-state Acquisition with RewoundGradient-Echo，射频扰相稳态采集反转梯度回波)序列的脉冲序列的说明图。
图6是本发明实施方式的磁场分布计算处理的流程图。
图7是本发明实施方式的磁化率分布计算处理的流程图。
图8是用于说明本发明实施方式的磁化率分布计算处理的流程的说明图。
图9是本发明实施方式的变形例的磁化率分布计算处理的流程图。
图10是用于说明本发明实施方式的变形例的磁化率分布计算处理的流程的说明图。
具体实施方式
对应用本发明的实施方式的一例进行说明。以下，在用于说明实施方式的所有附图中，只要没有特别限定，则对具有相同功能的部件赋予相同的符号，省略其重复的说明。此外，并不通过以下的记述来限定本发明。
[计算磁化率分布的原理]
在本实施方式中，使用MRI装置在预先决定的回波时间拍摄复图像I，根据复图像I计算捕捉到因生物体组织间的磁化率差而产生的磁场变化的磁场分布δ，并根据磁场分布δ计算出磁化率分布χ。通过所谓的QSM方法来计算出磁化率分布。说明实现这些的MRI装置前，说明在QSM方法中根据复图像I计算磁化率分布χ的方法。
在QSM方法中，根据通过梯度回波(Gradient Echo，GrE)方法拍摄到的相位图像，计算出因生物体组织间的磁化率差而产生的局部的磁场变化。
设位置向量为r时，用以下的式(1)表示因组织间的磁化率差而产生的相对的磁场变化(磁场分布)δ(r)。
[公式1]
Figure BDA0001289849760000041
其中，P(r)表示相位图像，γ表示质子的核磁旋转比，B0表示静磁场强度，TE表示回波时间。
此外，通过与静磁场相关的麦克斯韦方程式，使用生物体内的磁化率分布χ(r)，用以下的式(2)表示磁场分布δ(r)。
[公式2]
Figure BDA0001289849760000051
其中，α表示由向量(r’－r)和静磁场方向形成的角度，d(r)表示点偶极子磁场。
如式(2)所示，用磁化率分布χ(r)和点偶极子磁场d(r)的卷积积分来表示磁场分布δ(r)。因此，通过对式(2)的两边进行傅里叶变换，将式(2)变换为以下的式(3)。
Figure BDA0001289849760000052
其中，k＝(kx、ky、kz)表示k空间上的位置向量，Δ(k)、X(k)、D(k)表示磁场分布δ(r)、磁化率分布χ(r)、点偶极子磁场d(r)的傅里叶分量。
如式(3)所示，能够将磁场分布的傅里叶分量Δ(k)除以点偶极子磁场的傅里叶分量D(k)来计算出磁化率分布的傅里叶分量X(k)。然而，在式(3)中，在D(k)＝0的近旁区域，其倒数发散，因此无法直接计算出X(k)。
将该成为D(k)＝0的区域称作魔角(Magic Angle)，成为相对于磁场方向具有大致54.7°的2倍的顶角的反向双圆锥区域。因魔角的存在，根据磁场分布计算出磁化率分布的QSM方法回归到病态逆问题(Ill-conditioned inverse problem)。
作为解决该问题的有代表性的方法，如上所述，有TKD方法、迭代SWIM方法(以下，简称为SWIM方法)、MEDI方法。
在TKD方法中，取点偶极子磁场d(r)的傅里叶分量D(k)的倒数，将其魔角区域近旁舍位至0或任意值。之后，乘以所测量的磁场分布δ(r)的傅里叶分量Δ(k)，在实际空间进行傅里叶逆变换来计算出磁化率分布χ(r)。
在迭代SWIM方法中，对通过TKD方法计算出的磁化率分布，通过阈值处理提取微细结构。接着，计算出对微细结构以外的磁化率进行了掩码处理后的微细结构磁化率分布。在k空间上混合原本的磁化率分布和微细结构磁化率分布，并计算修正磁化率分布。对于修正磁化率分布，通过阈值处理再次提取微细结构。反复进行这样的一系列处理，来计算出磁化率分布。
在MEDI方法中，如以下的式(4)所示，以磁场和磁化率的关系式为制约条件，通过反复运算使称为正则化的项的给出平滑效果的项最小化，来计算出磁化率分布χ(r)。
[公式4]
Figure BDA0001289849760000061
其中，χ表示计算对象的磁化率分布向量，δ表示根据相位图像计算出的磁场分布向量，W表示在对角分量具有基于绝对值图像计算出的加权向量的矩阵，M表示在对角分量具有基于绝对值图像的空间梯度将边缘区域设为0，将其他区域设为1的二值边缘掩码向量的矩阵。此外，FD是表示卷积积分的算子，
Figure BDA0001289849760000062
(纳布拉，Nabla)是表示三维空间梯度算子，∥a∥n是表示任意向量a的Ln范数(norm)。此外，式(4)的第2项λ为取得平滑化效果和磁化率精度的平衡的正规化参数。
如式(4)所示，关于MEDI方法，作为正规化项，使用对计算过程中的磁化率分布的全变差(TV：Total Variation)
Figure BDA0001289849760000063
乘以基于绝对值图像的边缘掩码M而得的向量的L1范数。一般，任意向量a的带有L1范数正规化的最小二乘解成为稀疏(sparse：疏)的解。因此，在MEDI方法中，绝对值图像的变化较少的区域被平滑化，计算出保存有变化较大的边缘部分的磁化率分布。
可以将包含MEDI方法的带有正规化的最小二乘法解释为对与正规化参数的值对应的大小的信号值进行舍位的处理。因此，当正规化参数过大时噪声和微细结构都被去除，并且图像整体的磁化率减少，因此计算磁化率下降。此外，当正规化参数过小时，虽然维持微细结构或图像整体的磁化率，但噪声增加。
因此，在MEDI方法中，需要设定恰当的正规化参数。设定了恰当的正规化参数的情况下，式(4)的第1项的残差的标准偏差等于混有测量出的磁场分布的噪声的标准偏差。然而，通常无法提前得到混入到磁场分布中的噪声的信息，因此根据使正规化参数变化多次后计算出的残差的变化曲线计算出恰当的正规化参数。因此，在MEDI方法中，为了设定恰当的正规化参数，需要大量的计算时间。
在本实施方式中，根据上述公知的各方法，提供一种使用MRI计算出磁化率分布时，通过简单的处理降低伪影(Artifact)和噪声，并且能够计算出高准确度且高精度的磁化率分布的MRI装置以及QSM方法。
[MRI装置的外观]
首先，对本实施方式的磁共振拍摄装置(MRI装置)进行说明。图1的(a)～图1的(c)是本实施方式的MRI装置的外观图。图1的(a)是使用通过电磁线圈生成静磁场的隧道型磁铁的水平磁场方式的MRI装置100。图1的(b)是为了提高开放感将磁铁上下分离的汉堡包型(开放型)垂直磁场方式的MRI装置120。此外，在图1的(c)是使用与图1的(a)相同的隧道型磁铁，使磁铁的纵深变短，且倾斜而提高了开放感的MRI装置130。
在本实施方式中，也可以使用具有这些外观的MRI装置的任一个。
另外，这些为一例，本实施方式的MRI装置并不限定于这些方式。在本实施方式中，无论装置的方式、类型，可以使用各种的各种MRI装置。
以下，在不需要特别区别的情况下，以MRI装置100为代表。此外，在本实施方式中，例如使用设MRI装置100的静磁场方向为z方向，与之垂直的2个方向中与载置测定对称的被检体的床面平行的方向为x方向，另一方向为y方向的坐标系。此外，以下将静磁场简称为磁场。
[MRI装置的功能结构]
图2是本实施方式的MRI装置100的功能结构图。如本图所示，本实施方式的MRI装置100具备：静磁场线圈102，其由在放置有被检体101的空间生成静磁场的、例如静磁场线圈、静磁场生成磁铁等构成；匀场线圈104，其调整静磁场分布；发送用高频线圈105(以下，简称为发送线圈)，其对被检体101的测量区域发送高频磁场；接收用高频线圈106(以下，简称为接收线圈)，其接收从被检体101产生的核磁共振信号；梯度磁场线圈103，其为了向从被检体101产生的核磁共振信号附加位置信息，分别向x方向、y方向、z方向施加梯度磁场；发送机107；接收机108；计算机109；梯度磁场用电源部112；匀场用电源部113；以及顺序控制装置(sequence controller)114。
静磁场线圈102根据图1的(a)、图1的(b)、图1的(c)所示的各MRI装置100、120、130的结构，采用各种方式。梯度磁场线圈103及匀场线圈104分别由梯度磁场用电源部112及匀场用电源部113驱动。另外，在本实施方式中，以发送线圈105和接收线圈106分别使用不同的线圈的情况为例进行说明，但也可以由兼作发送线圈105和接收线圈106的功能的1个线圈构成。发送线圈105照射的高频磁场由发送机107生成。将接收线圈106检测出的核磁共振信号通过接收机108发送给计算机109。
顺序控制装置114控制梯度磁场线圈103的驱动用电源即梯度磁场用电源部112、匀场线圈104的驱动用电源即匀场用电源部113、发送机107以及接收机108的动作，并控制梯度磁场、高频磁场的施加定时以及核磁共振信号的接收定时。控制的时序图被称作脉冲时序图，根据测量被预先设定，存储在后述的计算机109所具备的存储装置等中。
计算机109控制MRI装置100整体的动作，并且对接收的核磁共振信号进行各种运算处理。在本实施方式中，生成至少一个以上的回波时间的复图像或相位图像、磁场分布、磁化率分布等。计算机109是具备CPU、存储器、存储装置等的信息处理装置，显示器110、外部存储装置111、输入装置115被连接到计算机109。
显示器110是将通过运算处理得到的结果等显示给操作员的接口。输入装置115是操作员用于输入在本实施方式中实施的测量或运算处理所需要的条件、参数等的接口。在本实施方式的输入装置115中，用户可以输入要测量的回波数、基准的回波时间、回波间隔等测量参数。外部存储装置111与存储装置保持计算机109执行的各种运算处理所使用的数据、通过运算处理得到的数据、所输入的条件、参数等。
[计算机的功能框图]
如上所述，在本实施方式中，以基于磁场和磁化率的关系式的制约条件，反复运算平滑化处理，从而能够简单且高精度地计算出磁化率分布。接着，对实现这些的本实施方式的计算机109的功能结构进行说明。
图3是表示本实施方式的计算机109的功能框图。本实施方式的计算机109具备：测量控制部310，其测量与高频磁场脉冲的照射对应地从被检体产生的核磁共振信号(回波信号)作为复信号；图像重构部320，其根据由测量控制部310测量出的复信号重构像素值为复数的复图像；磁场分布计算部330，其根据由图像重构部320重构的复图像生成相位图像，并根据该相位图像计算出磁场分布；以及磁化率分布计算部340，其根据由磁场分布计算部330计算出的磁场分布计算磁化率分布。
另外，计算机109具备CPU、存储器和存储装置。并且，CPU将存储装置所保持的程序加载到存储器中并执行，由此实现计算机109实现的各种功能。将各功能处理所使用的各种数据、处理中生成的各种数据存储到存储装置或外部存储装置111中。
此外，计算机109实现的各种功能中至少一个功能也可以由独立于MRI装置100的信息处理装置即可与MRI装置100收发数据的信息处理装置来实现。并且，全部或一部分功能也可以由ASIC(Application Specific Integrated Circuit，专用集成电路)、FPGA(Field-Programmable Gate Array，现场可编程门阵列)等硬件来实现。
[QSM处理的流程]
对基于本实施方式的计算机109的上述各功能的QSM方法的处理(QSM处理)的流程进行说明。图4是用于说明本实施方式的QSM处理的流程的流程图。
首先，测量控制部310按照预先决定的脉冲序列控制顺序控制装置114，测量预先决定的回波时间的回波信号(步骤S1001)。
之后，图像重构部320将得到的回波信号配置于k空间上并进行傅里叶变换，来对复图像I进行重构(步骤S1002)。
磁场分布计算部330进行根据复图像I计算磁场分布δ的磁场分布计算处理(步骤S1003)。
之后，磁化率分布计算部340根据计算出的磁场分布δ，进行计算磁化率分布χ的磁化率分布计算处理(步骤S1004)，并将算出的磁化率分布χ显示于显示器110(步骤S1005)。
另外，在将磁化率分布χ显示于显示器110时，根据需要，除在步骤S1003中计算出的磁场分布δ等外，还可以将在QSM处理的过程中计算出的图像显示于显示器110。
以下，对本实施方式的QSM处理的各处理的细节进行说明。
[测量：S1001]
测量控制部310按照根据用户经由输入装置115输入的参数设定的脉冲序列，使顺序控制装置114动作，进行得到预先决定的回波时间(TE)的核磁共振信号(回波信号)的测量。顺序控制装置114按照来自测量控制部310的指示，控制MRI装置100的各部来进行测量。在本实施方式中，得到至少一个回波时间的回波信号。
在此，对测量控制部310在测量中所使用的脉冲序列的例子进行说明。在本实施方式中，例如使用GrE(Gradient Echo，梯度回波)系的脉冲序列。该GrE系的脉冲序列对因生物体组织间的磁化率差而生成的局部的磁场变化敏锐。
图5中，作为GrE系的脉冲序列的一例，示出了RSSG(RF-spoiled-Steady-stateAcquisition with Rewound Gradient-Echo，射频扰相稳态采集反转梯度回波)序列510的例子。在该图中，RF、Gs、Gp、Gr分别表示高频磁场、切片梯度磁场、相位编码梯度磁场、读出梯度磁场。
在RSSG序列510中，施加切片梯度磁场脉冲501，并且照射高频磁场(RF)脉冲502，激发被检体101内预定切片的磁化。接着，向磁化相位施加用于附加切片方向和相位编码方向的位置信息的切片编码梯度磁场脉冲503和相位编码梯度磁场脉冲504。
施加使像素内的核磁化的相位分散的失相(dephase)用读出梯度磁场脉冲505后，边施加用于附加读出方向的位置信息的读出梯度磁场脉冲506，边测量一个核磁共振信号(回波)507。然后，最后施加使被切片编码梯度磁场脉冲503和相位编码梯度磁场脉冲504失相的核磁化的相位收敛的重相位用切片编码梯度磁场脉冲508和相位编码梯度磁场脉冲509。
测量控制部310边使切片编码梯度磁场脉冲503、508(切片编码数ks)和相位编码梯度磁场脉冲504、509(相位编码数kp)的强度、RF脉冲502的相位变化，边在反复时间TR反复执行以上的过程，为了得到1幅图像测量必要的回波。另外，此时使RF脉冲502的相位例如每次增加117度。此外，在图5中，连字符以下的数字表示是反复的第几次。
另外，测量的各回波被配置于以kr、kp、ks为坐标轴的三维k空间上。此时，一个回波在k空间上占与kr轴平行的一条线。关于通过该RSSG序列510得到的绝对值图像，在将TE设定得短时成为T1(纵衰减时间)强调图像，在将TE设定得长时成为反映像素内的相位分散的T2*强调图像。
另外，在此示例了与k空间的坐标轴平行地取得数据的笛卡儿(Cartesian)拍摄的一例即RSSG序列510，但在本实施方式中，可以使用任意序列来取得任意k空间区域的数据。
例如，也可以使用通过一次的TR取得多个回波，并将这些配置于不同的k空间的序列。此外，也可以是多次进行测量，通过各测量取得不同的k空间区域的数据，对这些数据进行整合，得到1幅图像的序列。此外，也可以使用在k空间旋转状地取得数据的径向扫描等非笛卡儿拍摄。
并且，也可以使用根据多个TE取得的不同的复图像、相位图像生成一幅复图像、相位图像的序列。例如，使用可实现根据通过多个TE取得的信号得到图像的Dixon方法的脉冲序列的情况下，在分离水和脂肪的过程得到磁场分布δ。因此，在该情况下不需要后述的生成相位图像P的处理。
[图像重构：S1002]
接着，对上述步骤S1002的图像重构处理进行说明。图像重构部320将在步骤S1001测量出的回波时间的回波信号配置于k空间，并进行傅里叶变换。由此，图像重构部320计算出各像素值为复数的复图像I。
[磁场分布计算处理：S1003]
接着，对上述步骤S1003的磁场分布计算部330进行的磁场分布计算处理进行说明。本实施方式的磁场分布计算部330根据图像重构部320重构的复图像I计算出磁场分布δ。计算出的磁场分布δ是取了因生物体组织间的磁化率差产生的局部的磁场变化的磁场分布。
磁场分布计算部330，首先根据复图像I计算出相位图像P，并对计算出的相位图像P进行用于去除背景相位的各种相位图像处理。然后，根据相位图像处理后的相位图像P取得磁场分布δ。
在相位图像处理中，去除源于磁化率以外的相位分量，提取源于磁化率的相位分量。具体地，在计算出的相位图像P中，进行去除超过－π至π的范围后折返的相位的相位展开(unwrap)处理、和去除因生物体形状产生的整体(global)的相位变化的整体相位变化去除处理。
为了实现这些，本实施方式的磁场分布计算部330，如图3所示，具备：根据复图像I计算出相位图像P的相位图像计算部331、进行相位展开处理的相位展开处理部332、进行整体相位变化去除处理的整体相位变化去除部333、和将相位处理后的相位图像P变换为磁场分布δ的磁场分布变换部334。
相位图像计算部331求出复图像I的各像素的像素值的偏角，并将该偏角设为各个像素值来计算出相位图像P。即，相位图像P是以复图像I的各像素的复数的相位分量为像素值的图像。此处，使用复图像I的像素值I(r)，通过式(5)来计算出像素r中的相位图像P的像素值P(r)。
P(r)＝arg{I(r)}···(5)
相位展开处理部332对相位图像P进行用于去除超过－π至π的范围后折返的相位的相位展开处理。在相位图像P中的部分区域，超过－π至π的范围的相位值在－π至π的范围内折返。因此，为了在拍摄部位(例如，头部等)的全部区域求出准确的相位，需要对这些折返进行修正。在本实施方式中，例如使用区域扩大法等，对在－π至π的范围折返的相位值进行修正。
整体相位变化去除部333在相位展开处理后的相位图像P中去除因生物体形状而产生的整体的相位变化。
相位图像P的各像素值是取决于摄像部位的形状等而产生的静磁场不均匀而引起的整体的相位变化、组织间的磁化率变化引起的局部的相位变化的和。整体的相位变化和局部的相位变化分别对应于相位图像P的k空间上的低频分量和高频分量。通过本处理，去除整体的相位变化，并提取组织间的磁化率变化引起的局部的相位变化。
本实施方式的整体相位变化去除部333，首先对拍摄的三维图像(原图像)，针对每个二维图像分别实施低通滤波处理，并计算出低分辨率图像，作为整体相位变化去除处理。然后，将原图像用低分辨率图像进行复数除法运算。由此，从原图像去除低分辨率图像所包含的整体的相位变化。
另外，去除整体的相位变化的方法包括各种方法。例如，有将三维图像以低阶多项式进行拟合来提取整体的相位变化，并将这些从原图像中减去的方法等。在整体相位变化去除中，也可以使用这样的其他的方法。
另外，相位图像处理的相位展开处理和整体相位变化去除处理仅为一例，并不限定于此。此外，也可以省略这2个处理中的任一个处理。此外，各处理的处理顺序是任意的。此外，也可以不进行这2个处理。
磁场分布变换部334使用式(1)，将相位图像处理后的相位图像P变换成磁场分布δ。
[磁场分布计算处理的流程]
对基于本实施方式的磁场分布计算部330的各功能的磁场分布计算处理的流程的一例进行说明。图6是本发明实施方式的磁场分布计算处理的处理流程。
相位图像计算部331根据复图像I计算出相位图像P(步骤S1101)。
然后，磁场分布计算部330对相位图像P进行各种相位图像处理，来得到相位处理后的相位图像。在本实施方式中，例如相位展开处理部332对计算出的相位图像P实施相位展开处理(步骤S1102)，整体相位变化去除部333对相位展开处理后的相位图像P进行整体相位变化去除处理(步骤S1103)。如上所述，无论这些相位图像处理的处理顺序。此外，相位图像处理可以进行任一个，也可以完全不进行。并且，并不限定于这些。
最终，磁场分布变换部334将位相处理后的相位图像P变换为磁场分布(步骤S1104)，结束处理。
[磁化率分布计算处理：S1004]
接着，对上述步骤S1004的磁化率分布计算部340进行的磁化率分布计算处理进行说明。本实施方式的磁化率分布计算部340使用由磁场分布计算部330计算出的磁场分布δ，来计算出磁化率分布χ。
磁化率分布计算部340利用通过反复运算使误差函数最小化的最小化处理来计算出磁化率分布χ，其中，误差函数是由预先决定的初始磁化率分布χ0和磁场分布δ定义的函数。此时，在最小化处理中，对于根据误差函数计算出的更新磁化率分布χUPD在每次算出该更新磁化率分布时进行平滑化。另外，更新磁化率分布χUPD是将初始磁化率分布χ0根据误差函数更新后的磁化率分布。
为了实现这些，本实施方式的磁化率分布计算部340具备：误差函数设定部341，其设定误差函数；更新磁化率分布计算部342，其使用所设定的误差函数，根据磁场分布δ和初始磁化率分布χ0计算更新磁化率分布χUPD；平滑化处理部343，其对更新磁化率分布χUPD进行平滑化，并生成合成磁化率分布χADD；以及判定部346，其进行反复运算的结束判定，并且将初始磁化率分布χ0更新为计算出的合成磁化率分布χADD。
此外，平滑化处理部343具备：滤波处理部344，其对更新磁化率分布χUPD应用预先决定的平滑化滤波器，来生成滤波处理后的磁化率分布(处理后磁化率分布)χSMT；以及合成部345，其对更新磁化率分布χUPD和处理后磁化率分布χSMT进行合成，生成合成磁化率分布χADD。
[误差函数设定部]
如式(2)所示，磁场分布δ是在空间上对磁化率分布χ进行卷积积分而得的分布。以相位图像P内的全部像素为对象，因此通过向量和矩阵表现式(2)时，表现为以下的式(6)。
δ＝FDχ···(6)
其中，δ为具有全部像素数N的大小的磁场分布的列向量，χ为计算出的磁化率分布的列向量，此外，FD为具有N×N大小，相当于对磁化率分布χ的卷积运算的矩阵。
本实施方式的更新磁化率分布计算部342根据式(6)，通过最小二乘法来求出磁化率分布χ。此时，导入误差函数，并使其最小化来求出解。
本实施方式的误差函数设定部341设定更新磁化率分布计算部342计算更新磁化率分布时所使用的误差函数。例如，通过以下的式(7)来表示所设定的误差函数e(χ0)。
[公式7]
Figure BDA0001289849760000141
其中，W表示对角分量具有根据复图像I的相位值分量计算出的加权向量的矩阵(以下，称为加权W),FD是表示卷积积分的算子，∥a∥2 2表示任意向量a的L2范数(欧几里得范数)的平方。
即，误差函数对由磁场分布δ和初始磁化率分布χ0进行卷积运算而得的结果决定的函数赋予加权W后的函数。如上所述根据复图像的相位值分量(相位图像P)，计算出该加权W。
具体地，误差函数设定部341根据相位图像P计算出各像素的相位偏差，根据所计算出的相位偏差的大小来计算出加权W。例如，相位偏差越大，则将加权W计算得越小。
对相位偏差例如使用标准偏差等统计指标，该标准偏差为使用计算相位偏差的对象像素的边缘区域的多个像素的像素值(相位值)而计算出的偏差。
另外，使用的相位图像P既可以是根据复图像I得到的图像，也可以是对根据复图像I得到的相位图像P实施上述各种相位图像处理而得到的图像。
[更新磁化率计算部]
如上所述，更新磁化率分布计算部342使用误差函数e(χ0)，从初始磁化率分布χ0计算出更新磁化率分布χUPD。在此，更新磁化率分布计算部342计算出误差函数e(χ0)的梯度
Figure BDA0001289849760000151
使用计算出的勾配
Figure BDA0001289849760000152
通过以下的式(8)计算出更新磁化率分布χUPD。
Figure BDA0001289849760000153
另外，在本实施方式中，最初设定的初始磁化率分布χ0可以被设定为任意的分布。例如，可以将全部值为0的分布设为磁化率分布χ0。
[平滑化处理部；滤波处理部]
平滑化处理部343使用平滑化滤波器对更新磁化率分布χUPD进行平滑化，得到处理后磁化率分布χSMT。其中，滤波处理部344使用边缘保存平滑化滤波器作为平滑化滤波器，与更新磁化率分布χUPD中的各像素值的局部的磁化率变动对应地适当地进行线性平滑化处理。
具体地，将更新磁化率分布χUPD中的任意的像素位置设为r时，若将其周边的n×n(n为1以上的整数)的局部区域的平均值设为μ(r)，将方差设为σ(r)2时，位置r处的边缘保存平滑化滤波器应用后的磁化率分布即处理后磁化率分布χSMT的像素值χSMT(r)用以下的式(9)来表示
[公式9]
Figure BDA0001289849760000161
其中，ν2为在磁化率分布χ混入了相当于噪声方差的参数。在本实施方式中，采用更新磁化率分布χUPD的各像素中的局部方差σ(r)2的平均值作为该参数。
根据在本实施方式中应用的边缘保存平滑化滤波器，如式(9)所示，在局部方差小的区域平滑化较强，成为接近局部平均值，因此能够维持磁化率的定量值。此外，在局部方差大的区域平滑化变弱，保存包含微细结构的边缘部分。因此，通过向更新磁化率分布χUPD应用该边缘保存平滑化滤波器，能够降低磁化率分布的噪声，边维持定量值，边自动提取并保存细微静脉或磁化率变化大的组织间的边缘。
另外，在本实施方式中，使用了上述的边缘保存平滑化滤波器，但在平滑化处理中使用的滤波器并不限定于此。例如，也可以使用中值滤波器或双边滤波器等公知的边缘保存平滑化滤波器。
[平滑化处理部；合成部]
合成部345对更新磁化率分布χUPD和处理后磁化率分布χSMT进行合成来计算出合成磁化率分布χADD。为了防止噪声的增加，且维持磁化率分布的微细结构而进行合成。例如，以向噪声的影响大的k空间的魔角近旁插入被边缘保存平滑化后的处理后磁化率分布χSMT的傅里叶分量的方式进行合成。
合成部345将更新磁化率分布χUPD和处理后磁化率分布χSMT变换为k空间数据，分别对变换后的更新磁化率分布χUPD和处理后磁化率分布χSMT的k空间数据加权并进行加法运算，来得到合成磁化率分布χADD。
具体地，合成部345首先将更新磁化率分布χUPD和处理后磁化率分布χSMT通过傅里叶变换变换为k空间数据。然后，对变换后的各k空间数据(傅里叶分量进行加权以便达成上述目的，对加权后的各k空间数据进行加法运算。然后，计算后，通过傅里叶逆变换得到合成磁化率分布χADD。
通过生成权重图像G，并对k空间数据乘以权重图像G来进行k空间的加权。此时，在与根据更新磁化率分布χUPD得到的k空间数据相乘的权重图像GUPD中，例如将向魔角区域近旁的数据赋予的权重设定得比向其他区域的数据赋予的权重小。此外，在与根据处理后磁化率分布χSMT得到的k空间数据相乘的权重图像GSMT中，相反，将向魔角区域近旁的数据赋予的权重设定得比向其他区域的数据赋予的权重大。
这是因为更新磁化率分布χUPD的噪声比较多，因此降低魔角区域近旁的数据的贡献，来防止噪声的增加。此外，处理后磁化率分布χSMT因噪声降低，使魔角区域近旁的数据的贡献增加，来维持微细结构。
作为满足这些条件的权重图像GUPD(k)、GSMT(k)，可以使用点偶极子磁场的傅里叶分量D(k)和调整魔角区域大小的调整参数b，例如可以使用由以下的式(10)定义的权重图像。另外，k表示k空间的数据位置。
[公式10]
Figure BDA0001289849760000171
在上述式(10)中，在更新磁化率分布χUPD中，将傅里叶分量的值小于预定值(b)的区域设为魔角区域。
另外，在调整参数b为0的情况下，GUPD(k)成为1，GSMT(k)成为0。因此，不实施边缘保存平滑化处理，而成为使式(6)的误差函数最小化的简单的最小二乘法。
[判定部]
在通过反复运算使误差函数e(χ0)最小化的最小化处理中，判定部346进行反复运算的结束判定。在判定为继续进行处理，即反复的情况下，将在该时间点得到的最新的合成磁化率分布χADD设定为初始磁化率分布χ0，来更新初始磁化率分布χ0。另外，在判定为结束时，输出该时间点的最新的合成磁化率分布χADD作为磁化率分布。
在本实施方式中，判定部346在计算更新磁化率分布χUPD后，使用该更新磁化率分布χUPD和该时间点的初始磁化率分布χ0，来判定是否结束。判定部346例如使用以下的式(11)所示的评价函数f(χUPD)，在其值小于预先决定的阈值ε的情况下，判定为结束。
[公式11]
Figure BDA0001289849760000181
[磁化率分布计算处理的流程]
接着，对本实施方式的磁化率分布计算部340进行的磁化率分布计算处理的流程进行说明。图7是本实施方式的磁化率分布计算处理的处理流程。此外，图8是用于说明磁化率分布计算处理的说明图。
在本实施方式中，通过反复运算进行使利用预先决定的初始磁化率分布χ0和磁场分布δ定义的误差函数最小化的最小化处理来计算出磁化率分布χ。此时，在最小化处理中，对根据误差函数计算出的更新磁化率分布进行平滑化处理，得到合成磁化率分布，通过该合成磁化率分布更新初始磁化率分布，反复进行处理。具体地，进行以下的处理。
磁化率分布计算部340设定初始磁化率分布χ0410(步骤S1201)。然后，误差函数设定部341利用初始磁化率分布χ0410和磁场分布δ400设定误差函数e(χ0)(步骤S1202)。
然后，更新磁化率分布计算部342计算出误差函数的梯度
Figure BDA0001289849760000182
(步骤S1203)，按照上述式(8)，对初始磁化率分布χ0进行更新，并计算出更新磁化率分布χUPD420(步骤S1204)。
判定部346使用上述式(11)的评价函数f(χUPD)来进行结束判定(步骤S1205)。
在步骤S1205中，判定为继续的情况下，滤波处理部344对更新磁化率分布χUPD420应用预先决定的平滑化滤波器，计算出处理后磁化率分布χSMT430(步骤S1206)。
合成部345对更新磁化率分布χUPD420和处理后磁化率分布χSMT430进行合成来得到合成磁化率分布χADD440(步骤S1207)。对更新磁化率分布χUPD420和处理后磁化率分布χSMT430分别进行傅里叶变换(FT)而得到的k空间数据421、431乘以上述式(10)的权重图像422、432并进行加法运算，对其加权相加的结果进行傅里叶逆变换(IFT)来进行合成。
判定部346将合成磁化率分布χADD440设定为初始磁化率分布χ0410(步骤S1208)，返回到步骤S1202，结束处理。
在此，在步骤S1205中，在判定为结束的情况下，判定部346将在该时间点得到的最新的合成磁化率分布χADD设为磁化率分布计算部340输出的磁化率分布χ(步骤S1209)，结束处理。
另外，在步骤S1205中，在没有算出合成磁化率分布χADD的情况下，即，通过第1次的运算判断为结束的情况下，通过上述方法根据该时间点的更新磁化率分布χUPD计算出合成磁化率分布χADD，并将计算出的合成磁化率分布χADD设为要输出的磁化率分布χ。
＜磁化率分布计算处理流程的变形例＞
另外，判定部346也可以使用合成磁化率分布χADD来判定是否结束。即，计算出合成磁化率分布χADD时，也可以使用该合成磁化率分布χADD和该时间点的初始磁化率分布χ0，来判定是否结束。
此时，作为判定结束的评价函数，可以使用以下的f(χADD)。
[公式12]
Figure BDA0001289849760000191
图9表示此时的磁化率分布计算部340的各部的磁化率分布计算处理流程。此外，图10是用于说明本处理的说明图。
磁化率分布计算部340设定初始磁化率分布χ0410(步骤S1301)。然后，误差函数设定部341利用初始磁化率分布χ0410和磁场分布δ400设定误差函数e(χ0)(步骤S1302)。
然后，更新磁化率分布计算部342计算出误差函数的梯度
Figure BDA0001289849760000192
(步骤S1303)，按照上述式(8)，计算出更新磁化率分布χUPD420(步骤S1304)。
滤波处理部344对更新磁化率分布χUPD420应用预先决定的平滑化滤波器，计算出处理后磁化率分布χSMT430(步骤S1305)。
合成部345对更新磁化率分布χUPD420和处理后磁化率分布χSMT430进行合成来得到合成磁化率分布χADD440(步骤S1306)。对更新磁化率分布χUPD420和处理后磁化率分布χSMT430分别进行傅里叶变换(FT)而得到的k空间数据421、431乘以上述式(10)的权重图像422、432并进行加法运算，对其加权相加的结果进行傅里叶逆变换(IFT)来进行合成。
接着，判定部346使用上述式(12)的评价函数f(χADD)来进行结束判定(步骤S1307)。
在步骤S1307中，在判定为处理继续的情况下，磁化率分布计算部340将合成磁化率分布χADD440设定为初始磁化率分布χ0410(步骤S1308)，返回到步骤S1302，重复处理。
在此，在步骤S1307中，在判定为结束的情况下，将在该时间点得到的最新的合成磁化率分布χADD设为磁化率分布计算部340输出的磁化率分布χ(步骤S1309)，结束处理。
通过合成磁化率分布χADD进行判定时，与通过更新磁化率分布χUPD进行判定时相比，能够得到滤波处理的效果大的图像。另一方面，通过更新磁化率分布χUPD进行判定时，与通过合成磁化率分布χADD进行判定时相比，能够以与实际的磁化率分布接近的数据进行判定，因此能够在更确切的范围内进行判定。
如上所述，本实施方式的MRI装置具备：测量控制部310，其根据高频磁场脉冲的照射将从被检体101发生的核磁共振信号测量为复信号；图像重构部320，其根据上述复信号重构像素值为复数的复图像I；磁场分布计算部330，其根据上述复图像I生成相位图像P，根据该相位图像P计算出磁场分布δ；以及磁化率分布计算部340，其根据上述磁场分布δ计算出磁化率分布χ，上述磁化率分布计算部340进行通过反复运算使利用预先决定的初始磁化率分布χ0和上述磁场分布δ定义的误差函数e(χ0)最小化的最小化处理来计算出上述磁化率分布χ，并且，在该最小化处理中，每次算出更新磁化率分布χUPD时，对根据上述误差函数e(χ0)计算出的该更新磁化率分布χUPD进行平滑化。
上述磁化率分布计算部340具备：更新磁化率分布342，其使用上述误差函数e(χ0)，根据上述初始磁化率分布χ0计算出上述更新磁化率分布χUPD；平滑化处理部343，其对上述更新磁化率分布χUPD进行平滑化来生成合成磁化率分布χADD；以及判定部346，其进行上述反复运算的结束判定，并且，将上述合成磁化率分布χADD设定为上述初始磁化率分布χ0，上述平滑化处理部343具备：滤波处理部344，其对上述更新磁化率分布χUPD应用预先决定的平滑化滤波器，得到处理后磁化率分布χSMT；以及合成部345，其对上述更新磁化率分布χUPD和上述处理后磁化率分布χSMT进行合成，来生成合成磁化率分布χADD，上述判定部346也可以将判定为结束时的上述合成磁化率分布χADD设为磁化率分布χ。
并且，上述误差函数e(χ0)被设定为向对根据上述磁场分布δ和上述初始磁化率分布χ0进行卷积运算而得的结果决定的函数赋予加权W，上述加权W也可以基于根据上述复图像I获得的相位图像P计算出。
这样，在本实施方式中，以基于磁场和磁化率的关系式的制约条件，对平滑化处理进行反复运算，从而能够计算出高精度的磁化率分布。
即，在本实施方式中，在使误差函数最小化的最小化处理的过程中，对所获得的更新磁化率分布χUPD进行平滑化，反复进行处理。在平滑化过程中，在k空间上进行加权来合成通过平滑化滤波器进行平滑化而得的处理后磁化率分布χSMT和使用误差函数的梯度计算出的更新磁化率分布χUPD，以便能够降低噪声且维持微细结构。
这样，根据本实施方式，通过平滑化滤波器进行的滤波，能够使噪声部分的数据平均化，因此能够抑制所计算出的磁化率分布χ中的噪声。并且，通过合成处理，在该磁化率分布χ中能够维持微细结构。
根据本实施方式，如现有方法那样，不使用正规化项，因此可以通过简单的处理获得抑制噪声且维持微细结构的磁化率分布χ。即，根据本实施方式，不需要进行用于设定恰当的正规化参数的庞大的计算。
此外，根据本实施方式，当解决不良条件问题时，不需要对魔角区域近旁的值进行舍位。因此，还能够抑制因舍位而产生的伪影，还可以提高计算的磁化率分布的精度。
并且，为了将更新磁化率分布χUPD设为能够抑制噪声且维持微细结构的磁化率分布而使用的处理后磁化率分布χSMT、以及进行反复运算时更新初始磁化率分布χ0的合成磁化率分布χADD均是根据更新磁化率分布χUPD计算出的。即，是取决于相位信息的数据。
在本实施方式中，使用这样的取决于相位信息的数据来更新在最小化处理的反复运算中计算出的磁化率分布，因此具有难以产生所测量的相位信息和所计算出的磁化率分布的不匹配的问题。
并且，根据本实施方式，根据相位图像计算出误差函数内的加权W。相位图像与绝对值图像相比相位偏差的反应精度高，因此能够获得高精度地反应了该相位偏差的加权W。然后，通过使用了该加权W的误差函数计算出磁化率分布χ，因此能够更高精度地计算出磁化率分布。
如以上说明的那样，根据本实施方式，通过简单的处理，边维持磁场和磁化率的关系，边反复进行保存微细结构等边缘信息的平滑化处理，从而能够降低伪影和噪声，并且提高计算的磁化率分布的准确度和精度。
＜合成时权重图像的变形例＞
另外，合成部345对更新磁化率分布χUPD和处理后磁化率分布χSMT进行合成时所使用的权重图像G并不限定于上述的权重图像。例如，也可以分别将GUPD(k)和GSMT(k)设为固定值。即，也可以对各k空间数据内的全部数据设定一样的权重。
通过设定一样的权重，容易控制最终得到的图像的SNR。
＜误差函数内的权重的变形例＞
此外，在上述实施方式中，根据复图像I的相位值分量计算出由误差函数设定部341设定的误差函数e(χ0)内的加权W，但并不限定于此。例如，可以根据复图像I的绝对值分量计算。根据复图像I的绝对值分量计算，能够以较少的处理简单地计算出权重。
附图标记说明
100：MRI装置、101：被检体、102：静磁场线圈、103：梯度磁场线圈、104：匀场线圈、105：发送线圈、106：接收线圈、107：发送机、108：接收机、109：计算机、110：显示器、111：外部存储装置、112：梯度磁场用电源部、113：匀场用电源部、114：顺序控制装置、115：输入装置、120：MRI装置、130：MRI装置、310：测量控制部、320：图像重构部、330：磁场分布计算部、331：相位图像计算部、332：相位展开处理部、333：整体相位变化去除部、334：磁场分布变换部、340：磁化率分布计算部、341：误差函数设定部、342：更新磁化率分布计算部、343：平滑化处理部、344：滤波处理部、345：合成部、346：判定部、400：磁场分布、410：初始磁化率分布、420：更新磁化率分布、421：更新磁化率分布的k空间数据、422：对更新磁化率分布进行乘法运算的权重图像、430：处理后磁化率分布、431：处理后磁化率分布的k空间数据、432：对处理后磁化率分布进行乘法运算的权重图像、440：合成磁化率分布、501：切片梯度磁场脉冲、502：RF脉冲、503：切片编码梯度磁场脉冲、504：相位编码梯度磁场脉冲、505：读出梯度磁场脉冲、506：读出梯度磁场脉冲、507：回波、508：切片编码梯度磁场脉冲、509：相位编码梯度磁场脉冲、510：RSSG序列。

Claims (11)

1.一种磁共振成像装置，其特征在于，具备：
测量控制部，其测量与高频磁场脉冲的照射对应地从被检体产生的核磁共振信号作为复信号；
图像重构部，其根据上述复信号，对像素值为复数的复图像进行重构；
磁场分布计算部，其根据上述复图像生成相位图像，并根据该相位图像计算出磁场分布；以及
磁化率分布计算部，其根据上述磁场分布计算出磁化率分布，
上述磁化率分布计算部进行通过反复运算使误差函数最小化的最小化处理来计算出上述磁化率分布，并且将在该最小化处理中根据上述误差函数计算出的更新磁化率分布在每次计算更新磁化率分布时进行平滑化，其中，误差函数是使用预先决定的初始磁化率分布和上述磁场分布定义的函数，
上述磁化率分布计算部具备：
更新磁化率分布计算部，其使用上述误差函数，根据上述初始磁化率分布计算出上述更新磁化率分布；
平滑化处理部，其对上述更新磁化率分布进行平滑化，并生成合成磁化率分布；以及
判定部，其进行上述反复运算的结束判定，并且将上述初始磁化率分布更新为上述合成磁化率分布，
上述平滑化处理部具备：
滤波处理部，其对上述更新磁化率分布应用预先决定的平滑化滤波器，并生成处理后磁化率分布；以及
合成部，其对上述更新磁化率分布和上述处理后磁化率分布进行合成，生成合成磁化率分布，
上述判定部将判定为结束时的上述合成磁化率分布设为上述磁化率分布。
2.根据权利要求1所述的磁共振成像装置，其特征在于，
将上述误差函数设定为向根据对上述磁场分布和上述初始磁化率分布进行卷积运算而得的结果决定的函数赋予权重，
根据从上述复图像得到的相位图像计算出上述权重。
3.根据权利要求1所述的磁共振成像装置，其特征在于，
上述更新磁化率分布计算部计算出上述误差函数的梯度，并使用计算出的梯度，来计算出上述更新磁化率分布。
4.根据权利要求1所述的磁共振成像装置，其特征在于，
上述合成部将上述更新磁化率分布和上述处理后磁化率分布分别变换为k空间数据，对上述变换后的k空间数据分别赋予权重并进行加法运算来生成上述合成磁化率分布。
5.根据权利要求4所述的磁共振成像装置，其特征在于，
在上述更新磁化率分布的k空间数据中，将向魔角区域近旁的数据赋予的权重设定为小于向其他区域的数据赋予的权重，
在上述处理后磁化率分布的k空间数据中，将向魔角区域近旁的数据赋予的权重设定为大于其他区域的数据赋予的权重。
6.根据权利要求4所述的磁共振成像装置，其特征在于，
对各k空间数据内的全部数据，将上述权重设定为相同。
7.根据权利要求1所述的磁共振成像装置，其特征在于，
上述滤波器为边缘保存平滑化滤波器。
8.根据权利要求1所述的磁共振成像装置，其特征在于，
上述初始磁化率分布的全部值为0。
9.根据权利要求1所述的磁共振成像装置，其特征在于，
当计算出上述更新磁化率分布时，上述判定部使用该更新磁化率分布和上述初始磁化率分布进行上述结束判定。
10.根据权利要求1所述的磁共振成像装置，其特征在于，
当计算出上述合成磁化率分布时，上述判定部使用该合成磁化率分布和上述初始磁化率分布进行上述结束判定。
11.一种定量性磁化率匹配方法，其特征在于，
从收集并重构根据高频磁场脉冲的照射而从被检体产生的核磁共振信号而得的复图像获得相位图像，
根据上述相位图像计算出磁场分布，
进行通过反复运算使由预先决定的初始磁化率分布和上述磁场分布定义的误差函数最小化的最小化处理，计算出磁化率分布，
在上述最小化处理中，
对根据上述误差函数计算出的、更新上述初始磁化率分布而得的更新磁化率分布应用预先决定的平滑化滤波器而生成处理后磁化率分布，并对上述更新磁化率分布和上述处理后磁化率分布进行合成，生成合成磁化率分布，
将上述初始磁化率分布更新为该合成磁化率分布，并将反复运算结束时的上述合成磁化率分布设为上述磁化率分布。
CN201580060823.6A 2014-11-11 2015-10-19 磁共振成像装置以及定量性磁化率匹配方法 Active CN107072592B (zh)

Applications Claiming Priority (3)

Application Number Priority Date Filing Date Title
JP2014-228843 2014-11-11
JP2014228843 2014-11-11
PCT/JP2015/079477 WO2016076076A1 (ja) 2014-11-11 2015-10-19 磁気共鳴イメージング装置および定量的磁化率マッピング方法

Publications (2)

Publication Number Publication Date
CN107072592A CN107072592A (zh) 2017-08-18
CN107072592B true CN107072592B (zh) 2020-06-30

Family

ID=55954168

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201580060823.6A Active CN107072592B (zh) 2014-11-11 2015-10-19 磁共振成像装置以及定量性磁化率匹配方法

Country Status (4)

Country Link
US (1) US10180474B2 (zh)
JP (1) JP6289664B2 (zh)
CN (1) CN107072592B (zh)
WO (1) WO2016076076A1 (zh)

Families Citing this family (15)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP6702691B2 (ja) * 2015-10-30 2020-06-03 キヤノンメディカルシステムズ株式会社 磁気共鳴イメージング及び医用画像処理装置
EP3391069B1 (en) * 2015-12-17 2022-02-09 Koninklijke Philips N.V. Segmentation of quantitative susceptibility mapping magnetic resonance images
EP3497463B1 (en) * 2016-08-09 2022-03-30 Koninklijke Philips N.V. Retrospective correction of field fluctuations in multiple gradient echo mri
JP7408271B2 (ja) * 2018-01-17 2024-01-05 キヤノンメディカルシステムズ株式会社 磁気共鳴イメージング装置および医用画像処理装置
CN108523890B (zh) * 2018-03-06 2020-07-31 厦门大学 基于磁共振偶极场空间分布的局部场估计方法
EP3550318A1 (de) * 2018-04-06 2019-10-09 Siemens Healthcare GmbH Verfahren zur erzeugung einer b0-karte mittels eines magnetresonanz-fingerabdruckverfahrens
CN108714028B (zh) * 2018-04-11 2022-02-25 上海联影医疗科技股份有限公司 一种磁共振成像方法、装置及医学成像系统
CN108896943B (zh) * 2018-05-10 2020-06-12 上海东软医疗科技有限公司 一种磁共振定量成像方法和装置
CN108665465B (zh) * 2018-05-14 2019-07-19 东软医疗系统股份有限公司 分割缺血半暗带的方法、装置、存储介质及设备
EP3591418A1 (en) * 2018-07-03 2020-01-08 Koninklijke Philips N.V. Mri method for b0-mapping
CN109239607A (zh) * 2018-07-13 2019-01-18 深圳先进技术研究院 电池检测方法、系统及电池分析装置
WO2020010624A1 (zh) * 2018-07-13 2020-01-16 深圳先进技术研究院 电池检测方法、系统及电池分析装置
CN111481827B (zh) * 2020-04-17 2023-10-20 上海深透科技有限公司 定量磁化率成像及用于dbs潜在刺激靶区定位的方法
CN114062988B (zh) * 2020-07-31 2023-09-22 上海联影医疗科技股份有限公司 磁共振波谱成像方法、装置、计算机设备和存储介质
JPWO2022190249A1 (zh) * 2021-03-10 2022-09-15

Family Cites Families (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7782051B2 (en) * 2008-04-18 2010-08-24 Mr Innovations, Inc. Geometry based field prediction method for susceptibility mapping and phase artifact removal
EP2283373B1 (en) 2008-04-28 2021-03-10 Cornell University Accurate quantification of magnetic susceptibility in molecular mri
JP5600946B2 (ja) * 2010-01-28 2014-10-08 株式会社島津製作所 断層撮影装置
US8422756B2 (en) 2010-04-27 2013-04-16 Magnetic Resonance Innovations, Inc. Method of generating nuclear magnetic resonance images using susceptibility weighted imaging and susceptibility mapping (SWIM)
CN103764025B (zh) * 2011-10-12 2015-12-09 株式会社日立制作所 磁共振成像装置以及磁化率强调图像生成方法
WO2014057716A1 (ja) * 2012-10-10 2014-04-17 株式会社日立製作所 磁気共鳴イメージング装置
WO2014076808A1 (ja) * 2012-11-16 2014-05-22 株式会社日立製作所 磁気共鳴イメージング装置および定量的磁化率マッピング法
JP6072723B2 (ja) * 2014-04-21 2017-02-01 株式会社日立製作所 磁気共鳴イメージング装置、及び画像撮像方法

Also Published As

Publication number Publication date
US10180474B2 (en) 2019-01-15
JPWO2016076076A1 (ja) 2017-07-27
US20180059198A1 (en) 2018-03-01
WO2016076076A1 (ja) 2016-05-19
JP6289664B2 (ja) 2018-03-07
CN107072592A (zh) 2017-08-18

Similar Documents

Publication Publication Date Title
CN107072592B (zh) 磁共振成像装置以及定量性磁化率匹配方法
JP5843876B2 (ja) 磁気共鳴イメージング装置および磁化率強調画像生成方法
JP5902317B2 (ja) 磁気共鳴イメージング装置および定量的磁化率マッピング法
CN108778116B (zh) 磁共振成像装置以及图像处理方法
US10765349B2 (en) Magnetic resonance imaging device and method for calculating oxygen extraction fractions
JP5394374B2 (ja) 磁気共鳴イメージング装置及び血管画像取得方法
US9709641B2 (en) Magnetic resonance imaging apparatus, image processing apparatus, and susceptibility map calculation method
JP5982494B2 (ja) 磁気共鳴イメージング装置
US10203387B2 (en) MR imaging with enhanced susceptibility contrast
WO2017056996A1 (ja) 磁気共鳴イメージング装置、および、画像処理装置
US11033199B2 (en) Echo-planar imaging magnetic resonance elastography pulse sequence
JP7230149B2 (ja) 磁気共鳴イメージングで得た画像の処理方法、画像処理プロブラム、及び、計算機
US10782375B2 (en) Multi-contrast images from a magnetic resonance imaging scan
Zhang Improved ultrashort time echo and dynamic contrast enhancement magnetic resonance imaging based on Stack-of-Stars golden angle radial sampling scheme
JP2023541775A (ja) ボリュームを網羅する拡散強調磁気共鳴画像のシーケンスを取得し再構成する方法及び装置
CN114076912A (zh) 抑制静态组织的方法、磁共振成像方法及系统

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
TR01 Transfer of patent right
TR01 Transfer of patent right

Effective date of registration: 20211123

Address after: Chiba County, Japan

Patentee after: Fujifilm medical health Co.,Ltd.

Address before: Tokyo, Japan

Patentee before: Hitachi, Ltd.