WO2019049443A1 - 磁気共鳴イメージング装置 - Google Patents

磁気共鳴イメージング装置 Download PDF

Info

Publication number
WO2019049443A1
WO2019049443A1 PCT/JP2018/021266 JP2018021266W WO2019049443A1 WO 2019049443 A1 WO2019049443 A1 WO 2019049443A1 JP 2018021266 W JP2018021266 W JP 2018021266W WO 2019049443 A1 WO2019049443 A1 WO 2019049443A1
Authority
WO
WIPO (PCT)
Prior art keywords
phase distribution
magnetic field
unit
distribution
local
Prior art date
Application number
PCT/JP2018/021266
Other languages
English (en)
French (fr)
Inventor
亨 白猪
良太 佐藤
久晃 越智
毅倫 村瀬
尾藤 良孝
Original Assignee
株式会社日立製作所
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 株式会社日立製作所 filed Critical 株式会社日立製作所
Priority to US16/646,033 priority Critical patent/US11051695B2/en
Publication of WO2019049443A1 publication Critical patent/WO2019049443A1/ja

Links

Images

Classifications

    • AHUMAN NECESSITIES
    • A61MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61BDIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00Measuring for diagnostic purposes; Identification of persons
    • A61B5/0033Features or image-related aspects of imaging apparatus classified in A61B5/00, e.g. for MRI, optical tomography or impedance tomography apparatus; arrangements of imaging apparatus in a room
    • A61B5/004Features or image-related aspects of imaging apparatus classified in A61B5/00, e.g. for MRI, optical tomography or impedance tomography apparatus; arrangements of imaging apparatus in a room adapted for image acquisition of a particular organ or body part
    • A61B5/0042Features or image-related aspects of imaging apparatus classified in A61B5/00, e.g. for MRI, optical tomography or impedance tomography apparatus; arrangements of imaging apparatus in a room adapted for image acquisition of a particular organ or body part for the brain
    • 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
    • G01R33/56Image enhancement or correction, e.g. subtraction or averaging techniques, e.g. improvement of signal-to-noise ratio and resolution
    • AHUMAN NECESSITIES
    • A61MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61BDIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00Measuring for diagnostic purposes; Identification of persons
    • A61B5/05Detecting, measuring or recording for diagnosis by means of electric currents or magnetic fields; Measuring using microwaves or radio waves 
    • A61B5/055Detecting, 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
    • 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/443Assessment of an electric or a magnetic field, e.g. spatial mapping, determination of a B0 drift or dosimetry
    • 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/4828Resolving the MR signals of different chemical species, e.g. water-fat imaging
    • 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
    • G01R33/543Control of the operation of the MR system, e.g. setting of acquisition parameters prior to or during MR data acquisition, dynamic shimming, use of one or more scout images for scan plane prescription
    • 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
    • G01R33/56Image enhancement or correction, e.g. subtraction or averaging techniques, e.g. improvement of signal-to-noise ratio and resolution
    • G01R33/5608Data 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

Definitions

  • the present invention relates to magnetic resonance imaging (MRI) technology.
  • the present invention relates to an image processing technique for calculating a local phase distribution caused by a magnetic susceptibility difference between living tissues using a complex image captured by the Gradient Echo method.
  • a magnetic resonance imaging apparatus is a non-invasive medical image that utilizes nuclear magnetic resonance (NMR) phenomenon in which hydrogen nuclei (protons) of a subject placed in a static magnetic field resonate with a high frequency magnetic field of a specific frequency. It is a diagnostic device.
  • the nuclear magnetic resonance signal generated from the subject changes with various physical property values such as proton density and relaxation time. Therefore, the magnetic resonance imaging apparatus is based on the nuclear magnetic resonance signal, the structure and composition of the living tissue, cell characteristics, etc. Various biological information can be depicted.
  • NMR nuclear magnetic resonance
  • the magnetic susceptibility is a physical property value representing the degree to which a substance in a static magnetic field is magnetically polarized (magnetized).
  • the magnetic resonance imaging apparatus can image the magnetic susceptibility difference between living tissues, so it has the potential for diagnosis of cerebral ischemic disease, prediction of radiotherapy effect of cancer, differentiation of neurodegenerative disease, etc. It is expected.
  • a quantitative susceptibility mapping (QSM) method, a susceptibility weighted imaging (SWI) method, etc. are known. There is.
  • the magnetic resonance imaging apparatus calculates a local phase distribution (distribution of local phase change) caused by a magnetic susceptibility difference between living tissues, and calculates a quantitative magnetic susceptibility distribution by the QSM method using the local phase distribution. Alternatively, the susceptibility-emphasized image is calculated by the SWI method using the local phase distribution.
  • the magnetic resonance imaging apparatus performs phase aliasing correction processing to calculate the total phase distribution, and removes the background phase distribution generated due to the shape of the object or the like from the total phase distribution. Perform phase removal processing. The following two methods have been proposed as representative methods in background phase removal processing.
  • the first method is a method called PDF (Projection onto Dipole Filed) method (see, for example, Patent Document 1).
  • PDF Project onto Dipole Filed
  • the PDF method estimates the magnetic susceptibility distribution outside the brain from the measured phase distribution inside the brain, calculates the phase distribution inside the brain generated by the magnetic susceptibility distribution as a background phase distribution, and removes the background phase distribution from the total phase distribution By doing this, the local phase distribution is calculated.
  • the second method is a method called the SHARP (Sophisticated Harmonic Artifact Reduction for Phase date) method (see, for example, Non-Patent Document 1).
  • the SHARP method utilizes the fact that the background phase distribution can be approximated by spherical harmonics.
  • Patent Document 1 has a problem that the background phase distribution of a portion (for example, a paranasal sinus etc.) which can not be reproduced by estimated magnetic susceptibility outside the brain can not be removed. Further, the method described in Non-Patent Document 1 has a problem that the calculation accuracy of the local phase distribution at the boundary between the inside of the brain and the outside of the brain (brain surface area) is lowered.
  • a portion for example, a paranasal sinus etc.
  • the present invention has been made to solve the above-described problems, and an object of the present invention is to provide a magnetic resonance imaging apparatus that accurately calculates the local phase distribution of the entire brain including the brain surface region.
  • a magnetic resonance imaging apparatus includes a static magnetic field generating magnet that generates a static magnetic field in a space in which a subject is disposed, and a transmission that transmits a high frequency magnetic field to the subject
  • a receiving unit for receiving a nuclear magnetic resonance signal generated from the subject by transmission of the high frequency magnetic field; a gradient magnetic field applying unit applying a gradient magnetic field for adding positional information to the nuclear magnetic resonance signal;
  • a control unit that controls the unit, the reception unit, and the gradient magnetic field application unit, and a display unit that displays an image calculated by the computer.
  • the computer includes: a target tissue extraction unit for calculating a target tissue extraction image in which a target tissue of the subject is extracted from a complex image constituted by the nuclear magnetic resonance signal; A phase distribution calculation unit that calculates a distribution, a phase reflection correction unit that corrects phase reflection of the phase distribution, and a phase distribution after correction based on the target tissue extraction image, the predetermined region of the target tissue, and the predetermined
  • the target tissue is obtained by combining an area dividing unit for dividing into an area other than the area, a background phase distribution calculating unit for calculating the background phase distribution of the predetermined area, and the phase distribution after the correction and the background phase distribution.
  • a local phase distribution calculating unit that calculates a local phase distribution in a predetermined region.
  • (A) is an external view of a perpendicular magnetic field type magnetic resonance imaging apparatus according to the present embodiment
  • (b) is an external appearance view of a horizontal magnetic field type magnetic resonance imaging apparatus according to the present embodiment
  • (c) is this It is an outline view of a magnetic resonance imaging device which raised an open feeling concerning an embodiment.
  • It is a block diagram showing a schematic structure of a magnetic resonance imaging device concerning this embodiment.
  • It is a functional block diagram of the computer concerning this embodiment. It is a flowchart explaining the flow of the arithmetic processing in the computer which concerns on this embodiment. It is a figure for demonstrating the pulse sequence of a RSSG (RF-spoiled-Steady-state Acquisition with Rewound Gradient-Echo) sequence.
  • RSSG RF-spoiled-Steady-state Acquisition with Rewound Gradient-Echo
  • the weight is a diagram for explaining the weighting image W 1 for use in all areas local phase distribution calculating unit according to the present embodiment, (b), to be used in all areas local phase distribution calculating unit according to the present embodiment it is a diagram for explaining an image W 2.
  • (A) is a figure which shows the image to which arithmetic processing was performed by the computer which concerns on this embodiment, (b) is the enlarged view in the partial area
  • FIG. 1 is an external view of an MRI apparatus according to the present embodiment.
  • FIG. 1A shows a horizontal magnetic field type MRI apparatus 100 using a tunnel type magnet which generates a static magnetic field by a solenoid coil.
  • FIG. 1 (b) shows an open-type vertical magnetic field type MRI apparatus 120 in which magnets are vertically separated to enhance the sense of openness.
  • FIG. 1 (c) is an MRI apparatus 130 which uses the same tunnel type magnet as FIG. 1 (a), shortens the depth of the magnet, and inclines obliquely to enhance the sense of openness.
  • the configuration of the MRI apparatus is not particularly limited, and any configuration shown in FIG. 1 may be used, or a known configuration may be used.
  • the MRI apparatus 100 is representative.
  • the static magnetic field direction of the MRI apparatus 100 is the z direction, and of the two directions perpendicular to the z direction, the direction parallel to the bed surface on which the subject 101 to be measured is placed is x direction, y direction.
  • FIG. 2 is a block diagram showing a schematic configuration of the MRI apparatus according to the present embodiment.
  • the MRI apparatus 100 generates a static magnetic field generating magnet 102 that generates a static magnetic field in the space where the subject 101 is disposed, and a transmission high-frequency coil (transmission unit) that transmits a high frequency magnetic field to the subject 101 105, a high frequency reception coil (reception unit) 106 for receiving a nuclear magnetic resonance signal generated from the subject 101 by transmission of a high frequency magnetic field, and a gradient magnetic field application for applying a gradient magnetic field for adding positional information to the nuclear magnetic resonance signal
  • a computer 109 that controls the unit 103, the transmitter 105, the receiver 106, and the gradient magnetic field application unit 103, and a computer 109 that performs arithmetic processing on nuclear magnetic resonance signals, and a display that displays an image calculated by the computer 109.
  • the MRI apparatus 100 further includes a shim coil 104, a transmitter 107, a receiver 108, an external storage device 111, a gradient power supply 112, a shim power supply 113, a sequence control device 114, an input device 115, and the like.
  • the static magnetic field generating magnet 102 is a magnet that generates a static magnetic field in the space in which the subject 101 is disposed, and various forms can be adopted according to the configuration of the MRI apparatus.
  • the gradient magnetic field application unit 103 applies a gradient magnetic field in each of the x direction, y direction, and z direction in the space where the object 101 is disposed, in order to add spatial positional information to the nuclear magnetic resonance signal. It is a coil.
  • the gradient magnetic fields (x direction gradient magnetic field, y direction gradient magnetic field, z direction gradient magnetic field) applied in each direction are, for example, a phase encoding gradient magnetic field, a readout gradient magnetic field, and a slice gradient magnetic field.
  • the shim coil 104 is a current coil used to adjust the static magnetic field distribution, such as canceling out and removing the inhomogeneity of the static magnetic field generated in the space in which the subject 101 is disposed.
  • the transmission unit 105 is a transmission high-frequency coil that transmits a high-frequency magnetic field that excites hydrogen nuclei of the subject 101 to the subject 101.
  • the receiving unit 106 is a high frequency receiving coil for receiving a nuclear magnetic resonance signal generated from the subject 101 by transmitting a high frequency magnetic field.
  • the transmitting unit 105 and the receiving unit 106 may be configured separately or may be combined.
  • the transmitter 107 is controlled by the sequence control device 114 and generates a high frequency magnetic field to be transmitted to the transmission unit 105.
  • the receiver 108 is controlled by the sequence control device 114, detects a nuclear magnetic resonance signal received by the receiving unit 106, and transmits a complex signal to the computer 109.
  • the transmitting unit 105 and the receiving unit 106 may be configured separately or may be combined.
  • the computer 109 is an information processing apparatus including a central processing unit (CPU), a memory, a storage device, and the like.
  • the display unit 110, the external storage device 111, the input device 115, and the like are connected to the computer 109.
  • the computer 109 reads out a control program stored in the storage device, expands it in a work area, and executes the control program, thereby controlling the operation of the entire MRI apparatus 100 and performing various arithmetic processing.
  • the computer 109 may be a part of the MRI apparatus 100 or an information processing apparatus independent of the MRI apparatus 100, and may be an information processing apparatus capable of transmitting and receiving data with the MRI apparatus 100. .
  • all or part of the various functions realized by the computer 109 may be realized by hardware such as an application specific integrated circuit (ASIC) or a field-programmable gate array (FPGA).
  • ASIC application specific integrated circuit
  • FPGA field-programmable gate array
  • the display unit 110 is configured of, for example, a liquid crystal display, an organic EL display, or the like.
  • the display unit 110 displays the calculation result, various images, and the like on the display screen based on the data acquired from the computer 109.
  • the display unit 110 displays, for example, a local phase distribution described later, a quantitative magnetic susceptibility distribution, a magnetic susceptibility weighted image, and the like on the display screen.
  • the external storage device 111 is, together with the storage device, data used for various arithmetic processing performed by the computer 109, data generated during the arithmetic processing, data acquired by the various arithmetic processing, and input via the input device 115. Various conditions, various parameters, etc. are stored.
  • the gradient magnetic field power supply unit 112 is a power supply for driving the gradient magnetic field application unit 103, and controls the drive of the gradient magnetic field application unit 103 based on a control instruction from the sequence control device 114.
  • the shim power supply unit 113 is a power supply for driving the shim coil 104, and controls driving of the shim coil 104 based on a control instruction from the sequence control device 114.
  • the sequence control device 114 controls the operation of each unit of the MRI apparatus 100 such as the transmitter 107, the receiver 108, the gradient power supply unit 112, and the shim power supply unit 113.
  • the sequence control device 114 applies the high frequency magnetic field application timing and the gradient magnetic field (x direction gradient magnetic field) based on setting information (for example, a time chart called a pulse sequence described later) stored in the external storage device 111 or the storage device.
  • setting information for example, a time chart called a pulse sequence described later
  • Y-direction gradient magnetic field, z-direction gradient magnetic field) application timing, nuclear magnetic resonance signal reception timing, etc. are controlled.
  • the input device 115 has a function of receiving an input operation by the operator, and includes, for example, a keyboard, a camera, a mouse, a touch panel, a microphone that receives an audio input, and the like.
  • the input device 115 is an interface for the operator to input various conditions, various parameters, and the like necessary for various arithmetic processing and measurement processing performed by the computer 109. Note that various conditions, various parameters, and the like can be set in advance, or can be appropriately set by the operator via the input device 115.
  • the MRI apparatus calculates a phase distribution P from a complex image captured at an arbitrary echo time (time from transmission of a high frequency magnetic field to reception of a nuclear magnetic resonance signal) using the Gradient Echo method. Then, the MRI apparatus performs a phase aliasing correction process for correcting the phase which is folded in the range of ⁇ to + ⁇ occurring in the phase distribution P.
  • the phase distribution (total phase distribution) P total after the phase aliasing correction process uses a background phase distribution P bkgr attributed to the shape of the target tissue and the like, and a local phase distribution P local attributed to the magnetic susceptibility difference between the living tissues. It is expressed by equation (1). Then, the MRI apparatus uses the property of an average value that the pixel value at an arbitrary pixel and the average of the pixel values at pixels in a sphere (kernel) having a predetermined radius centered on that pixel are substantially equal. The background phase distribution P bkgr and the local phase distribution P local are separated from the phase distribution P total . The background phase distribution P bkgr satisfies the equation (2) according to the nature of the mean value.
  • is a delta function
  • r r is a spherical kernel of radius r
  • A is a convolution integral
  • M shrink represents a mask reduced according to the radius r of the spherical kernel r r from the target tissue extraction image M in which the target tissue of the subject 101 is extracted.
  • the convolution integral A represented by the equation (2) can be converted to the equation (3) by using the Fourier transform matrix F and the inverse Fourier transform matrix F ⁇ 1 .
  • C F ( ⁇ - ⁇ r )
  • the MRI apparatus calculates a local phase distribution P local satisfying the relational expression of equation (4) from the total phase distribution P total, and the local phase distribution P ′ local represented by equation (5) by least squares processing.
  • the MRI apparatus 100 divides the total phase distribution into the brain surface area and the area inside the brain, and combines the total phase distribution and the background phase distribution in the brain surface area to obtain the brain surface. Calculate the local phase distribution of the whole brain including the region. As a result, the calculation accuracy of the local phase distribution in the brain surface area can be improved.
  • the functional configuration of the computer 109 for realizing this will be described below.
  • the computer 109 includes a measurement control unit 310, an image reconstruction unit 320, a local phase distribution calculation unit 330, and a desired distribution calculation unit 340.
  • the local phase distribution calculating unit 330 includes a target tissue extracting unit 331, a phase distribution calculating unit 332, a phase aliasing correcting unit 333, an area dividing unit 334, a background phase distribution calculating unit 335, and an all area local phase distribution calculating unit 336; Equipped with
  • the measurement control unit 310 measures a nuclear magnetic resonance signal of any echo time generated from the subject 101 as a complex signal by transmitting a high frequency magnetic field from the transmission unit 105 to the subject 101 using, for example, the Gradient Echo method. Do.
  • the measurement control unit 310 controls the sequence control device 114 according to a preset pulse sequence to control the application timing of the high frequency magnetic field, the application timing of the gradient magnetic field, the reception timing of the nuclear magnetic resonance signal, and the like.
  • the measurement control unit 310 may apply a flow compensation gradient magnetic field pulse that compensates for the influence of flow such as blood flow in each of the x direction, the y direction, and the z direction.
  • the image reconstruction unit 320 arranges the complex signal (nuclear magnetic resonance signal) measured by the measurement control unit 310 in, for example, a three-dimensional k space having the kr axis, the kp axis, and the ks axis as coordinate axes to perform Fourier transform By doing this, a complex image in which pixel values of arbitrary pixels are represented by complex numbers is reconstructed.
  • the image reconstruction unit 320 uses the method of filling the nuclear magnetic resonance signal one row at a time in parallel with the coordinate axis (for example, the kr axis) of the k space.
  • Cartesian filling method a method may be used in which nuclear magnetic resonance signals are radially filled while changing the angle around the origin of k space (Radial Scan filling method).
  • the local phase distribution calculating unit 330 calculates the local phase distribution in the entire region of the target tissue based on the complex image reconstructed by the image reconstructing unit 320.
  • the target tissue extraction unit 331 generates a target tissue extraction image M in which the target tissue of the subject 101 is extracted from the complex image.
  • the phase distribution calculation unit 332 calculates the phase distribution P from the complex image.
  • the phase folding correction unit 333 corrects the phase folding of the phase distribution P, and calculates the total phase distribution P total .
  • the region division unit 334 divides the total phase distribution P total into at least two or more regions (for example, a predetermined region, a region other than the predetermined region).
  • the background phase distribution calculating unit 335 calculates the background phase distribution P edge in the predetermined area. Then, the all-region local phase distribution calculation unit 336 combines the total phase distribution P total and the background phase distribution P edge in the predetermined region to generate the local phase distribution P ′ local in the entire region of the target tissue including the predetermined region. Calculate
  • the desired distribution calculation unit 340 uses the local phase distribution P ′ local calculated by the local phase distribution calculation unit 330 to calculate a desired distribution such as a quantitative magnetic susceptibility distribution or a susceptibility-emphasized image, and displays the desired distribution. Output to unit 110.
  • the computer 109 extracts only the target tissue from the complex image, divides the total phase distribution into the predetermined region of the target tissue and a region other than the predetermined region, and obtains the background phase in the predetermined region.
  • the distribution is calculated, and the local phase distribution over the entire region of the target tissue is calculated using the background phase distribution.
  • FIG. 4 is a flowchart for explaining the flow of arithmetic processing in the computer 109 according to the present embodiment.
  • FIG. 5 is a diagram for explaining a pulse sequence of the RSSG sequence.
  • FIG. 6 is a flow chart for explaining the flow of the local phase distribution calculation process in the local phase distribution calculation unit 330 according to this embodiment.
  • FIG. 7 is a diagram for explaining a weight image used in the all-region local phase distribution calculation unit 336 according to the present embodiment.
  • FIG. 8 is a view showing an example of an image displayed on the display unit 110 according to the present embodiment.
  • step S1001 the computer 109 performs measurement processing.
  • the measurement control unit 310 controls the sequence controller 114 according to a preset pulse sequence. Then, the measurement control unit 310 measures a nuclear magnetic resonance signal of any echo time as a complex signal.
  • the RSSG (RF-spoiled-steady-state acquisition with Rewound Gradient-Echo) sequence 510 shown in FIG. 5 is a pulse sequence that is sensitive to local magnetic field changes caused by differences in magnetic susceptibility between living tissues.
  • the sequence used by the measurement control unit 310 for measurement is not limited to the RSSG sequence 510.
  • RF indicates a high frequency magnetic field
  • Gs indicates a slice gradient magnetic field
  • Gp indicates a phase encoding gradient magnetic field
  • Gr indicates a readout gradient magnetic field.
  • TE indicates an echo time
  • TR indicates a repetition time.
  • 501 is a slice gradient magnetic field pulse 502, a high frequency magnetic field pulse 502, a slice encode gradient magnetic field pulse 503, 508, a phase encode gradient magnetic field pulse 504, 509, a readout gradient magnetic field pulse 505, 506, nuclear magnetic resonance 507 Signal is shown.
  • the numbers below the hyphen indicate how many times the measurement control unit 310 performs repeated measurements, for example, “-2” is the second time that the measurement control unit 310 performs repeated measurements. It is shown that.
  • the measurement control unit 310 controls the sequence control device 114 according to the RSSG sequence 510 to make the slice gradient magnetic field pulse 501, the high frequency magnetic field pulse 502, the slice encode gradient magnetic field pulse 503, the phase encode gradient magnetic field pulse 504, and the readout gradient magnetic field By appropriately controlling the application timing of the pulse 505, readout gradient magnetic field pulse 506, slice encode gradient magnetic field pulse 508, phase encode gradient magnetic field pulse 509, etc., data necessary for image reconstruction is measured.
  • one nuclear magnetic resonance signal 507 is measured within the repetition time TR.
  • the hydrogen nucleus in the object 101 has a predetermined frequency in accordance with the static magnetic field strength.
  • a high frequency magnetic field pulse 502-1 matching the frequency is transmitted to excite the hydrogen nucleus of the object 101.
  • the gradient magnetic field applying unit 103 applies the slice gradient magnetic field pulse 501-1 in the z direction in the measurement region of the subject 101, and the magnetization in a predetermined slice (imaging position of a three-dimensional image) in the subject 101. Selectively excite.
  • the measurement control unit 310 changes the phase of the high frequency magnetic field pulse 502 by a predetermined angle (for example, 117 degrees, 122 degrees). For example, if the first measurement is repeated, the predetermined angle is 117 degrees. If the second measurement is repeated, the angle is 117 ⁇ 2 degrees. If the third measurement is repeated, the angle is 117 ⁇ 3 degrees. I assume. Thereby, when applying the next high frequency magnetic field pulse, it becomes possible to perform appropriate measurement without leaving an extra signal component.
  • a predetermined angle for example, 117 degrees, 122 degrees.
  • the gradient magnetic field application unit 103 applies a slice encoding gradient magnetic field pulse 503-1 in the z direction in the measurement region of the subject 101 to add positional information in the slice encoding direction to the phase of magnetization. Further, the gradient magnetic field application unit 103 applies the phase encoding gradient magnetic field pulse 504-1 in the x direction in the measurement region of the object 101, and adds positional information in the phase encoding direction to the phase of magnetization. That is, the gradient magnetic field application unit 103 adds positional information to the phase of magnetization in the z direction and the x direction prior to reading of the nuclear magnetic resonance signal.
  • the magnetic field strength changes every repetition time TR.
  • seven horizontal lines shown in FIG. 5 indicate gradient magnetic fields of different strengths, and indicate that the positive gradient magnetic field becomes stronger according to the upward arrow.
  • the magnetic field strength changes every repetition time TR.
  • seven horizontal lines shown in FIG. 5 indicate gradient magnetic fields of different strengths, and indicate that the positive gradient magnetic field becomes stronger according to the upward arrow.
  • FIG. 5 illustrates the case where the magnetic field strengths of the slice encoding gradient magnetic field pulse 503-1 and the phase encoding gradient magnetic field pulse 504-1 change in seven ways, the present invention is limited thereto. It is not a thing.
  • the gradient magnetic field application unit 103 applies a readout gradient magnetic field pulse for dephasing 505-1 in the y direction in the measurement region of the object 101 to disperse the phase of the magnetization in the pixel. Then, the gradient magnetic field applying unit 103 applies the readout gradient magnetic field pulse 506-1 in the y direction in the measurement region of the object 101, and adds positional information in the readout direction to the phase of the magnetization.
  • the measurement control unit 310 can measure one nuclear magnetic resonance signal 507 within the repetition time TR.
  • the gradient magnetic field application unit 103 applies the slice encoding gradient magnetic field pulse 508-1 for rephasing in the z direction in the measurement region of the subject 101, and further rephases in the x direction in the measurement region of the subject 101.
  • a phase encoding gradient magnetic field pulse 159-1 is applied to converge the phase of the magnetization in the dephased pixel.
  • the slice encoding gradient magnetic field pulse 508-1 changes in magnetic field strength every repetition time TR.
  • seven horizontal lines shown in FIG. 5 indicate gradient magnetic fields of different strengths, and indicate that the negative gradient magnetic field becomes stronger according to the downward arrow.
  • the magnetic field strength changes every repetition time TR.
  • seven horizontal lines shown in FIG. 5 indicate gradient magnetic fields of different strengths, and indicate that the negative gradient magnetic field becomes stronger according to the downward arrow.
  • the measurement control unit 310 repeats the above procedure while changing the phase of the high frequency magnetic field pulse 502, the intensity of the slice encoding gradient magnetic field pulses 503 and 508, the intensity of the phase encoding gradient magnetic field pulses 504 and 509, etc.
  • the nuclear magnetic resonance signal 507 required for the configuration is measured.
  • the absolute value component of the complex image obtained by the RSSG sequence 510 becomes a T1 (longitudinal relaxation time) emphasized image when the echo time TE is set short, and the phase dispersion in the pixel is set when the echo time TE is set long. It becomes a T2 * (apparent lateral relaxation time) weighted image reflected.
  • the echo time TE is set long, and imaging of a T2 * -weighted image is performed.
  • step S1002 the computer 109 performs an image reconstruction process.
  • the image reconstruction unit 320 reconstructs a complex image by arranging the nuclear magnetic resonance signals measured in step S1001 on k-space and performing Fourier transform. It does not specifically limit as an imaging method, Well-known methods, such as a Cartesian imaging, a non-Cartesian imaging, a multi echo echo planar imaging method, are applicable.
  • step S1003 the computer 109 performs local phase distribution calculation processing.
  • the local phase distribution calculating unit 330 calculates the local phase distribution in the entire region of the target tissue based on the complex image reconstructed in step S1002.
  • the target tissue extraction unit 331 extracts the target tissue of the object 101 based on the magnitude of the pixel value in the absolute value component (absolute value image) of the complex image configured by the image reconstruction unit 320.
  • a target tissue extraction image M is calculated.
  • the target tissue extraction unit 331 separates the noise component and the signal component from the histogram of the pixel values of the absolute value image using, for example, a known method such as the discriminant analysis method, and extracts the target tissue of the subject 101
  • a target tissue extraction image M (for example, an image in which the pixel value inside the brain is 1 and the pixel value outside the brain is 0) is calculated.
  • the target tissue extraction unit 331 removes a region (for example, a region such as subcutaneous fat) which is unnecessary for calculation of the magnetic susceptibility distribution from the complex image, using a known method such as, for example, a region growing method. Then, the target tissue extraction image M may be calculated by extracting only the substantial region of the target tissue, which is a region necessary for calculation of the magnetic susceptibility distribution.
  • a region for example, a region such as subcutaneous fat
  • step S1102 the phase distribution calculating unit 332 calculates the phase distribution P based on the declination components of the complex image reconstructed by the image reconstructing unit 320.
  • step S1103 the phase aliasing correction unit 333 corrects the phase aliasing in the range of - ⁇ to + ⁇ occurring in the phase distribution P calculated in step S1102 using, for example, a known method such as the region expansion method. And calculate the total phase distribution P total .
  • the region dividing unit 334 determines the total phase distribution P total into at least two or more regions, that is, a predetermined region of the target tissue, and the target tissue. Divided into areas other than the predetermined area of The predetermined region of the target tissue is, for example, a boundary portion between a region where the pixel value is 1 and a region where the pixel value is 0 in an image in which the pixel value inside the brain is 1 and the pixel value outside the brain is 0. It is a certain brain surface area. Further, the area other than the predetermined area of the target tissue is, for example, an area inside the brain excluding the brain surface area.
  • the region defined as the brain surface region is the radius r of the spherical kernel r r used when the all-region local phase distribution calculation unit 336 calculates the local phase distribution in the entire region of the target tissue in step S1106 described later. Determined by (kernel size).
  • the background phase distribution calculating unit 335 calculates the background phase distribution P edge in a predetermined region (for example, brain surface region) of the target tissue from the total phase distribution P total calculated in step S1103.
  • a predetermined region for example, brain surface region
  • a method of calculating the background phase distribution P edge in a predetermined region of the target tissue will be described.
  • the background phase distribution calculating unit 335 calculates the background phase distribution P edge in a predetermined region of the target tissue from the total phase distribution P total by local polynomial approximation processing.
  • the background phase distribution calculating unit 335 sets (x, y, z) the position of an arbitrary pixel in the local region, and calculates polynomial coefficients Ak, l, n that satisfy Expression (6).
  • K represents the order in the x-direction polynomial
  • L represents the order in the y-direction polynomial
  • N represents the order in the z-direction polynomial.
  • the background phase distribution P edge in the predetermined region of the target tissue can be calculated by Expression (7).
  • the background phase distribution calculating unit 335 calculates the background phase distribution P edge by performing Expression (7) on all the pixels in the predetermined region of the target tissue.
  • step S1106 the all-region local phase distribution calculation unit 336 combines the total phase distribution P total and the background phase distribution P edge in the predetermined region of the target tissue to obtain the local phase distribution P ′ local in the entire region of the target tissue. calculate.
  • a method of calculating the local phase distribution P ′ local in the entire region of the target tissue will be described.
  • the all-region local phase distribution calculation unit 336 calculates an estimated value P ′ local of the local phase distribution in the entire region of the target tissue by the constrained least squares process.
  • W 1 represents a weight image with a small weight in the brain surface region
  • W 2 represents a weight image with a large weight in the brain surface region.
  • the first term in Equation (8) represents a constraint term based on the nature of the mean value in the region inside the brain, as in the SHARP method. For the constraint terms that satisfy the condition of the nature of the mean value, the weight in the region inside the brain is increased and the weight in the brain surface region is decreased.
  • the second term in equation (8) represents a constraint term based on the background phase distribution in the brain surface area. For the constraint terms based on the background phase distribution in the brain surface area, the weight in the brain surface area is increased and the weight in the area inside the brain is decreased.
  • FIGS. 7 (a) shows a one-dimensional profile of the weighted image W 1.
  • FIG. 7 (b) shows a one-dimensional profile of the weighted image W 2.
  • the horizontal axis shows the position, and the vertical axis shows the weight.
  • the combined area of the area 701 and the area 702 indicates the brain surface area
  • the area 703 indicates the area inside the brain excluding the brain surface area (an area other than the brain surface area).
  • the combined area of the area 704 and the area 705 indicates the brain surface area
  • the area 706 indicates the area inside the brain excluding the brain surface area (an area other than the brain surface area).
  • the area 701 (area 704) and the area 702 (area 705) are set to be larger than the radius r of the spherical kernel r r .
  • the all-region local phase distribution calculation unit 336 sets the weight in the region 701 to be 0, sets the weight in the region 703 to be 1, and The weight in the transition region 702 is set to be a value between 0 and 1 (linearly changing value). Further, the whole region local phase distribution calculation unit 336 sets the weight in the region 704 to be 1 and sets the weight in the region 706 to be 0, and makes a transition from the region 704 to the region 706 The weight at 705 is set to be a value between 0 and 1 (linearly changing value).
  • the all-region local phase distribution calculation unit 336 increases the weight in the region inside the brain and reduces the weight in the brain surface region, thereby achieving equation (8).
  • the constraint term based on the nature of the average value shown in can be calculated with the same accuracy as the SHARP method (see equation (5)).
  • the total area local phase distribution calculating unit 336 by using the weighted image W 1 and the weighted image W 2, to increase the weight of the brain surface area, to reduce the weight of the brain inside the area, the formula (8)
  • the constraint term based on the background phase distribution in the brain surface area shown in can be approximated to the background phase distribution in the brain surface area calculated by the local polynomial approximation process.
  • the whole region local phase distribution calculation unit 336 linearly changes the weight in the region 702 and the region 706 between 0 and 1 using the weight image W 1 and the weight image W 2 to obtain a brain surface region.
  • the combined area of the area 701 and the area 702, the combined area of the area 704 and the area 705) and the area inside the brain (the area 703 and the area 706) can be connected smoothly.
  • the all-region local phase distribution calculation unit 336 performs not only the constrained least squares processing as shown in equation (8) but also the least squares processing with regularization as shown in equation (9) (regularization of local phase distribution A condition that minimizes a term may be used to calculate the local phase distribution in the entire region of the target tissue (entire brain including the brain surface region).
  • the all-region local phase distribution calculation unit 336 inserts a regularization term in consideration of variations in pixel values of the local phase distribution P local , and implements equation (9) to implement locality in the entire region of the target tissue. It is also possible to calculate the phase distribution.
  • ⁇ 1 is a regularization parameter for adjusting the size of the constraint in a predetermined region of the target tissue.
  • ⁇ 2 is a regularization parameter for allowing the local phase distribution P local to some extent.
  • the first term in Equation (9) represents a constraint term based on the nature of the average value in the region inside the brain, as in the SHARP method.
  • the second term in equation (9) represents a constraint term based on the background phase distribution in the brain surface area.
  • the third term in the equation (9) represents a regularization term having a role of preventing excessive removal of the local phase change component from the phase distribution.
  • the calculation accuracy of the local phase distribution in the region inside the brain maintains the calculation accuracy substantially the same as the conventional SHARP method
  • the calculation accuracy of the local phase distribution in the brain surface area can be significantly improved as compared with the conventional case. That is, according to the arithmetic processing of the computer 109 according to the present embodiment, the brain including the brain surface area from the restriction term based on the property of the average value in the area inside the brain and the restriction term based on the background phase distribution in the brain surface area It has been shown that the overall local phase distribution can be calculated with high accuracy.
  • the desired distribution calculation unit 340 calculates a desired distribution such as a quantitative magnetic susceptibility distribution or a magnetic susceptibility weighted image by the QSM method or the SWI method using the local phase distribution P local calculated in step S1003.
  • the desired distribution calculating unit 340 includes a local magnetic field calculating unit and a magnetic susceptibility distribution calculating unit, and calculates a quantitative magnetic susceptibility distribution by the QSM method.
  • the desired distribution calculation unit 340 includes a mask generation unit and a magnetic susceptibility weighted image generation unit, and calculates a magnetic susceptibility weighted image by the SWI method.
  • Equation 10 the relative magnetic field distribution generated by the difference in magnetic susceptibility between the living tissues is expressed by equation (10).
  • g represents the nuclear gyromagnetic ratio of protons
  • B 0 represents the static magnetic field strength.
  • the local magnetic field distribution is expressed by equation (11) according to the Maxwell equation for the static magnetic field.
  • d (r) is the local magnetic field distribution
  • c (r) is the magnetic susceptibility distribution in the living body
  • a is the angle between the vector (r'-r) and the static magnetic field direction
  • d (r) is the point dipole magnetic field Represent each.
  • equation (11) can be transformed to equation (12) by Fourier transforming both sides of equation (11).
  • k (k x , k y , k z ) represents a position vector in k space.
  • D (k) represents the Fourier component of the local magnetic field distribution d (r).
  • D (k) represents the Fourier component of the point dipole magnetic field d (r).
  • X (k) represents the Fourier component of the susceptibility distribution c (r) in the living body.
  • SWIM Seceptibility Weighted Imaging and Mapping
  • the desired distribution calculating unit 340 may use any of these methods.
  • the desired distribution calculation unit 340 creates a predetermined susceptibility emphasizing mask, and multiplies the susceptibility emphasizing mask a predetermined number of times. Thereafter, the desired distribution calculation unit 340 calculates the susceptibility-emphasized image by multiplying the absolute value image of any echo time by the susceptibility-emphasized mask.
  • the susceptibility-emphasized image is an image in which a portion with high susceptibility is dark and a portion with low susceptibility is bright. The operator can determine whether the magnetic susceptibility is high or low based on the image.
  • the desired distribution calculating unit 340 can calculate the susceptibility-emphasized image using various known methods, and is not limited to the SWI method.
  • the computer 109 causes the quantitative susceptibility distribution and the susceptibility-emphasized image calculated in step S1004 to be displayed on the display screen of the display unit 110.
  • the computer 109 can store various data calculated in step S 1004 in the external storage device 111 or the storage device, and can display a quantitative magnetic susceptibility distribution or a susceptibility-emphasized image on a desired display device. .
  • FIG. 8A is an image (magnetic susceptibility weighted image) subjected to arithmetic processing by the computer 109 according to the present embodiment
  • FIG. 8B is an enlarged view of a region 801 in FIG. 8A. is there.
  • FIG. 8 (c) is an image (magnetic susceptibility weighted image) subjected to arithmetic processing by a conventional computer
  • FIG. 8 (d) is an enlarged view of a region 802 of FIG. 8 (c).
  • FIGS. 8 (a) and 8 (b) Comparing the images shown in FIGS. 8 (a) and 8 (b) with the images shown in FIGS. 8 (c) and 8 (d), the images shown in FIGS. 8 (a) and 8 (b) show: It can be seen that while the brain surface area is clearly imaged, in the images shown in FIGS. 8C and 8D, the brain surface area is imaged as surrounded by a thick solid line. That is, it is understood that the computer 109 according to the present embodiment can calculate more information by the area surrounded by the thick solid line, as compared with the conventional computer.
  • the operator diagnoses the structure, composition, and the like of the living tissue in the entire brain, it does not use the images shown in FIGS. 8 (c) and 8 (d), but does not use FIGS. By using the image shown in), more appropriate diagnosis can be performed. Furthermore, when diagnosing a disease that has occurred on the surface of the brain, such as subarachnoid hemorrhage or cerebral infarction, the operator performs an accurate diagnosis by using the images shown in FIGS. 8 (a) and 8 (b). be able to. That is, when the computer 109 according to the present embodiment performs predetermined arithmetic processing, the operator can acquire a highly reliable image in the brain surface area.
  • the computer 109 extracts, for example, a brain area from the complex image, divides the total phase distribution into a brain surface area and an area inside the brain, and outputs the background phase in the brain surface area.
  • the distribution is calculated, and the local phase distribution of the whole brain is calculated by combining the total phase distribution and the background phase distribution.
  • the MRI apparatus it is possible to calculate a quantitative magnetic susceptibility distribution or a magnetic susceptibility-emphasized image by using a highly accurate local phase distribution. As a result, the operator can acquire an image in which the clear area is expanded, and the reliability of the MRI apparatus can be enhanced.

Landscapes

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

Abstract

磁気共鳴イメージング装置(100)は、被検体(101)が配置される空間に静磁場を生成する静磁場生成磁石(102)と、前記被検体に高周波磁場を送信する送信部(105)と、前記高周波磁場の送信により前記被検体から発生する核磁気共鳴信号を受信する受信部(106)と、前記核磁気共鳴信号に位置情報を付加する傾斜磁場を印加する傾斜磁場印加部(103)と、前記送信部、前記受信部、及び前記傾斜磁場印加部を制御するとともに、前記核磁気共鳴信号に対して演算処理を行う計算機(109)と、前記計算機により演算処理された画像を表示する表示部(110)と、を備え、前記計算機は、所定の演算処理を行う。

Description

磁気共鳴イメージング装置
 本発明は、磁気共鳴イメージング(MRI:magnetic resonance imaging)技術に関する。特に、Gradient Echo法で撮像された複素画像を用いて、生体組織間の磁化率差に起因した局所位相分布を算出する画像処理技術に関する。
 磁気共鳴イメージング装置は、静磁場内に配置される被検体の水素原子核(プロトン)が、特定周波数の高周波磁場に共鳴する核磁気共鳴(NMR:nuclear magnetic resonance)現象を利用した非侵襲な医用画像診断装置である。被検体から発生する核磁気共鳴信号は、プロトン密度や緩和時間など様々な物性値によって変化するため、磁気共鳴イメージング装置は、核磁気共鳴信号に基づいて、生体組織の構造や組成、細胞性状などの様々な生体情報を描出できる。
 近年、磁気共鳴イメージング装置により測定可能な物性値として、生体組織間の磁化率差が注目されている。磁化率とは、静磁場中の物質が磁気分極(磁化)する度合いを表した物性値である。生体内には、静脈血中のデオキシヘモグロビンや鉄たんぱく質などの常磁性体と、生体組織の大部分を占める水や石灰化の基となるカルシウムなどの反磁性体とが存在する。磁気共鳴イメージング装置は、生体組織間の磁化率差を画像化することができるため、脳虚血疾患の診断やがんの放射線治療効果の予測、神経変性疾患の鑑別、などへの可能性が期待されている。生体組織間の磁化率差を画像化する方法としては、例えば、定量的磁化率マッピング(QSM:quantitative susceptibility mapping)法、磁化率強調画像化(SWI:Susceptibility Weighted Imaging)法、などが知られている。
 磁気共鳴イメージング装置は、生体組織間の磁化率差に起因する局所位相分布(局所的な位相変化の分布)を算出し、当該局所位相分布を用いてQSM法により定量的磁化率分布を算出する、或いは、当該局所位相分布を用いてSWI法により磁化率強調画像を算出する。磁気共鳴イメージング装置は、当該局所位相分布を算出する際、位相折り返し補正処理を行って総位相分布を算出し、総位相分布から被検体の形状などに起因して生じる背景位相分布を除去する背景位相除去処理を行う。背景位相除去処理における代表的な方法として、以下の二つの方法が提案されている。
 一つ目の方法は、PDF(Projection onto Dipole Fieled)法と呼ばれる方法である(例えば、特許文献1参照)。PDF法は、計測した脳内部の位相分布から脳外部の磁化率分布を推定し、該磁化率分布によって生じる脳内部の位相分布を背景位相分布として算出し、総位相分布から背景位相分布を除去することで局所位相分布を算出する。
 二つ目の方法は、SHARP(Sophisticated Harmonic Artifact Reduction for Phase date)法と呼ばれる方法である(例えば、非特許文献1参照)。SHARP法は、背景位相分布が球面調和関数近似できることを利用する。即ち、任意の点における値が、その点を中心とした任意の半径を有する球(カーネル)内の点における値の平均に等しい、という平均値の性質(Spherical mean value property)を用いて、総位相分布から背景位相分布と局所位相分布とを分離する。
 しかしながら、特許文献1に記載の方法では、脳外部の推定磁化率で再現できない部分(例えば、副鼻腔など)の背景位相分布を除去することができないという問題がある。また、非特許文献1に記載の方法では、脳内部と脳外部との境界部分(脳表領域)における局所位相分布の算出精度が低下してしまうという問題がある。
 本発明は、上記した課題を解決するためになされたものであり、脳表領域を含めた脳全体の局所位相分布を高精度に算出する磁気共鳴イメージング装置を提供することを課題とする。
 前記課題を解決するために、本開示の実施形態に係る磁気共鳴イメージング装置は、被検体が配置される空間に静磁場を生成する静磁場生成磁石と、前記被検体に高周波磁場を送信する送信部と、前記高周波磁場の送信により前記被検体から発生する核磁気共鳴信号を受信する受信部と、前記核磁気共鳴信号に位置情報を付加する傾斜磁場を印加する傾斜磁場印加部と、前記送信部、前記受信部、及び前記傾斜磁場印加部を制御するとともに、前記核磁気共鳴信号に対して演算処理を行う計算機と、前記計算機により演算処理された画像を表示する表示部と、を備え、前記計算機は、前記核磁気共鳴信号によって構成された複素画像から前記被検体の対象組織を抽出した対象組織抽出画像を算出する対象組織抽出部と、前記複素画像から位相分布を算出する位相分布算出部と、前記位相分布の位相折り返しを補正する位相折り返し補正部と、前記対象組織抽出画像に基づいて、補正後の位相分布を、前記対象組織の所定領域と前記所定領域以外の領域とに分割する領域分割部と、前記所定領域の背景位相分布を算出する背景位相分布算出部と、前記補正後の位相分布と前記背景位相分布とを組み合わせて、前記対象組織の所定領域における局所位相分布を算出する局所位相分布算出部と、を備えることを特徴とする。
 脳表領域を含めた脳全体の局所位相分布を高精度に算出する磁気共鳴イメージング装置を提供できる。
(a)は、本実施形態に係る垂直磁場方式の磁気共鳴イメージング装置の外観図、(b)は、本実施形態に係る水平磁場方式の磁気共鳴イメージング装置の外観図、(c)は、本実施形態に係る開放感を高めた磁気共鳴イメージング装置の外観図である。 本実施形態に係る磁気共鳴イメージング装置の概略構成を示すブロック図である。 本実施形態に係る計算機の機能ブロック図である。 本実施形態に係る計算機における演算処理の流れを説明するフローチャートである。 RSSG(RF-spoiled-Steady-state Acquisition with Rewound Gradient-Echo)シーケンスのパルスシーケンスを説明するための図である。 本実施形態に係る局所位相分布算出部における局所位相分布算出処理の流れを説明するフローチャートである。 (a)は、本実施形態に係る全領域局所位相分布算出部で用いる重み画像Wを説明するための図、(b)は、本実施形態に係る全領域局所位相分布算出部で用いる重み画像Wを説明するための図である。 (a)は、本実施形態に係る計算機によって演算処理が施された画像を示す図、(b)は、(a)の一部領域における拡大図、(c)は、従来の計算機によって演算処理が施された画像を示す図、(d)は、(c)の一部領域における拡大図である。
 本発明を適用する実施形態の一例を説明する。以下、実施形態を説明するための全図において、特に断らない限り、同一機能を有するものは同一符号を付し、繰り返しの説明は省略する。また、以下の記述により本発明が限定されるものではない。
 ≪MRI装置の全体構成≫
 図1は、本実施形態に係るMRI装置の外観図である。図1(a)は、ソレノイドコイルで静磁場を生成するトンネル型磁石を用いた水平磁場方式のMRI装置100である。図1(b)は、開放感を高めるために磁石を上下に分離したオープン型の垂直磁場方式のMRI装置120である。また、図1(c)は、図1(a)と同じトンネル型磁石を用い、磁石の奥行を短くし、かつ、斜めに傾けることによって、開放感を高めたMRI装置130である。なお、本実施形態において、MRI装置の構成は、特に限定されるものではなく、図1に示す何れの構成であってもよいし、公知の構成であってもよい。
 以下、特に区別する必要がない場合は、MRI装置100で代表する。MRI装置100において、MRI装置100の静磁場方向をz方向、z方向に垂直な2方向のうち、測定対称である被検体101を載置するベッド面に平行な方向をx方向、他方向をy方向とする。
 図2は、本実施形態に係るMRI装置の概略構成を示すブロック図である。図2に示すように、MRI装置100は、被検体101が配置される空間に静磁場を生成する静磁場生成磁石102と、被検体101に高周波磁場を送信する送信用高周波コイル(送信部)105と、高周波磁場の送信により被検体101から発生する核磁気共鳴信号を受信する受信用高周波コイル(受信部)106と、核磁気共鳴信号に位置情報を付加する傾斜磁場を印加する傾斜磁場印加部103と、送信部105、受信部106、及び傾斜磁場印加部103を制御するとともに、核磁気共鳴信号に対して演算処理を行う計算機109と、計算機109により演算処理された画像を表示する表示部110と、を備える。
 更に、MRI装置100は、シムコイル104、送信機107、受信機108、外部記憶装置111、傾斜磁場用電源部112、シム用電源部113、シーケンス制御装置114、入力装置115、などを備える。
 静磁場生成磁石102は、被検体101が配置される空間に静磁場を生成する磁石であり、MRI装置の構成に応じて、種々の形態を採用することができる。傾斜磁場印加部103は、核磁気共鳴信号に空間的な位置情報を付加するために、被検体101が配置される空間において、x方向、y方向、z方向の各方向に傾斜磁場を印加するコイルである。各方向に印加される傾斜磁場(x方向傾斜磁場、y方向傾斜磁場、z方向傾斜磁場)は、例えば、位相エンコード傾斜磁場、リードアウト傾斜磁場、スライス傾斜磁場である。
 シムコイル104は、被検体101が配置される空間に生成される静磁場の不均一性を打ち消して除去するなど、静磁場分布を調整するために使用される電流コイルである。送信部105は、被検体101の水素原子核を励起させる高周波磁場を、被検体101に送信する送信用高周波コイルである。受信部106は、高周波磁場の送信により被検体101から発生する核磁気共鳴信号を受信する受信用高周波コイルである。なお、送信部105と受信部106とは、別々の構成としてもよいし、兼用する構成としてもよい。
 送信機107は、シーケンス制御装置114によって制御され、送信部105へと送信する高周波磁場を生成する。受信機108は、シーケンス制御装置114によって制御され、受信部106によって受信された核磁気共鳴信号を検波し、複素信号を計算機109へと送信する。なお、送信部105と受信部106とは、別々の構成としてもよいし、兼用する構成としてもよい。
 計算機109は、CPU(Central Processing Unit)、メモリ、記憶装置などを備える情報処理装置であり、計算機109には、表示部110、外部記憶装置111、入力装置115などが接続される。計算機109は、記憶装置に記憶される制御プログラムを読み出して、ワークエリアに展開し、当該制御プログラムを実行することで、MRI装置100全体の動作を制御するとともに、各種演算処理を行う。なお、計算機109は、MRI装置100の一部であってもよいし、MRI装置100から独立した情報処理装置であって、MRI装置100とデータの送受信が可能な情報処理装置であってもよい。更に、計算機109が実現する各種機能のうち、全部または一部の機能は、ASIC(Application Specific Integrated Circuit)、FPGA(field-programmable gate array)などのハードウェアによって実現されていてもよい。
 表示部110は、例えば、液晶ディスプレイ、有機ELディスプレイ、等で構成される。表示部110は、計算機109から取得したデータに基づいて、ディスプレイ画面上に、演算結果、各種画像などを表示する。表示部110は、例えば、後述する局所位相分布、定量的磁化率分布、磁化率強調画像などをディスプレイ画面上に表示する。
 外部記憶装置111は、記憶装置とともに、計算機109が行う各種演算処理に使用されるデータ、演算処理中に生成されるデータ、各種演算処理により取得されるデータ、入力装置115を介して入力される各種条件、各種パラメータ、などを記憶する。
 傾斜磁場用電源部112は、傾斜磁場印加部103の駆動用電源であり、シーケンス制御装置114からの制御指示に基づいて、傾斜磁場印加部103の駆動を制御する。シム用電源部113は、シムコイル104の駆動用電源であり、シーケンス制御装置114からの制御指示に基づいて、シムコイル104の駆動を制御する。
 シーケンス制御装置114は、送信機107、受信機108、傾斜磁場用電源部112、シム用電源部113、などMRI装置100の各部の動作を制御する。例えば、シーケンス制御装置114は、外部記憶装置111や記憶装置に記憶される設定情報(例えば、後述のパルスシーケンスと呼ばれるタイムチャート)に基づいて、高周波磁場の印加タイミング、傾斜磁場(x方向傾斜磁場、y方向傾斜磁場、z方向傾斜磁場)の印加タイミング、核磁気共鳴信号の受信タイミング、などを制御する。
 入力装置115は、操作者による入力操作を受け付ける機能を備えており、例えば、キーボード、カメラ、マウス、タッチパネル、音声入力を受け付けるマイク、などで構成される。入力装置115は、計算機109が行う各種演算処理や計測処理などに必要な各種条件、各種パラメータ、などを操作者が入力するためのインタフェースである。なお、各種条件、各種パラメータ、などは、予め設定しておくことも可能であるし、操作者が入力装置115を介して適宜設定することも可能である。
 ≪計算機の構成≫
 次に、図3を参照して、本実施形態に係る計算機109の構成について説明する。本実施形態では対象組織として被検体101の脳を適用する場合を一例に挙げて説明する。
 ここで、計算機109の構成の説明に先立ち、SHARP法を用いた一般的な背景位相除去処理の概要について説明する。
 まず、MRI装置は、Gradient Echo法を用いて、任意のエコー時間(高周波磁場を送信してから核磁気共鳴信号を受信するまでの時間)で撮像された複素画像から位相分布Pを算出する。そして、MRI装置は、位相分布Pに生じている-π~+πの範囲に折り返された位相を補正する位相折り返し補正処理を行う。
 位相折り返し補正処理後の位相分布(総位相分布)Ptotalは、対象組織の形状などに起因する背景位相分布Pbkgr、生体組織間の磁化率差に起因する局所位相分布Plocalを用いて、式(1)で表される。
Figure JPOXMLDOC01-appb-M000001
 
 そして、MRI装置は、任意画素における画素値と、その画素を中心とした所定の半径を有する球(カーネル)内の画素における画素値の平均とが略等しいという平均値の性質を用いて、総位相分布Ptotalから、背景位相分布Pbkgrと局所位相分布Plocalとを分離する。
 背景位相分布Pbkgrは、平均値の性質によって、式(2)を満たす。
Figure JPOXMLDOC01-appb-M000002
 
 δはデルタ関数、ρは半径rの球状カーネル、Aは畳み込み積分をそれぞれ表す。また、Mshrinkは、被検体101の対象組織を抽出した対象組織抽出画像Mから、球状カーネルρの半径rに応じて縮小されたマスクを表す。
 式(2)で表される畳み込み積分Aは、フーリエ変換行列F、逆フーリエ変換行列F-1を用いることで、式(3)に変換できる。
Figure JPOXMLDOC01-appb-M000003
 但し、C=F(δ-ρ)
 式(3)に式(1)を代入することで、式(4)の関係式が得られる。
Figure JPOXMLDOC01-appb-M000004
 
 そして、MRI装置は、総位相分布Ptotalから、式(4)の関係式を満たす局所位相分布Plocalを算出し、最小二乗処理によって、式(5)で表される局所位相分布P’localを算出する。
Figure JPOXMLDOC01-appb-M000005
 
 即ち、SHARP法を用いた一般的な背景位相除去処理では、マスクMshrinkの範囲における局所位相分布しか算出することができず、脳内部と脳外部との境界部分(脳表領域)における局所位相分布の算出精度が低下するという問題が生じていた。
 そこで、本実施形態に係るMRI装置100は、総位相分布を、脳表領域と脳内部の領域とに領域分割し、総位相分布と脳表領域における背景位相分布とを組み合わせることで、脳表領域を含めた脳全体の局所位相分布を算出する。これにより、脳表領域における局所位相分布の算出精度を向上させることができる。以下、これを実現する計算機109の機能構成について説明する。
 計算機109は、計測制御部310と、画像再構成部320と、局所位相分布算出部330と、所望分布算出部340と、を備える。更に、局所位相分布算出部330は、対象組織抽出部331、位相分布算出部332、位相折り返し補正部333、領域分割部334、背景位相分布算出部335、全領域局所位相分布算出部336と、を備える。
 計測制御部310は、例えば、Gradient Echo法を用いて、送信部105から被検体101への高周波磁場の送信により、被検体101から発生する任意のエコー時間の核磁気共鳴信号を複素信号として計測する。計測制御部310は、予め設定されたパルスシーケンスに従って、シーケンス制御装置114を制御することで、高周波磁場の印加タイミング、傾斜磁場の印加タイミング、核磁気共鳴信号の受信タイミング、などを制御する。なお、計測制御部310は、血流などの流れの影響を補償する流速補正(Flow Compensation)傾斜磁場パルスを、x方向、y方向、z方向の各方向に印加してもよい。
 画像再構成部320は、計測制御部310によって計測された複素信号(核磁気共鳴信号)を、例えば、kr軸、kp軸、ks軸を座標軸とする3次元のk空間に配置してフーリエ変換することにより、任意画素の画素値が複素数で表される複素画像を再構成する。画像再構成部320は、k空間に核磁気共鳴信号を充填する際、k空間の座標軸(例えば、kr軸)に平行に、1行ずつ核磁気共鳴信号を充填していく方法を用いてもよいし(Cartesian充填法)、k空間の原点を中心に、角度を変化させながら放射状に核磁気共鳴信号を充填していく方法を用いてもよい(Radial Scan充填法)。
 局所位相分布算出部330は、画像再構成部320によって再構成された複素画像に基づいて、対象組織の全領域における局所位相分布を算出する。具体的には、対象組織抽出部331が、複素画像から被検体101の対象組織を抽出した対象組織抽出画像Mを生成する。そして、位相分布算出部332が、複素画像から位相分布Pを算出する。そして、位相折り返し補正部333が、位相分布Pの位相折り返しを補正し、総位相分布Ptotalを算出する。そして、領域分割部334が、対象組織抽出画像Mに基づいて、総位相分布Ptotalを少なくとも2つ以上の領域(例えば、所定領域、所定領域以外の領域)に分割する。そして、背景位相分布算出部335が、所定領域における背景位相分布Pedgeを算出する。そして、全領域局所位相分布算出部336が、総位相分布Ptotalと、所定領域における背景位相分布Pedgeとを組み合わせることで、所定領域を含めた対象組織の全領域における局所位相分布P’localを算出する。
 所望分布算出部340は、局所位相分布算出部330によって算出された局所位相分布P’localを用いて、定量的磁化率分布や磁化率強調画像などの所望分布を算出し、当該所望分布を表示部110へと出力する。
 上述のように、本実施形態に係る計算機109は、複素画像から対象組織のみを抽出し、総位相分布を対象組織の所定領域と所定領域以外の領域とに分割して、所定領域における背景位相分布を算出し、当該背景位相分布を利用して、対象組織の全領域における局所位相分布を算出する。これにより、カーネルの半径に応じて脳表領域で誤差が大きくなり、脳表領域における局所位相分布の算出精度が低下するという従来のSHARP法に生じていた問題点を回避し、脳表領域を含めた脳全体における局所位相分布を高精度に算出することができる。
 ≪計算機における演算処理の流れ≫
 次に、図4乃至図8を参照して、本実施形態に係る計算機109が行う演算処理の流れについて説明する。図4は、本実施形態に係る計算機109における演算処理の流れを説明するフローチャートである。図5は、RSSGシーケンスのパルスシーケンスを説明するための図である。図6は、本実施形態に係る局所位相分布算出部330における局所位相分布算出処理の流れを説明するフローチャートである。図7は、本実施形態に係る全領域局所位相分布算出部336で用いる重み画像を説明するための図である。図8は、本実施形態に係る表示部110に表示される画像の一例を示す図である。
 [計測:S1001]
 ステップS1001において、計算機109は計測処理を行う。
 操作者が、入力装置115を介して各種パラメータ(例えば、核磁気共鳴信号の数、エコー時間、エコー時間の間隔など)を設定し、計算機109に撮像開始指示が入力されると、計測制御部310は、予め設定されたパルスシーケンスに従って、シーケンス制御装置114を制御する。そして、計測制御部310は、任意のエコー時間の核磁気共鳴信号を複素信号として計測する。
 ここで、図5を参照して、計測制御部310が計測に使用するGradient Echo型のパルスシーケンスについて説明する。図5に示すRSSG(RF-spoiled-Steady-state Acquisition with Rewound Gradient-Echo)シーケンス510は、生体組織間の磁化率差によって生じる局所的な磁場変化に鋭敏なパルスシーケンスである。なお、計測制御部310が計測に使用するシーケンスは、RSSGシーケンス510に限定されるものではない。
 図5において、RFは高周波磁場を示し、Gsはスライス傾斜磁場を示し、Gpは位相エンコード傾斜磁場を示し、Grはリードアウト傾斜磁場を示している。また、TEはエコー時間を示し、TRは繰り返し時間を示している。更に、501はスライス傾斜磁場パルス、502は高周波磁場パルス、503、508はスライスエンコード傾斜磁場パルス、504、509は位相エンコード傾斜磁場パルス、505、506はリードアウト傾斜磁場パルス、507は核磁気共鳴信号を示している。なお、ハイフン以下の数字は、計測制御部310が行う繰り返しの計測が何回目であるかを示しており、例えば、「-2」は、計測制御部310が行う繰り返しの計測が2回目であることを示している。
 計測制御部310は、RSSGシーケンス510に従って、シーケンス制御装置114を制御することで、スライス傾斜磁場パルス501、高周波磁場パルス502、スライスエンコード傾斜磁場パルス503、位相エンコード傾斜磁場パルス504、リードアウト傾斜磁場パルス505、リードアウト傾斜磁場パルス506、スライスエンコード傾斜磁場パルス508、位相エンコード傾斜磁場パルス509、などの印加タイミングを適切に制御することで、画像再構成に必要なデータを計測する。
 図5に示すRSSGシーケンス510では、繰り返し時間TR内に、1つの核磁気共鳴信号507を計測する。
 まず、静磁場内に被検体101が配置されると、静磁場強度に応じて被検体101内の水素原子核が、所定の周波数を持つため、送信部105は、被検体101の計測領域に、当該周波数と一致するような高周波磁場パルス502-1を送信し、被検体101の水素原子核を励起させる。この際、傾斜磁場印加部103は、被検体101の計測領域におけるz方向に、スライス傾斜磁場パルス501-1を印加し、被検体101における所定のスライス(3次元画像の撮像位置)での磁化を選択励起させる。
 なお、計測制御部310は、繰り返しの計測を行う際、高周波磁場パルス502の位相を、所定角度(例えば、117度、122度)ずつ変更する。例えば、繰り返しの計測が1回目の場合、所定角度を117度とし、繰り返しの計測が2回目の場合、角度を117×2度とし、繰り返しの計測が3回目の場合、角度を117×3度とする。これにより、次の高周波磁場パルスを印加する際、余分な信号成分を残すことなく適切な計測を行うことが可能になる。
 続いて、傾斜磁場印加部103は、被検体101の計測領域におけるz方向に、スライスエンコード傾斜磁場パルス503-1を印加して、磁化の位相に、スライスエンコード方向の位置情報を付加する。また、傾斜磁場印加部103は、被検体101の計測領域におけるx方向に、位相エンコード傾斜磁場パルス504-1を印加して、磁化の位相に、位相エンコード方向の位置情報を付加する。即ち、傾斜磁場印加部103は、核磁気共鳴信号の読み出しに先だって、z方向及びx方向における磁化の位相に、位置情報を付加する。
 スライスエンコード傾斜磁場パルス503-1は、繰り返し時間TR毎に、磁場強度が変化する。例えば、図5に示す7本の横線は、異なる強さの傾斜磁場を示しており、上方向の矢印に従って、正の傾斜磁場が強くなることを示している。同様に、位相エンコード傾斜磁場パルス504-1も、繰り返し時間TR毎に、磁場強度が変化する。例えば、図5に示す7本の横線は、異なる強さの傾斜磁場を示しており、上方向の矢印に従って、正の傾斜磁場が強くなることを示している。
 なお、図5では、スライスエンコード傾斜磁場パルス503-1及び位相エンコード傾斜磁場パルス504-1における磁場強度が、7通りに変化する場合を一例に挙げて説明しているが、これに限定されるものではない。
 続いて、傾斜磁場印加部103は、被検体101の計測領域におけるy方向に、ディフェーズ用のリードアウト傾斜磁場パルス505-1を印加して、画素内における磁化の位相を分散させる。そして、傾斜磁場印加部103は、被検体101の計測領域におけるy方向に、リードアウト傾斜磁場パルス506-1を印加して、磁化の位相に、リードアウト方向の位置情報を付加する。これにより、計測制御部310は、繰り返し時間TR内に、1つの核磁気共鳴信号507を計測することが可能になる。
 最後に、傾斜磁場印加部103は、被検体101の計測領域におけるz方向に、リフェーズ用のスライスエンコード傾斜磁場パルス508-1を印加し、更に、被検体101の計測領域におけるx方向に、リフェーズ用の位相エンコード傾斜磁場パルス509-1を印加して、ディフェーズされた画素内における磁化の位相を収束させる。
 スライスエンコード傾斜磁場パルス508-1は、繰り返し時間TR毎に、磁場強度が変化する。例えば、図5に示す7本の横線は、異なる強さの傾斜磁場を示しており、下方向の矢印に従って、負の傾斜磁場が強くなることを示している。同様に、位相エンコード傾斜磁場パルス509-1も、繰り返し時間TR毎に、磁場強度が変化する。例えば、図5に示す7本の横線は、異なる強さの傾斜磁場を示しており、下方向の矢印に従って、負の傾斜磁場が強くなることを示している。
 計測制御部310は、上述の手順を、高周波磁場パルス502の位相、スライスエンコード傾斜磁場パルス503、508の強度、位相エンコード傾斜磁場パルス504、509の強度、などを変化させながら繰り返し行い、画像再構成に必要とされる核磁気共鳴信号507を計測する。なお、RSSGシーケンス510により得られる複素画像の絶対値成分は、エコー時間TEが短く設定されるとT1(縦緩和時間)強調画像となり、エコー時間TEが長く設定されると画素内の位相分散を反映したT2(見かけの横緩和時間)強調画像となる。本実施形態では、エコー時間TEを長く設定し、T2強調画像の撮像を行う。
 [画像再構成:S1002]
 ステップS1002において、計算機109は画像再構成処理を行う。
 画像再構成部320は、ステップS1001で計測された核磁気共鳴信号を、k空間上に配置してフーリエ変換することにより、複素画像を再構成する。撮像方法としては、特に限定されるものではなく、カーテシアン撮像、ノンカーテシアン撮像、マルチエコーエコープラナーイメージング法、などの公知の方法を適用することができる。
 [局所位相分布算出処理:S1003]
 ステップS1003において、計算機109は局所位相分布算出処理を行う。
 ここで、図6を参照して、局所位相分布算出処理の流れについて詳細に説明する。局所位相分布算出部330は、ステップS1002で再構成された複素画像に基づいて、対象組織の全領域における局所位相分布を算出する。
 ステップS1101において、対象組織抽出部331は、画像再構成部320によって構成された複素画像の絶対値成分(絶対値画像)における画素値の大きさに基づいて、被検体101の対象組織を抽出した対象組織抽出画像Mを算出する。対象組織抽出部331は、例えば、判別分析法などの公知の方法を用いて、絶対値画像の画素値のヒストグラムから、ノイズ成分と信号成分とを分離し、被検体101の対象組織を抽出した対象組織抽出画像M(例えば、脳内部の画素値を1、脳外部の画素値を0とした画像)を算出する。なお、対象組織抽出部331は、例えば、領域拡大(Region Growing)法などの公知の方法を用いて、複素画像から磁化率分布の算出に不要な領域(例えば、皮下脂肪などの領域)を除去して、磁化率分布の算出に必要な領域である対象組織の実質領域のみを抽出した対象組織抽出画像Mを算出してもよい。
 ステップS1102において、位相分布算出部332は、画像再構成部320によって再構成された複素画像の偏角成分に基づいて、位相分布Pを算出する。
 ステップS1103において、位相折り返し補正部333は、例えば、領域拡大法など公知の方法を用いて、ステップS1102で算出された位相分布Pに生じている-π~+πの範囲に折り返された位相を補正し、総位相分布Ptotalを算出する。
 ステップS1104において、領域分割部334は、ステップS1101で算出された対象組織抽出画像Mに基づいて、総位相分布Ptotalを、少なくとも2つ以上の領域、即ち、対象組織の所定領域と、対象組織の所定領域以外の領域と、に分割する。対象組織の所定領域とは、例えば、脳内部の画素値を1、脳外部の画素値を0とした画像において、画素値が1である領域と画素値が0である領域との境界部分である脳表領域である。また、対象組織の所定領域以外の領域とは、例えば、脳表領域を除いた脳内部の領域である。なお、脳表領域として定義される領域は、後述のステップS1106において、全領域局所位相分布算出部336が、対象組織の全領域における局所位相分布を算出する際に用いる球状カーネルρの半径r(カーネルサイズ)によって決定される。
 ステップS1105において、背景位相分布算出部335は、ステップS1103で算出された総位相分布Ptotalから、対象組織の所定領域(例えば、脳表領域)における背景位相分布Pedgeを算出する。以下、対象組織の所定領域における背景位相分布Pedgeの算出法について説明する。
 背景位相分布算出部335は、局所多項式近似処理により、総位相分布Ptotalから、対象組織の所定領域における背景位相分布Pedgeを算出する。
 背景位相分布算出部335は、局所領域内の任意画素の位置を(x、y、z)とし、式(6)を満たす多項式係数Ak、l、nを算出する。
Figure JPOXMLDOC01-appb-M000006
 
 
 Kはx方向の多項式における次数、Lはy方向の多項式における次数、Nはz方向の多項式における次数をそれぞれ表す。
 局所領域内の中心画素の位置を(x、y、z)とすると、対象組織の所定領域における背景位相分布Pedgeは、式(7)により算出できる。
Figure JPOXMLDOC01-appb-M000007
 
 背景位相分布算出部335は、対象組織の所定領域における画素の全てに対して、式(7)を実施することにより、背景位相分布Pedgeを算出する。
 ステップS1106において、全領域局所位相分布算出部336は、総位相分布Ptotalと、対象組織の所定領域における背景位相分布Pedgeとを組み合わせて、対象組織の全領域における局所位相分布P’localを算出する。以下、対象組織の全領域における局所位相分布P’localの算出法について説明する。
 全領域局所位相分布算出部336は、制約付き最小二乗処理により、対象組織の全領域における局所位相分布の推定値P’localを算出する。
Figure JPOXMLDOC01-appb-M000008
 
 
 Wは脳表領域における重みが小さい重み画像、Wは脳表領域における重みが大きい重み画像をそれぞれ表す。
 式(8)における第一項は、SHARP法と同様、脳内部の領域における、平均値の性質に基づく制約項を表している。平均値の性質における条件を満たす制約項に対しては、脳内部の領域における重みを大きくし、脳表領域における重みを小さくする。
 式(8)における第二項は、脳表領域における背景位相分布に基づく制約項を表している。脳表領域における背景位相分布に基づく制約項に対しては、脳表領域における重みを大きくし、脳内部の領域における重みを小さくする。
 ここで、図7を参照して、重み画像W及び重み画像Wと式(8)との関係について説明する。図7(a)は、重み画像Wの一次元プロファイルを示している。図7(b)は、重み画像Wの一次元プロファイルを示している。横軸は位置を示し、縦軸は重みを示している。
 図7において、領域701と領域702とを合わせた領域は、脳表領域を示し、領域703は、脳表領域を除いた脳内部の領域(脳表領域以外の領域)を示している。同様に、領域704と領域705とを合わせた領域は、脳表領域を示し、領域706は、脳表領域を除いた脳内部の領域(脳表領域以外の領域)を示している。なお、領域701(領域704)及び領域702(領域705)は、球状カーネルρの半径rよりも大きくなるように設定される。
 図7に示すように、全領域局所位相分布算出部336は、領域701での重みが0となるように設定し、領域703での重みが1となるように設定し、領域701から領域703へと遷移する領域702での重みが0から1の間の値(線形に変化する値)となるように設定する。また、全領域局所位相分布算出部336は、領域704での重みが1となるように設定し、領域706での重みが0となるように設定し、領域704から領域706へと遷移する領域705での重みが0から1の間の値(線形に変化する値)となるように設定する。
 即ち、全領域局所位相分布算出部336は、重み画像W及び重み画像Wを用いて、脳内部の領域における重みを大きくし、脳表領域における重みを小さくすることで、式(8)に示す平均値の性質に基づく制約項を、SHARP法と同等の精度で算出することができる(式(5)参照)。
 また、全領域局所位相分布算出部336は、重み画像W及び重み画像Wを用いて、脳表領域における重みを大きくし、脳内部の領域における重みを小さくすることで、式(8)に示す脳表領域における背景位相分布に基づく制約項を、局所多項式近似処理により算出された脳表領域における背景位相分布に近づけることができる。
 更に、全領域局所位相分布算出部336は、重み画像W及び重み画像Wを用いて、領域702及び領域706における重みを、0から1の間で線形に変化させることで、脳表領域(領域701と領域702とを合わせた領域、領域704と領域705とを合わせた領域)と、脳内部の領域(領域703、領域706)とを、滑らかに接続することができる。
 なお、全領域局所位相分布算出部336は、式(8)に示すような制約付き最小二乗処理のみならず、式(9)に示すような正則化付最小二乗処理(局所位相分布の正則化項を最小化する条件)を用いて、対象組織の全領域(脳表領域を含めた脳全体)における局所位相分布を算出してもよい。
 例えば、全領域局所位相分布算出部336は、局所位相分布Plocalの画素値のばらつきを考慮した正則化項を挿入して、式(9)を実施することにより、対象組織の全領域における局所位相分布を算出することも可能である。
Figure JPOXMLDOC01-appb-M000009
 
 λは、対象組織の所定領域における制約の大きさを調整するための正則化パラメータである。λは、局所位相分布Plocalのばらつきを、ある程度許容するための正則化パラメータである。
 式(9)における第一項は、SHARP法と同様、脳内部の領域における、平均値の性質に基づく制約項を表している。
 式(9)における第二項は、脳表領域における背景位相分布に基づく制約項を表している。
 式(9)における第三項は、位相分布から局所位相変化成分を過剰に除去してしまうことを防ぐ役割を持つ正則化項を表している。
 式(8)及び式(9)より、本実施形態に係る計算機109の演算処理によれば、脳内部の領域における局所位相分布の算出精度は、従来のSHARP法と略同様の算出精度を維持しつつ、脳表領域における局所位相分布の算出精度を、従来と比較して格段に向上させることができることが示された。即ち、本実施形態に係る計算機109の演算処理によれば、脳内部の領域における平均値の性質に基づく制約項及び脳表領域における背景位相分布に基づく制約項から、脳表領域を含めた脳全体の局所位相分布を高精度に算出できることが示された。
 [所望分布算出:S1004]
 S1004において、計算機109は所望分布算出処理を行う。
 所望分布算出部340は、ステップS1003で算出された局所位相分布Plocalを用いて、QSM法やSWI法などにより、定量的磁化率分布、磁化率強調画像などの所望分布を算出する。例えば、所望分布算出部340は、局所磁場算出部と磁化率分布算出部とを備え、QSM法により、定量的磁化率分布を算出する。また、例えば、所望分布算出部340は、マスク作成部と磁化率強調画像作成部とを備え、SWI法により、磁化率強調画像を算出する。以下、所望分布算出部340が行う具体的な処理例として、QSM法を用いた定量的磁化率分布の算出法と、SWI法を用いた磁化率強調画像の算出法について説明する。
 ≪QSM法を用いた定量的磁化率分布の算出法≫
 まず、QSM法による定量的磁化率分布の算出法について説明する。QSM法は、局所位相分布Plocalが、生体組織間の磁化率差を反映することを利用して、局所位相分布Plocalから、局所的な磁場分布を算出し、磁場と磁化率の関係式を用いて定量的磁化率分布を算出する手法である。
 位置ベクトルをrとすると、生体組織間の磁化率差によって生じる相対的な磁場分布は、式(10)で表される。
Figure JPOXMLDOC01-appb-M000010
 
 gはプロトンの核磁気回転比、Bは静磁場強度をそれぞれ表す。
 ここで、局所的な磁場分布は、静磁場に関するマクスウェル方程式より、式(11)で表される。
Figure JPOXMLDOC01-appb-M000011
 
 ここで、d(r)は局所磁場分布、c(r)は生体内の磁化率分布、aはベクトル(r’-r)と静磁場方向との為す角度、d(r)は点ダイポール磁場をそれぞれ表す。
 即ち、式(11)に示すように、局所磁場分布d(r)は、磁化率分布c(r)と点ダイポール磁場d(r)との畳み込み積分で表すことができる。従って、式(11)の両辺をフーリエ変換することにより、式(11)は、式(12)に変換できる。
Figure JPOXMLDOC01-appb-M000012
 
 k=(k、k、k)は、k空間上の位置ベクトルを表す。
 D(k)は、局所磁場分布d(r)のフーリエ成分を表す。
 D(k)は、点ダイポール磁場d(r)のフーリエ成分を表す。
 X(k)は、生体内の磁化率分布c(r)のフーリエ成分を表す。
 式(12)に示すように、生体内の磁化率分布c(r)のフーリエ成分X(k)は、局所磁場分布d(r)のフーリエ成分D(k)を、点ダイポール磁場d(r)のフーリエ成分D(k)で除算することによって算出できる。
 しかしながら、式(12)は、D(k)=0近傍の領域において、その逆数が発散してしまうため、直接的にX(k)を算出することができない。このD(k)=0近傍の領域はマジックアングル領域と呼ばれ、静磁場方向に対しておよそ54.7度の2倍の頂角を持つ逆双円錐領域となる。
 QSM法は、マジックアングル領域の存在により不良条件逆問題(ill-conditioned inverse problem)に帰着されるため、不良条件逆問題を解くための幾つかの解法が提案されている。不良条件逆問題を解く代表的な方法としては、例えば、磁場と磁化率との関係式に基づく制約条件下で、磁場分布から算出した磁化率分布に対して平滑化処理を行うことを繰り返す方法、磁場分布及び点ダイポール磁場のk空間上の演算により磁化率分布を算出するTKD(Truncated-based K-space Division)法、TKD法で算出した磁化率分布と閾値処理により微細構造を抽出した磁化率分布とを繰り返し演算により合成するIterative SWIM(Susceptibility Weighted Imaging and Mapping)法、正則化付最小二乗法を用いたMEDI(Morphology enabled dipole inversion)法、などが挙げられる。なお、所望分布算出部340は、これらの何れの方法を用いてもよい。
 ≪SWI法を用いた磁化率強調画像の算出法≫
 次に、SWI法による磁化率強調画像の算出法について説明する。SWI法は、局所位相分布Plocalから、磁化率分布を強調する磁化率強調マスクを作成し、計測した強度画像(絶対値画像)に、磁化率強調マスクを乗算することで、磁化率強調画像を算出する手法である。
 所望分布算出部340は、所定の磁化率強調マスクを作成し、この磁化率強調マスクを、所定回数乗算する。その後、所望分布算出部340は、任意のエコー時間の絶対値画像に、当該磁化率強調マスクを乗算することで、磁化率強調画像を算出する。磁化率強調画像とは、磁化率の高い所が暗くなり、磁化率の低い所が明るくなる画像である。操作者は、当該画像に基づいて、磁化率が高くなっているか、或いは低くなっているかを判断することができる。なお、所望分布算出部340は、公知のさまざまな手法を用いて、磁化率強調画像を算出することが可能であり、SWI法に限定されるものではない。
 [画像表示:S1005]
 S1005において、計算機109は画像表示処理を行う。
 計算機109は、ステップS1004で算出された定量的磁化率分布や磁化率強調画像を、表示部110のディスプレイ画面上に表示させる。また、計算機109は、ステップS1004で算出された各種データを外部記憶装置111或いは記憶装置に記憶させ、所望の表示装置に、定量的磁化率分布や磁化率強調画像を表示させることも可能である。
 ここで、図8を参照して、MRI装置100の表示部110に表示される画像の一例について説明する。
 図8(a)は、本実施形態に係る計算機109によって演算処理が施された画像(磁化率強調画像)であり、図8(b)は、図8(a)の領域801における拡大図である。図8(c)は、従来の計算機によって演算処理が施された画像(磁化率強調画像)であり、図8(d)は、図8(c)の領域802における拡大図である。
 図8(a)及び図8(b)に示す画像と図8(c)及び図8(d)に示す画像とを比較すると、図8(a)及び図8(b)に示す画像では、脳表領域が鮮明に撮像されている一方で、図8(c)及び図8(d)に示す画像では、脳表領域が太い実線で囲まれたように撮像されていることがわかる。即ち、本実施形態に係る計算機109は、従来の計算機と比較して、太い実線で囲まれた領域の分だけ、より多くの情報を算出できることがわかる。
 従って、操作者は、脳全体における生体組織の構造や組成などを診断する場合、図8(c)及び図8(d)に示す画像を用いるよりも、図8(a)及び図8(b)に示す画像を用いることで、より適切な診断を行うことができる。更に、操作者は、くも膜下出血や脳梗塞など、脳の表面に発生した病気を診断する場合、図8(a)及び図8(b)に示す画像を用いることで、正確な診断を行うことができる。即ち、本実施形態に係る計算機109が所定の演算処理を行うことで、操作者は、脳表領域において信頼性の高い画像を取得することができる。
 以上、上述の各処理を行うことで、本実施形態に係る計算機109が行う演算処理が終了する。なお、上述の図4及び図6に示すフローは一例であり、その一部を省略することや、別の処理を追加することが可能であることは勿論である。
 本実施形態に係るMRI装置によれば、計算機109が、複素画像から、例えば脳領域を抽出し、総位相分布を脳表領域と脳内部の領域とに分割して、脳表領域における背景位相分布を算出し、総位相分布と該背景位相分布とを組み合わせることで、脳全体の局所位相分布を算出する。これにより、脳表領域を含めた脳全体の局所位相分布を高精度に算出するMRI装置を実現することができる。
 更に、本実施形態に係るMRI装置によれば、高精度な局所位相分布を用いて、定量的磁化率分布や磁化率強調画像を算出することができる。これにより、操作者は、鮮明な領域が拡張された画像を取得することができるため、MRI装置の信頼性を高められる。
 100,120,130 磁気共鳴イメージング装置
 101         被検体
 102         静磁場生成磁石
 103         傾斜磁場印加部
 105         送信部
 106         受信部
 109         計算機
 110         表示部
 331         対象組織抽出部
 332         位相分布算出部
 333         位相折り返し補正部
 334         領域分割部
 335         背景位相分布算出部
 336         全領域局所位相分布算出部(局所位相分布算出部)

Claims (6)

  1.  被検体が配置される空間に静磁場を生成する静磁場生成磁石と、
     前記被検体に高周波磁場を送信する送信部と、
     前記高周波磁場の送信により前記被検体から発生する核磁気共鳴信号を受信する受信部と、
     前記核磁気共鳴信号に位置情報を付加する傾斜磁場を印加する傾斜磁場印加部と、
     前記送信部、前記受信部、及び前記傾斜磁場印加部を制御するとともに、前記核磁気共鳴信号に対して演算処理を行う計算機と、
     前記計算機により演算処理された画像を表示する表示部と、を備え、
     前記計算機は、
     前記核磁気共鳴信号によって構成された複素画像から前記被検体の対象組織を抽出した対象組織抽出画像を算出する対象組織抽出部と、
     前記複素画像から位相分布を算出する位相分布算出部と、
     前記位相分布の位相折り返しを補正する位相折り返し補正部と、
     前記対象組織抽出画像に基づいて、補正後の位相分布を、前記対象組織の所定領域と前記所定領域以外の領域とに分割する領域分割部と、
     前記所定領域の背景位相分布を算出する背景位相分布算出部と、
     前記補正後の位相分布と前記背景位相分布とを組み合わせて、前記対象組織の所定領域における局所位相分布を算出する局所位相分布算出部と、
     を備えることを特徴とする磁気共鳴イメージング装置。
  2.  前記局所位相分布算出部は、
     任意画素の画素値と、前記任意画素を中心とした所定の半径を有する球内における画素の画素値の平均とが略等しいという条件に基づいて、前記局所位相分布を算出する、
     ことを特徴とする請求項1に記載の磁気共鳴イメージング装置。
  3.  前記局所位相分布算出部は、
     前記局所位相分布の正則化項を最小化する条件に基づいて、前記局所位相分布を算出する、
     ことを特徴とする請求項2に記載の磁気共鳴イメージング装置。
  4.  前記計算機は、
     前記局所位相分布を用いて、局所磁場分布を算出し、
     前記局所磁場分布及び磁場と磁化率との関係式を用いて、磁化率分布を算出する、
     ことを特徴とする請求項1から請求項3のいずれか一項に記載の磁気共鳴イメージング装置。
  5.  前記計算機は、
     前記局所位相分布を用いて、前記磁化率分布を強調する磁化率強調マスクを作成し、
     前記複素画像の絶対値成分に、前記磁化率強調マスクを乗じて磁化率強調画像を算出する、
     ことを特徴とする請求項4に記載の磁気共鳴イメージング装置。
  6.  前記領域分割部は、
     前記補正後の位相分布を、脳表領域と脳内部の領域とに分割する、
     ことを特徴とする請求項1から請求項5のいずれか一項に記載の磁気共鳴イメージング装置。
PCT/JP2018/021266 2017-09-11 2018-06-01 磁気共鳴イメージング装置 WO2019049443A1 (ja)

Priority Applications (1)

Application Number Priority Date Filing Date Title
US16/646,033 US11051695B2 (en) 2017-09-11 2018-06-01 Magnetic resonance imaging apparatus

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
JP2017174412A JP6987339B2 (ja) 2017-09-11 2017-09-11 磁気共鳴イメージング装置
JP2017-174412 2017-09-11

Publications (1)

Publication Number Publication Date
WO2019049443A1 true WO2019049443A1 (ja) 2019-03-14

Family

ID=65633633

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2018/021266 WO2019049443A1 (ja) 2017-09-11 2018-06-01 磁気共鳴イメージング装置

Country Status (3)

Country Link
US (1) US11051695B2 (ja)
JP (2) JP6987339B2 (ja)
WO (1) WO2019049443A1 (ja)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112515653B (zh) * 2020-10-09 2024-03-26 天津大学 一种基于核磁共振图像的脑网络构建方法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20130221961A1 (en) * 2012-02-27 2013-08-29 Medimagemetric LLC System, Process and Computer-Accessible Medium For Providing Quantitative Susceptibility Mapping
WO2014076808A1 (ja) * 2012-11-16 2014-05-22 株式会社日立製作所 磁気共鳴イメージング装置および定量的磁化率マッピング法
JP2017184935A (ja) * 2016-04-04 2017-10-12 株式会社日立製作所 磁気共鳴イメージング装置、及び、画像処理方法

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5539316A (en) * 1995-08-25 1996-07-23 Bruker Instruments, Inc. Shimming method for NMR magnet having large magnetic field inhomogeneities
US9448289B2 (en) 2010-11-23 2016-09-20 Cornell University Background field removal method for MRI using projection onto dipole fields
JP6640530B2 (ja) * 2015-10-30 2020-02-05 キヤノンメディカルシステムズ株式会社 磁気共鳴イメージング装置、医用画像処理装置及び画像処理方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20130221961A1 (en) * 2012-02-27 2013-08-29 Medimagemetric LLC System, Process and Computer-Accessible Medium For Providing Quantitative Susceptibility Mapping
WO2014076808A1 (ja) * 2012-11-16 2014-05-22 株式会社日立製作所 磁気共鳴イメージング装置および定量的磁化率マッピング法
JP2017184935A (ja) * 2016-04-04 2017-10-12 株式会社日立製作所 磁気共鳴イメージング装置、及び、画像処理方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
FANG, JINSHENG ET AL.: "Background field removal using a region adaptive kernel for quantitative susceptibility mapping of human brain", J. MAGN. RESON., vol. 281, 10 May 2017 (2017-05-10), pages 130 - 140, XP055582441 *
KIM, HAHNSUNG ET AL.: "Robust Susceptibility Weighted Imaging using Single-Slab 3D GRASE with Removal of Background Phase Variation", PROC. INTL. SOC. MAG. RESON. MED., vol. 21, 26 April 2013 (2013-04-26) *

Also Published As

Publication number Publication date
JP2022002762A (ja) 2022-01-11
JP6987339B2 (ja) 2021-12-22
US20200260958A1 (en) 2020-08-20
JP7230149B2 (ja) 2023-02-28
JP2019047978A (ja) 2019-03-28
US11051695B2 (en) 2021-07-06

Similar Documents

Publication Publication Date Title
CN107072592B (zh) 磁共振成像装置以及定量性磁化率匹配方法
JP5843876B2 (ja) 磁気共鳴イメージング装置および磁化率強調画像生成方法
US9766316B2 (en) Magnetic resonance imaging device and quantitative susceptibility mapping method
CN108778116B (zh) 磁共振成像装置以及图像处理方法
JP5982494B2 (ja) 磁気共鳴イメージング装置
JP6679467B2 (ja) 磁気共鳴イメージング装置および酸素摂取率算出方法
US10552955B2 (en) Imaging acceleration methods for MRI parameter mapping
US9709641B2 (en) Magnetic resonance imaging apparatus, image processing apparatus, and susceptibility map calculation method
US11796617B2 (en) System and method for reconstruction of magnetic resonance images acquired with partial Fourier acquisition
US20190076049A1 (en) Magnetic resonance imaging device, magnetic resonance imaging method and susceptibility calculation program
JP7230149B2 (ja) 磁気共鳴イメージングで得た画像の処理方法、画像処理プロブラム、及び、計算機
JP5650724B2 (ja) 磁気共鳴イメージング装置
US10818047B2 (en) Iterative reconstruction of quantitative magnetic resonance images

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 18853252

Country of ref document: EP

Kind code of ref document: A1

NENP Non-entry into the national phase

Ref country code: DE

122 Ep: pct application non-entry in european phase

Ref document number: 18853252

Country of ref document: EP

Kind code of ref document: A1