WO2017204251A1 - 磁気共鳴イメージング装置、磁気共鳴イメージング方法及び磁化率算出プログラム - Google Patents

磁気共鳴イメージング装置、磁気共鳴イメージング方法及び磁化率算出プログラム Download PDF

Info

Publication number
WO2017204251A1
WO2017204251A1 PCT/JP2017/019381 JP2017019381W WO2017204251A1 WO 2017204251 A1 WO2017204251 A1 WO 2017204251A1 JP 2017019381 W JP2017019381 W JP 2017019381W WO 2017204251 A1 WO2017204251 A1 WO 2017204251A1
Authority
WO
WIPO (PCT)
Prior art keywords
image
region
mask
susceptibility
magnetic
Prior art date
Application number
PCT/JP2017/019381
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/097,284 priority Critical patent/US11141078B2/en
Publication of WO2017204251A1 publication Critical patent/WO2017204251A1/ja

Links

Images

Classifications

    • 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/24Arrangements or instruments for measuring magnetic variables involving magnetic resonance for measuring direction or magnitude of magnetic fields or magnetic flux
    • G01R33/243Spatial mapping of the polarizing magnetic field
    • 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/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
    • A61B2576/00Medical imaging apparatus involving image processing or analysis
    • A61B2576/02Medical imaging apparatus involving image processing or analysis specially adapted for a particular organ or body part
    • A61B2576/026Medical imaging apparatus involving image processing or analysis specially adapted for a particular organ or body part for the brain
    • 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
    • G01R33/565Correction of image distortions, e.g. due to magnetic field inhomogeneities
    • G01R33/56536Correction of image distortions, e.g. due to magnetic field inhomogeneities due to magnetic susceptibility variations
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/10Image acquisition modality
    • G06T2207/10072Tomographic images
    • G06T2207/10088Magnetic resonance imaging [MRI]
    • GPHYSICS
    • G16INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16HHEALTHCARE INFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR THE HANDLING OR PROCESSING OF MEDICAL OR HEALTHCARE DATA
    • G16H30/00ICT specially adapted for the handling or processing of medical images
    • G16H30/40ICT specially adapted for the handling or processing of medical images for processing medical images, e.g. editing

Definitions

  • the present invention relates to a magnetic resonance imaging (MRI) technique.
  • the present invention relates to an image processing technique for a reconstructed image.
  • MRI magnetic resonance imaging
  • a magnetic resonance imaging apparatus (hereinafter also referred to as an MRI apparatus) applies a high-frequency magnetic field and a gradient magnetic field to a subject placed in a static magnetic field, and measures and images a signal generated from the subject by nuclear magnetic resonance.
  • This is a medical image diagnostic apparatus.
  • MRI apparatus in addition to an image reflecting the density of protons (hydrogen nuclei) and relaxation times (T1, T2) (an image having absolute values as pixel values), static magnetic field inhomogeneity, magnetic susceptibility difference between living tissues, and the like An image reflecting the change in the magnetic field caused by the above (image having phase as pixel value: phase image) is obtained.
  • Non-Patent Documents 1 and 2 a quantitative susceptibility mapping (QSM) method has been proposed in which a phase image reflects a magnetic susceptibility difference between tissues, and a susceptibility distribution in a living body is estimated from the phase image.
  • a magnetic susceptibility image From the magnetic susceptibility distribution estimated by the QSM method (hereinafter referred to as a magnetic susceptibility image), the iron concentration in the basal ganglia and the oxygen uptake rate in the cerebral vein can be calculated. Therefore, the QSM method is expected to be useful for early diagnosis of neurodegenerative diseases and cerebrovascular diseases that are brain diseases.
  • application of the QSM method to the trunk is also being studied. Typical examples include detection of breast calcification and estimation of liver iron concentration (Non-patent Document 3). In the prostate, application to hypoxia region and calcification detection is expected.
  • the prostate is a male-specific organ located in the lower abdomen and is adjacent to subcutaneous fat such as the buttocks. Due to the presence of this subcutaneous fat, the calculation accuracy of the prostate susceptibility image may be reduced.
  • One factor that reduces the calculation accuracy in the QSM method is a lack of MRI signal in the peripheral region of the target tissue.
  • Subcutaneous fat in the vicinity of the prostate is prone to a reduction in accuracy due to signal loss because there are signal-free areas such as outside the body and bones and low-signal areas in the vicinity.
  • the calculation accuracy fall of a magnetic susceptibility propagates also to an adjacent area
  • Non-Patent Document 1 discloses a method of separately calculating and integrating the magnetic susceptibility of the bleeding region and other regions in the brain.
  • Non-Patent Document 2 discloses a technique for calculating the magnetic susceptibility of the liver by strongly smoothing the fat region. In any method, the magnetic susceptibility is calculated under the constraint of smoothing the pixel value of the magnetic susceptibility image. However, although noise and artifacts can be reduced by smoothing, it is considered difficult to reduce accuracy reduction due to signal loss.
  • Non-Patent Document 3 describes a method of calculating the magnetic susceptibility in the brain under the constraint that the magnetic susceptibility in the region outside the brain is zero after performing background magnetic field removal considering the region outside the brain as the background. It is disclosed.
  • a susceptibility image can be calculated with high accuracy by constraining the calculated susceptibility in a region outside the brain.
  • this method is applied to the calculation of subcutaneous fat, the magnetic field is removed considering the water region, which is a high-signal region, as the background, and the magnetic susceptibility cannot be calculated using the magnetic field information of the water region, resulting in a decrease in accuracy. To do.
  • the present invention has been made in view of the above circumstances, and an object thereof is to calculate a magnetic susceptibility with high accuracy by reducing a decrease in accuracy due to signal loss.
  • the present invention measures a nuclear magnetic resonance signal (echo) reflecting spatial magnetic field inhomogeneity, calculates a complex image from the echo, and calculates a magnetic field image. Also, a plurality of masks (for example, first masks) divided by a predetermined threshold value of the discrimination image, including three or more masks from the absolute value image calculated from the complex image, a low signal region mask representing the low signal region, and a high signal region. 1 high signal area mask, second high signal area mask). Then, using a plurality of masks including a low signal region mask, the magnetic susceptibility is calculated for each of a plurality of regions that are high signal regions and have different magnetic susceptibility. At that time, for the magnetic susceptibility of one region, a magnetic susceptibility image is calculated under a restriction (specific value restriction) in which the magnetic susceptibility of another region is set to a predetermined specific value.
  • a restriction specific value restriction
  • the MRI apparatus of the present invention is generated from a static magnetic field magnet, a magnetic field generator that applies a gradient magnetic field and a high frequency magnetic field to a space formed by the static magnetic field magnet, and a subject placed in the space.
  • a reception unit that receives a nuclear magnetic resonance signal; and a computer that performs an operation on the nuclear magnetic resonance signal, the computer controlling operations of the magnetic field generation unit and the reception unit, and effects of magnetic field inhomogeneity
  • a measurement unit that measures a nuclear magnetic resonance signal, an image reconstruction unit that performs image reconstruction using the nuclear magnetic resonance signal and creates a complex image, and a magnetic susceptibility of the subject using the complex image
  • the image conversion unit uses a magnetic field calculation unit that calculates a magnetic field reflection image from the complex image and an image created by the image reconstruction unit, and a low signal region mask corresponding to a low signal region and a high signal region Using a mask calculation unit for calculating a plurality of region masks corresponding to a plurality of regions, a magnetic field reflected image calculated by the magnetic field calculation unit, the low signal region mask and the plurality of region masks calculated by the mask calculation unit A magnetic susceptibility calculation unit that calculates a magnetic susceptibility for each of the plurality of regions, the magnetic susceptibility calculation unit including a constraint with the low signal region as a background, and a specific region of the plurality of regions The magnetic susceptibility of a region different from the specific region is calculated under the restriction that the specific magnetic susceptibility is a specific value.
  • the magnetic field reflected image is an image having a frequency or a magnetic field as a pixel value and reflecting spatial magnetic field inhomogeneity.
  • FIG. 1A is an external view of a vertical magnetic field type magnetic resonance imaging apparatus
  • FIG. 1B is an external view of a horizontal magnetic field type magnetic resonance imaging apparatus
  • FIG. 2C is an external view of the magnetic resonance imaging apparatus with a sense of openness.
  • It is a block diagram which shows schematic structure of the MRI apparatus in 1st embodiment. It is a functional block diagram of a computer in a first embodiment. It is a flowchart of the imaging process in 1st embodiment. It is explanatory drawing for demonstrating an example of the gradient echo type
  • FIG. 1 It is a figure which shows the relationship between the image used for mask calculation and magnetic susceptibility calculation in 1st embodiment, and a mask. It is a flowchart of the water fat separation process in 1st embodiment. It is a flowchart of the mask calculation process in 1st embodiment. It is a flowchart of the magnetic susceptibility calculation process in 1st embodiment. It is a functional block diagram of the image conversion part in the computer of 2nd embodiment. It is a flowchart of the image conversion process of 2nd embodiment. It is a flowchart of the mask calculation process in 2nd embodiment. It is a functional block diagram of the image conversion part in the computer of 3rd embodiment. It is a flowchart of the image conversion process of 3rd embodiment.
  • FIG. 1 is an external view of an MRI apparatus.
  • FIG. 1 (a) shows a hamburger type (open type) vertical magnetic field type MRI apparatus (vertical magnetic field MRI apparatus) 100 in which magnets are separated into upper and lower parts in order to enhance the feeling of opening.
  • FIG. 1B shows a horizontal magnetic field type MRI apparatus (horizontal magnetic field MRI apparatus) 101 using a tunnel magnet that generates a static magnetic field with a solenoid coil.
  • FIG. 1C shows an MRI apparatus 102 that uses the same tunnel-type magnet as in FIG. 1B and has a feeling of openness by shortening the depth of the magnet and tilting it obliquely.
  • 1 is an example of the vertical magnetic field method and the horizontal magnetic field method, respectively, and the present invention is not limited to these.
  • the horizontal magnetic field MRI apparatus 101 will be described as an example.
  • FIG. 2 is a block diagram showing a schematic configuration of the MRI apparatus 101.
  • the MRI apparatus 101 includes a magnet 201 that generates a static magnetic field in a direction parallel to a subject, a gradient magnetic field coil 202 that generates a gradient magnetic field, a sequencer 204, a gradient magnetic field power source 205, a high-frequency magnetic field generator 206, A probe 207 that irradiates a high-frequency magnetic field and detects a nuclear magnetic resonance signal (echo), a receiver 208, a calculator 209, a display device 210, and a storage device 211 are provided.
  • a subject (for example, a living body) 203 is placed on a bed (table) or the like and placed in a static magnetic field space generated by the magnet 201.
  • the sequencer 204 sends commands to the gradient magnetic field power source 205 and the high frequency magnetic field generator 206 to generate a gradient magnetic field and a high frequency magnetic field, respectively.
  • the generated high frequency magnetic field is applied to the subject 203 through the probe 207.
  • the echo generated from the subject 203 is received by the probe 207 and detected by the receiver 208.
  • the nuclear magnetic resonance frequency (detection reference frequency f0) as a reference for detection is set by the sequencer 204.
  • the detected signal is sent to a computer 209 where signal processing such as image reconstruction is performed.
  • the result is displayed on the display device 210.
  • the detected signal, measurement conditions, image information after signal processing, and the like may be stored in the storage device 211.
  • the sequencer 204 performs control so that each unit operates at a timing and intensity programmed in advance.
  • a program that particularly describes a high-frequency magnetic field, a gradient magnetic field, and timing and intensity of signal reception is called a pulse sequence.
  • the purpose is to calculate the magnetic susceptibility based on the spatial magnetic field change obtained from the phase image.
  • a pulse sequence that acquires at least one echo reflecting the non-uniformity of the distribution, that is, an echo whose phase of each spin is shifted due to a spatial magnetic field change is used.
  • a GrE (Gradient Echo) pulse sequence can be used.
  • the GrE pulse sequence includes, for example, an RSSG (RF-soiled-Steady-state Acquisition with Rebound Gradient-Echo) sequence.
  • the computer 209 of the present embodiment instructs the sequencer 204 to measure echo according to the measurement parameter and pulse sequence set (or preset) via the input device 212, measures the echo, and the obtained echo Are arranged in the k space, the susceptibility image is calculated by performing an operation on the echo arranged in the k space, and the obtained image is displayed on the display device 210. Further, if necessary, an ROI is set on the image, and a statistical value of a pixel in the ROI is calculated.
  • Each of these functions of the computer 209 is realized when the CPU of the computer 209 loads the program stored in the storage device 211 to the memory and executes it. Note that some or all of the functions of the computer 209 may be realized by hardware such as ASIC (Application Specific Integrated Circuit) or FPGA (Field Programmable Gate Array).
  • ASIC Application Specific Integrated Circuit
  • FPGA Field Programmable Gate Array
  • the program of the present embodiment is a measurement program that executes a predetermined pulse sequence, an image reconstruction program that creates a complex image or an absolute value image using the measured echo, and three or more using an absolute value image Calculation program for calculating different area masks, magnetic field calculation program for calculating frequency image or magnetic field image from complex image, susceptibility calculation for calculating magnetic susceptibility image using frequency image or magnetic field image and three or more different area masks Contains the program.
  • the magnetic susceptibility distribution is calculated and integrated for each of the water region (region where the water proton signal is dominant) and the fat region (region where the fat proton signal is dominant).
  • a restriction for giving the magnetic susceptibility a specific value for a specific region is given.
  • the specific area is, for example, a water area, and the specific value is zero.
  • FIG. 3 shows a functional block diagram of the computer 209 of the present embodiment that performs the above processing.
  • the computer 209 includes a measurement unit 300, an image reconstruction unit 400, an image conversion unit 500, and a display processing unit 600 as main elements.
  • the image conversion unit 500 includes a water fat separation processing unit 510, a magnetic field calculation unit 520, a mask calculation unit 530, and a magnetic susceptibility calculation unit 540. Details of the mask calculation unit 530 and the magnetic susceptibility calculation unit 540 will be described later.
  • the measurement unit 300 When various measurement parameters are set and an instruction to start imaging is received, the measurement unit 300 performs measurement (S1101). Here, measurement unit 300 instructs sequencer 204 according to a predetermined pulse sequence, acquires an echo signal, and arranges it in k-space. The sequencer 204 sends instructions to the gradient magnetic field power supply 205 and the high frequency magnetic field generator 206 in accordance with the instructions to generate a gradient magnetic field and a high frequency magnetic field, respectively. The echo received by the probe 207 and detected by the receiver 208 is received as a complex signal.
  • a pulse sequence for measuring at least one echo in which a phase shift for each spin occurs corresponding to the spatial magnetic field inhomogeneity for example, a GrE pulse sequence is used.
  • An example of a measurement sequence executed by the measurement unit 300 in the measurement of step S1101 is shown in FIG.
  • RF represents an RF pulse
  • Gs represents a slice selection gradient magnetic field
  • Gp represents a phase encoding gradient magnetic field
  • Gr represents a readout gradient magnetic field.
  • the echo indicates the acquisition timing of the echo signal.
  • echo signals are measured in the following procedure within one repetition time TR.
  • a water image and a fat image are calculated using a frequency difference between water and fat, two or more echoes are acquired within one TR.
  • echo signals are acquired at four different echo times. It is assumed that the first echo time is t 1 and the subsequent echo time interval (echo interval) is ⁇ t.
  • the RF pulse 711 is irradiated to excite the hydrogen nuclear spin of the subject 203.
  • a slice selection gradient magnetic field (Gs) 712 is applied simultaneously with the RF pulse 711 in order to select a specific slice of the subject 203.
  • a phase encoding gradient magnetic field (Gp) 713 for phase encoding is applied to the echo signal.
  • a read gradient magnetic field (Gr) 721 is applied after time t 1 from the first RF pulse 711 irradiation, and an echo signal (first echo signal) 731 is measured.
  • the readout gradient magnetic field (Gr) 722 with reversed polarity is applied to measure the echo signal (second echo signal) 732.
  • the readout gradient magnetic field (Gr) 723 with reversed polarity is applied to measure the echo signal (third echo signal) 733.
  • the third echo signal time t 4 measured from after the time ⁇ t of 733, to measure the echo signal (fourth echo signal) 734 by applying a readout gradient (Gr) 724 that polarity reversal.
  • At least one of the echo times t 1 , t 2 , t 3 , t 4 when measuring the first echo signal 731, the second echo signal 732, the third echo signal 733, and the fourth echo signal 734 is water and fat.
  • the echo time t 1 and the echo interval ⁇ t are set so that the time is such that the phase difference between the two does not become zero.
  • t an In is m / f wf. Note that m is an integer.
  • an echo time t 1 , t 2 , t 3 , t 4 , or an echo time t 1 and an echo interval ⁇ t satisfying the above conditions are selected.
  • the echo time, the echo interval, and the number of echo acquisitions are set by the user via the input device 212. Alternatively, it is set in advance.
  • the measurement unit 300 irradiates the measurement sequence 710 with irradiation of the RF pulse 711 to a predetermined imaging region of the subject 203 and echo signals 731, 732 from the region while changing the intensity of the phase encoding gradient magnetic field 713. Measurements 733 and 734 are repeated a predetermined number of times. The number of repetitions is, for example, 128 times, 256 times, or the like. As a result, the number of echo signals necessary for image reconstruction of the imaging region is repeatedly acquired.
  • first original image is formed by the first echo signal 731 for the number of repetitions, and each of the second echo signal 732, the third echo signal 733, and the fourth echo signal 734 for the number of repetitions, A second original image, a third original image, and a fourth original image are formed. These are stored and used in a storage device or the like as an original image for calculation for calculating a water image and a fat image described later.
  • the number of different echo times that is, the number of original images is not limited to four and is arbitrary.
  • non-Cartesian imaging such as a radial scan that acquires data in a rotational manner in the k space may be used.
  • the image reconstruction unit 400 performs image reconstruction processing for reconstructing an image from the echo signals of the measured echo times t 1 , t 2 , t 3 , and t 4 (step S1102).
  • each echo signal is arranged in the k space and Fourier transformed.
  • the original image I 1 first original image
  • the original image I 2 second original image
  • the original image I 3 third original image
  • the original image I 4 fourth original image
  • each original image is a complex image in which each pixel value is a complex number.
  • the image conversion unit 500 performs various image conversion processes to be described later on the obtained complex image (step S1103).
  • the image conversion unit 500 calculates the complex image obtained by the image reconstruction unit 400 into a water fat separation processing unit 510 and a water fat separation processing unit 510 that convert the complex image into a water image, a fat image, and a frequency image.
  • a mask calculation unit 530 that calculates three masks: a low signal region mask, a fat mask, and a water mask, from the magnetic field image and the three masks It has the magnetic susceptibility calculation part 540 which calculates a magnetic susceptibility image. Details of the image conversion processing of this embodiment will be described later.
  • the display processing unit 600 displays the obtained magnetic susceptibility image on the display device 210 as a grayscale image (step S1104).
  • the magnetic susceptibility image may be displayed by integrating a plurality of spatially continuous image information using a method such as a maximum value projection process or a minimum value projection process.
  • image processing may be performed on the magnetic susceptibility image to create an image having a contrast different from that of the magnetic susceptibility image and displayed on the display device 210.
  • an enhancement mask that emphasizes the magnetic susceptibility difference between tissues may be created from the magnetic susceptibility image, and a magnetic susceptibility difference enhanced image obtained by multiplying it by an absolute value image may be displayed.
  • a background magnetic field removal process is performed to calculate a local magnetic field due to a magnetic susceptibility difference between living tissues, excluding a global magnetic field change due to a magnetic susceptibility difference inside and outside the body. . Thereafter, the magnetic susceptibility is calculated based on the relational expression between the magnetic field change and the magnetic susceptibility distribution.
  • the background magnetic field removal method there are a SHARP (sophisticated harmonic artifact reduction for phase data) method and a PDF (projection on dipole field) method.
  • SHARP statistical harmonic artifact reduction for phase data
  • PDF projection on dipole field
  • the magnetic susceptibility is calculated using the relational expression (formula (1)) between the relative magnetic field change ⁇ and the magnetic susceptibility distribution ⁇ expressed by the following formula.
  • represents an angle formed by the vector (r′ ⁇ r) and the static magnetic field direction
  • d (r) represents a point dipole magnetic field
  • r represents a position vector corresponding to the coordinates of the pixel.
  • the magnetic field distribution ⁇ (r) is represented by a convolution integral of the magnetic susceptibility distribution ⁇ (r) and the point dipole magnetic field d (r). Therefore, the formula (1) is converted into the following formula (2) by Fourier-transforming both sides of the formula (1).
  • k (k x , k y , k z ) is a position vector in k space, ⁇ (k), X (k), D (k), magnetic field distribution ⁇ (r), magnetic susceptibility
  • the distribution ⁇ (r) and the Fourier component of the point dipole magnetic field d (r) are respectively represented.
  • the Fourier component X (k) of the magnetic susceptibility distribution can be calculated by dividing the Fourier component ⁇ (k) of the magnetic field distribution by the Fourier component D (k) of the point dipole magnetic field.
  • X (k) cannot be directly calculated.
  • the QSM method for estimating the magnetic susceptibility distribution from the magnetic field distribution due to the presence of the magic angle has been reduced to an ill-conditioned inverse problem, and several solutions have been proposed.
  • an error function e ( ⁇ ) shown in the following equation (3) is used to obtain a magnetic susceptibility image that minimizes this.
  • is a column vector of a magnetic field image
  • is a column vector of a magnetic susceptibility image candidate
  • C is a matrix corresponding to a convolution operation on ⁇
  • W is a weight
  • is an arbitrary constant
  • G is a differential operator.
  • the second term is a regularization term, which represents a constraint for smoothing the magnetic susceptibility value and reduces noise and artifacts.
  • the above is a general magnetic susceptibility calculation method.
  • the calculation accuracy may decrease in the water region around the subcutaneous fat.
  • the decrease in calculation accuracy occurs as uneven magnetic susceptibility values in the water region around the subcutaneous fat. For example, when the magnetic susceptibility is calculated under the restriction of smoothing as in equation (3), noise and artifacts are reduced, but it is difficult to reduce the accuracy of the fat region and the unevenness of the water region due to signal loss. is there.
  • the susceptibility of the fat region and the water region is separately calculated so that the calculated susceptibility decrease of fat does not propagate to the surroundings.
  • the fat content image is calculated by the water fat separation process, the region is divided into the fat region and the water region using the fat content image, and the magnetic susceptibility of the fat region and the water region is calculated separately. Finally, integrate them.
  • the background magnetic field is removed with the air region and the fat region as the background, and the magnetic susceptibility is calculated.
  • the magnetic susceptibility is calculated under the restriction that the magnetic region does not have uneven magnetic susceptibility.
  • FIG. 6 is a diagram showing a processing flow
  • FIG. 7 is a diagram showing a relationship between a created image and a mask.
  • the image conversion process in the present embodiment includes four steps of a water / fat separation process S1201, a magnetic field calculation process S1203, a mask image calculation process S1203, and a magnetic susceptibility image calculation process S1204.
  • the water / fat separation processing unit 510 of the present embodiment is obtained by measurement with a signal model represented by water, fat, R 2 * (apparent transverse magnetization relaxation rate common to water and fat), and offset frequency distribution.
  • a water image, a fat image, and a frequency image are calculated by the water / fat separation processing (processing within a dotted square in FIG. 7).
  • the offset frequency is an offset amount of the resonance frequency that varies spatially due to non-uniformity of the static magnetic field
  • the frequency image is an image having the offset frequency as a pixel value.
  • the “water image” and “fat image” mean an image in which the signal from the water proton is dominant and an image in which the signal from the fat proton is dominant, respectively.
  • FIG. 8 is a flowchart of the water-fat separation process of the present embodiment.
  • a signal model representing a change in signal value due to echo time is set (step S1301).
  • various initial values for fitting processing to be described later are set (step S1302).
  • the fitting process of fitting the measured signal to the signal model by non-linear least square method Perform (step S1303) to calculate a water / fat / frequency image.
  • a fat content image is calculated from the obtained water image and fat image (step S1304). Details of each process will be described below.
  • a signal model is set in which the water signal intensity, the fat signal intensity, the offset frequency, and the apparent relaxation rate R 2 * are fitting variables.
  • t n is the n-th echo time
  • is the offset frequency
  • ⁇ w and ⁇ f are the complex signal intensities of water and fat
  • K n is the phase change amount (complex number) of the fat signal at time t n .
  • R 2 * represents the apparent transverse magnetization relaxation rate common to water and fat, respectively.
  • the signal model in step S1301 is not limited to the form of equation (4). For example, when the number of echoes obtained is small, a signal model using only the water signal intensity, fat signal intensity, and offset frequency as a fitting variable may be set.
  • the fat signal is known to have a plurality of spectral peaks due to its molecular structure. Therefore, if (the P to an integer of 1 or more) P number considered fat signals having a peak of the phase change amount K n fat signal is represented by the following formula (5).
  • a p and f p represent the relative signal intensity of the p-th fat peak (p is an integer satisfying 1 ⁇ p ⁇ P) and the frequency difference with water, respectively. Note that ap satisfies the condition of the following formula (6).
  • the initial values to be set are the complex signal strength of water and the fat complex signal strength of each pixel, the offset frequency distribution, and the apparent transverse magnetization relaxation rate R 2 * .
  • the initial value of each pixel of the complex signal strength ⁇ w of water and the complex signal strength ⁇ f of fat is the maximum value in the time direction of the absolute value
  • max is used.
  • the initial value of the offset frequency distribution is 0 for all pixels.
  • the initial value of the apparent transverse magnetization relaxation rate R 2 * is 1 for all pixels.
  • the initial value of each pixel of water and fat of the complex signal strength [rho w and [rho f, the initial value of the offset frequency distribution, transverse magnetization relaxation rate R 2 * initial value of the apparent may not necessarily be the above value
  • the calculation result by the non-linear least square method may be a value that does not diverge or vibrate.
  • the variables for estimating the true value by fitting are the water signal strength ⁇ w , the fat signal strength ⁇ f , the apparent transverse magnetization relaxation rate R 2 * , and the offset frequency ⁇ of each pixel.
  • the estimated values of the variables are ⁇ w ′, ⁇ f ′, R 2 * ′, and ⁇ ′, and the differences between the true values and the estimated values are ⁇ w , ⁇ f , ⁇ R 2 * , and ⁇ , respectively.
  • the signal value obtained by actual measurement was obtained by substituting s n , estimated values ⁇ w ′, ⁇ f ′, R 2 * ′, and ⁇ ′ into the signal model represented by Expression (4).
  • the estimated signal is s ′ n
  • the difference ⁇ s n between the measurement signal s n and the estimated signal s ′ n is expressed by the following expression (7) in matrix notation.
  • the vector x can be calculated by the following formula (8) using the pseudo inverse matrix of the matrix A.
  • a H represents a complex conjugate transpose matrix of A.
  • ⁇ w , ⁇ f , ⁇ R 2 * , ⁇ which are the elements of the difference vector x calculated by Expression (8), are estimated values ⁇ w ′, ⁇ f ′, R 2 *.
  • the difference vector x is again calculated using Equation (8). By repeating this procedure, the difference vector x is minimized and the estimated value is brought closer to the true value.
  • an arbitrary threshold value is set, and the above procedure is repeated until the norm of the difference vector x becomes equal to or less than the threshold value.
  • the finally obtained water signal strength ⁇ w 'fat signal strength ⁇ f ', apparent transverse magnetization relaxation rate R 2 * ', and offset frequency ⁇ ' of each pixel are respectively obtained as a water image and a fat.
  • An image, an R 2 * image, and a frequency image are used.
  • any known method such as the Levenberg-Marquardt method can be used for the nonlinear least square method.
  • a fat content image is calculated from the water image and fat image obtained in the fitting process S1303 (step S1304).
  • the fat content rate image is an image in which the pixel value of each pixel is the fat content rate.
  • Each pixel value of the fat content rate in the present embodiment is calculated by dividing the pixel value of the fat image by the sum of the pixel values of the water image and the fat image.
  • the fat content rate image is used as a determination image when calculating a fat mask in a mask calculation process described later.
  • each pixel value ⁇ of the frequency image is obtained by the following equation (9) using each pixel value I 1 (r) of the original image I 1 and the echo time t 1. 'Can be calculated.
  • arg (c) represents the argument of the complex number c.
  • a frequency image can be calculated by the following equation (10). Therefore, a process of acquiring a frequency image may be performed separately from the water / fat separation process by the fitting described above.
  • Magnetic field calculation processing (S1202) Next, a magnetic field calculation process described below is performed on the calculated frequency image ⁇ ′ to calculate a magnetic field image ⁇ (step S1202).
  • the magnetic field image is an image representing a relative magnetic field change normalized by the static magnetic field intensity.
  • the magnetic field calculation unit 520 calculates the magnetic field image ⁇ (r) by converting the frequency image ⁇ ′ (r) by the following equation (11).
  • is the proton gyromagnetic ratio
  • B 0 is the static magnetic field strength
  • the mask calculation unit 530 of the present embodiment includes a low signal region mask representing a low signal region such as air from the complex image calculated by the image reconstruction unit 400 and the water image / fat image calculated by the water / fat separation processing unit 510, A water mask representing a region where water is dominant in the voxel and a fat mask representing a region where fat is dominant in the voxel are calculated. Therefore, the mask calculation unit 530 includes a low signal region mask calculation unit 531, a fat mask calculation unit 532, and a water mask calculation unit 533, as shown in FIG.
  • FIG. 9 shows a process flowchart of the mask calculation unit 530 of the present embodiment.
  • a mask absolute value image is calculated from a complex image, and a low signal region mask and a high signal region mask are calculated from the mask absolute value image (step S1401).
  • a fat mask is calculated from the high signal region mask and the fat content rate image calculated in step S1304 (step S1402).
  • a water mask is calculated from the high signal area mask and the fat mask (step S1403). Details of each process will be described below.
  • the mask absolute value image is an image in which each pixel is represented by a real number of 0 or more.
  • the absolute value image for masking in the present embodiment is a single absolute value image obtained by calculating the square sum square in the time direction for each pixel from the absolute value image of all echoes. Note that the mask absolute value image is calculated by an arbitrary method. For example, an absolute value image of the first echo may be used.
  • a low signal area mask is calculated from the absolute value image for masking.
  • the low signal area mask Ma is calculated from the mask absolute value image by threshold processing.
  • the low signal area mask uses a predetermined threshold value, and an area where the pixel value of the mask absolute value image is smaller than the threshold value is 1 and a large area is 0.
  • the threshold value may be obtained from a pixel value distribution of all pixels of the absolute value image using a method such as a discriminant analysis method. Further, after threshold processing, image processing such as isolated point removal processing may be performed.
  • image processing such as isolated point removal processing may be performed.
  • There are various methods for calculating the low signal region mask For example, a boundary part inside and outside the body may be detected, and the internal region may be set to 0 and the external region may be set to 1 based on the detection result. Such other methods may be used in the noise mask process.
  • a value obtained by subtracting the pixel value of the low signal area mask from 1 for each pixel is calculated as the high signal area mask Mh.
  • the fat mask calculation method in step S1402 in order to divide the high signal area mask Mh into the fat mask Mf and the water mask Mw, an image (discrimination image) for discriminating both is prepared, and a discrimination mask is created.
  • the discrimination mask Md is calculated by using the fat content rate image calculated in the water / fat separation process S1201 (S1304) as the discrimination image and performing threshold processing using an arbitrary threshold. For example, in the fat content rate image, a region where the fat content rate is larger than the threshold value is 1, and a region where the fat content rate is smaller than the threshold value is 0. In this embodiment, the threshold is set to 0.8.
  • the fat mask Mf is calculated by multiplying the discrimination mask Md thus calculated and the high signal region mask Mh.
  • the method for calculating the fat mask is arbitrary, and for example, the fat mask may be calculated by a method such as threshold processing of a fat image.
  • the fat content image calculation process (FIG. 8, S1304) in the water / fat separation process (FIG. 6, S1201) can be omitted.
  • processing for removing isolated points or processing for setting fat regions below an arbitrary threshold to 0 may be added to form a fat mask. By adding this processing, it is possible to remove the region determined to be a fat region by noise.
  • the water mask Mw is calculated (step S1403).
  • the water mask Mw is calculated by subtracting the fat mask Mf from the high signal region mask Mh.
  • the calculation method of the water mask is arbitrary, and for example, the water mask may be calculated by a method such as thresholding a water image.
  • the pixel values of the low signal area mask Ma, the fat mask Mf, and the water mask Mw are set to 0 or 1, but can be set to arbitrary values.
  • the fat mask may be calculated by multiplying the fat content image by the high signal region mask, and the water mask by multiplying the image obtained by subtracting the fat content rate from 1 for each pixel by the high signal region mask.
  • the magnetic susceptibility calculation unit 540 of the present embodiment uses the low signal area mask, water mask, and fat mask calculated by the mask calculation unit 530 and the magnetic field image calculated by the magnetic field calculation unit 520 to magnetize. A rate image is calculated. Therefore, the magnetic susceptibility calculation unit 540 includes a background magnetic field removal unit 541, a fat magnetic susceptibility calculation unit 542, a background magnetic field removal unit 543, a water magnetic susceptibility calculation unit 544, and an integrated magnetic susceptibility calculation unit 545, as shown in FIG. .
  • Step S1501 to S1505 in FIG. 10 correspond to the processes of the units 541 to 545 of the magnetic susceptibility calculation unit 540, respectively.
  • background removal processing is performed on a magnetic field image, assuming that an area defined by a low signal area mask is a background area (step S1501).
  • a magnetic susceptibility calculation process is performed from the image whose background has been removed in S1501, and a fat magnetic susceptibility image is calculated under the constraint that the water region is 0 (step S1502).
  • a background removal process is performed on the magnetic field image in which the area defined by the low signal area mask and the fat mask is regarded as the background area (step S1503).
  • a magnetic susceptibility calculation process is performed from the image with the background removed in S1503, and a water susceptibility image is calculated (step S1504).
  • an integrated magnetic susceptibility image is calculated by integrating the fat magnetic susceptibility image and the water magnetic susceptibility image (step S1505). Details of each process will be described below.
  • step S1501 background magnetic field removal is performed on the magnetic field image ⁇ calculated in the magnetic field calculation process S1202 with the region defined by the low signal region mask regarded as the background region.
  • the background magnetic field can be removed using any method, but the SHARP method is used in this embodiment.
  • the local magnetic field is calculated by solving the above equation by the singular value decomposition method or the least square method.
  • M is an orthogonal matrix having a binary mask (high signal mask) with a region of interest of 1 as a diagonal component
  • H is a matrix representing a convolution operation
  • ⁇ local is a local magnetic field image (magnetic field ⁇ after background removal) ')
  • ⁇ total is the column vector of the magnetic field image before the background removal (the magnetic field image ⁇ calculated in the magnetic field calculation processing S1202).
  • the region of interest is a region defined by the diagonal component of M
  • the background region is a region defined by the diagonal component of IM where I is a unit matrix.
  • step S1502 from the magnetic field image from which the background magnetic field has been removed in step S1501, a smoothing constraint that places a constraint on smoothing the magnetic susceptibility and a specific value constraint that places a restriction on the magnetic susceptibility of a specific region are set.
  • a fat susceptibility image ⁇ f is calculated.
  • the specific area is an area defined by a water mask, and the specific value is 0.
  • C represents an operator representing a convolution operation
  • G represents a differential operator
  • M is a high signal area mask
  • ⁇ ′ is a magnetic field whose background is removed in S1501
  • Mw is a water mask
  • Mf is a fat mask.
  • the first term in equation (13) is a term that minimizes the error in the relational expression between the magnetic susceptibility and frequency with the mask M as the region of interest
  • the second term is a term that places a constraint that the water region is 0,
  • the third term is It is a term that smoothes the fat.
  • ⁇ w and ⁇ f define the size of each constraint and are arbitrary constants.
  • ⁇ w is called a specific value parameter
  • ⁇ f is called a smoothing parameter.
  • the specific value parameter ⁇ w that defines the size of the specific value constraint and the smoothing parameter ⁇ f that defines the size of the smoothing constraint are:
  • the fat susceptibility image is set to a value that maximizes the ratio between the average value and the standard deviation in the fat mask.
  • Non-Patent Document 3 discloses a method of calculating the magnetic susceptibility in the brain under the constraint that the magnetic susceptibility outside the brain, which is the background region, is zero after removing the background magnetic field that regards the region outside the brain as the background region. Has been.
  • regularization parameter of regularization term to be used for the constraint is as large as 10 20.
  • the water area to which the restriction is applied is not the background area.
  • ⁇ w the value of ⁇ w to be equal to or less than a predetermined value ( ⁇ wth )
  • ⁇ wth a predetermined value
  • Expression (13) can be minimized by iterative calculation using a conjugate gradient method or the like.
  • the initial vector is arbitrary, but for example, a zero vector is used.
  • the end condition of the repetitive calculation is arbitrary. For example, there is a method of ending the repetitive calculation when the number of repetitions exceeds a predetermined arbitrary threshold.
  • the calculation accuracy of fat magnetic susceptibility can be improved by calculating the magnetic susceptibility using the second term under the constraint that unevenness does not occur in the water region.
  • the weight W may be used in place of the mask M in the first item in the equation (13).
  • the weight W is normally defined within the region of interest excluding the background region, and is calculated based on, for example, a phase image (the method described in Japanese Patent Laid-Open No. 2015-62637 by the present inventors). Specifically, the phase variation of each pixel is calculated from the phase image P, and the weight W is calculated based on the calculated phase variation. For example, the weight W is calculated to be smaller as the phase variation is larger.
  • a statistical index such as a standard deviation calculated using pixel values (phase values) of a plurality of pixels in a peripheral region of the target pixel for which the phase variation is calculated is used.
  • the following formula (14) may be used in addition to the formula (13).
  • Expression (14) by using Mf for the first term, it is possible to place a specific value constraint that makes the magnetic susceptibility of the water region zero.
  • the second term is the same smoothing constraint as the third term of Equation (13).
  • Mf ′ a weight function Mf ′ in which the water region has an arbitrary constant c w (for example, 0.2) and the fat region has an arbitrary constant c f (for example, 0.8) may be used. Good (but c w ⁇ c f ).
  • a method of alternately updating a solution based on an error function defined by a predetermined initial susceptibility distribution ⁇ f0 and a magnetic field distribution ⁇ and a smoothing process (according to the present inventors)
  • the magnetic susceptibility image may be calculated based on the method described in Japanese Patent Application No. 2014-228843: International Publication No. 2014076076 (hereinafter referred to as “new magnetic susceptibility calculation method”).
  • an updated magnetic susceptibility distribution ⁇ UPD is calculated from the magnetic field distribution ⁇ and the initial magnetic susceptibility distribution ⁇ f0 , the updated magnetic susceptibility distribution ⁇ UPD is smoothed, and the filtered magnetic susceptibility distribution (post-processed magnetic susceptibility) is calculated.
  • Distribution) ⁇ SMT is calculated, the updated magnetic susceptibility distribution ⁇ UPD and the processed magnetic susceptibility distribution ⁇ SMT are combined to calculate a composite magnetic susceptibility distribution ⁇ ADD, and it is determined whether or not the repetitive operation is completed, and the initial magnetic susceptibility distribution The process of updating ⁇ f0 to the composite magnetic susceptibility distribution ⁇ ADD is repeatedly performed.
  • the updated susceptibility distribution ⁇ UPD is calculated from the initial susceptibility distribution ⁇ f0 using an error function e ( ⁇ f0 ) shown in the following equation.
  • the gradient ⁇ e ( ⁇ f0 ) of the error function e ( ⁇ f0 ) is calculated from the equation (15), and the updated magnetic susceptibility is calculated by the following equation (16) using the calculated gradient ⁇ e ( ⁇ f0 ).
  • the initial susceptibility distribution ⁇ f0 to be set first can be set arbitrarily.
  • an initial magnetic susceptibility distribution ⁇ f0 may be obtained by setting all values to 0.
  • Mf in equation (15) it is possible to place a constraint that the susceptibility of the water region is zero.
  • a weighting function Mf ′ in which the water region has an arbitrary constant c w (for example, 0.2) and the fat region has an arbitrary constant c f (for example, 0.8) is used. (But c w ⁇ c f ).
  • the post-processing magnetic susceptibility distribution ⁇ SMT is calculated by applying an edge-preserving smoothing filter (filter described in International Publication No. 20116077606) to the updated magnetic susceptibility distribution ⁇ UPD and performing smoothing (processing corresponding to smoothness constraint) ).
  • This edge-preserving smoothing filter is a smoothing filter that has the property that smoothing works strongly in a region where local dispersion is small and smoothing works weakly in a region where local dispersion is large, and can maintain a quantitative value of magnetic susceptibility.
  • the edge portion including the fine section image is preserved. Note that the smoothing process may be performed only on the fat region.
  • the combined magnetic susceptibility distribution ⁇ ADD is calculated by combining the updated magnetic susceptibility distribution ⁇ UPD and the processed magnetic susceptibility distribution ⁇ SMT .
  • each of the updated susceptibility distribution X UPD and the processed susceptibility distribution X SMT is first converted into k-space data by Fourier transform. Then, each k-space data (Fourier component) after conversion is weighted, and each weighted k-space data is added. Then, after addition, a composite magnetic susceptibility distribution ⁇ ADD is obtained by inverse Fourier transform.
  • Weighting in the k space is performed by creating a weighted image G (k) (k indicates a data position in the k space; hereinafter the same) and multiplying the weighted image G by the k space data.
  • the weighted image G UPD for multiplying the k-space data obtained from the updated magnetic susceptibility distribution X UPD for example, the weight given to the data in the vicinity of the magic angle area is set smaller than the weight given to the data in other areas Is done.
  • a region where the Fourier component value is smaller than the predetermined value (b) is defined as a magic angle region.
  • the minimization process for minimizing the error function e ( ⁇ f0 ) by iterative computation it is determined whether or not the iterative computation has ended. Continuing the process, i.e., if it is determined that repeated, updating the latest synthetic magnetic susceptibility distribution chi ADD obtained at that time, by setting the initial magnetic susceptibility distribution chi f0, the magnetic susceptibility distribution chi f0 To do.
  • the latest combined magnetic susceptibility distribution ⁇ ADD at that time is output as the magnetic susceptibility distribution (here, fat susceptibility distribution).
  • the evaluation function f ( ⁇ UPD ) shown in the following equation (18) is used, and when the value is smaller than a predetermined threshold ⁇ , it is determined that the process is finished.
  • step S1502 The specific example of the intra-fat susceptibility calculation processing in step S1502 has been described above. However, in step S1502, a smoothing constraint for smoothing the susceptibility and a constraint for setting the susceptibility of a specific region as a specific value are specified. Any method that calculates the fat susceptibility image ⁇ f under the value constraint can be adopted without being limited to these specific examples.
  • the resolution of the input image (magnetic field image after background removal, high signal mask image, fat mask image, and water mask image) such as ⁇ ′, M, Mf, and Mw is calculated by a method such as interpolation.
  • the magnetic susceptibility image ⁇ f may be calculated from the end, and the resolution of the magnetic susceptibility image ⁇ f may be returned to the original size by a method such as interpolation after the calculation.
  • ⁇ f may be calculated after halving the image size in each direction of each input image, and the image size in each direction may be returned to the original size after the calculation.
  • step S1503 background magnetic field removal is performed on the magnetic field image ⁇ calculated in the magnetic field calculation process S1202 with the region defined by the sum of the low signal region mask and the fat mask regarded as the background region (step S1503).
  • the background magnetic field removal can be performed by an arbitrary method such as the SHARP method or the PDF method.
  • the background magnetic field removal in step S1503 is performed by using the background magnetic field removal described in step S1501 except that a region defined by the sum of the low signal region mask and the fat mask is used as the background region. It is the same.
  • step S1504 a magnetic susceptibility calculation process is performed using the region defined by the water mask as a region of interest from the magnetic field image from which the background magnetic field has been removed in S1503, thereby calculating a water susceptibility image ⁇ w.
  • the magnetic susceptibility can be calculated by an arbitrary method in the same manner as in step S1502 of the fat magnetic susceptibility calculation process.
  • the magnetic susceptibility is calculated by a new magnetic susceptibility calculation method (method described in International Publication No. 20106077606). Is calculated.
  • the constraint condition based on the relational expression between the magnetic field and the magnetic susceptibility is defined in a region where the water mask is 1. That is, the following equation (19) is used instead of the above-described equation (15) used for calculating the fat susceptibility.
  • the other calculation method is the same as the method described in the fat susceptibility calculation method, and repeats the process of alternately updating the predetermined initial susceptibility distribution ⁇ w0 and the smoothing process for a predetermined number of times.
  • the combined magnetic susceptibility distribution ⁇ ADD obtained after repeated calculations is taken as the water magnetic susceptibility.
  • step S1505 the fat magnetic susceptibility (intrafat magnetic susceptibility image) obtained in S1502 and the water magnetic susceptibility (external fat susceptibility image) obtained in step S1504 are integrated by equation (20).
  • a magnetic susceptibility image ⁇ is calculated (step S1505).
  • the magnetic susceptibility calculation processing (FIG. 6: S1204) is completed, and the image conversion processing step S1103 in FIG. 4 is completed.
  • Predetermined image processing is performed on the integrated magnetic susceptibility image as necessary to create a magnetic susceptibility image (absolute value image) or a magnetic susceptibility-enhanced image and display it on the display unit 600 (step S1104). It is as follows.
  • the MRI apparatus of the present embodiment includes a static magnetic field magnet, a magnetic field generation unit that applies a gradient magnetic field and a high-frequency magnetic field to a space formed by the static magnetic field magnet, and a subject placed in the space.
  • a reception unit that receives a nuclear magnetic resonance signal generated from the computer, and a computer that performs an operation on the nuclear magnetic resonance signal.
  • the computer controls operations of the magnetic field generation unit and the reception unit, and includes a A measurement unit that measures a nuclear magnetic resonance signal including a uniform influence, an image reconstruction unit that performs image reconstruction using the nuclear magnetic resonance signal and creates a complex image, and the subject using the complex image
  • An image conversion unit uses a magnetic field calculation unit that calculates a magnetic field reflection image from the complex image and an image created by the image reconstruction unit, and a low signal region mask corresponding to a low signal region and a plurality of signals in the high signal region Using a mask calculation unit that calculates a plurality of region masks corresponding to the region, a magnetic field reflection image calculated by the magnetic field calculation unit, the low signal region mask and a plurality of region masks calculated by the mask calculation unit, A susceptibility calculating unit that calculates a susceptibility for each of the plurality of regions.
  • the magnetic susceptibility calculation unit is configured to determine a magnetic susceptibility of a region different from the specific region under a restriction with the low signal region as a background and a restriction with a specific value of the magnetic susceptibility of the specific region among the plurality of regions. Is calculated.
  • the magnetic susceptibility calculation unit calculates a magnetic susceptibility of the specific area under a restriction with a background different from the low signal area and the specific area, and a magnetic susceptibility of the area different from the specific area. Integrate with.
  • the image conversion unit of the present embodiment uses the image created by the image reconstruction unit to separate a water image in which a signal from water protons is dominant and a fat image in which a signal from fat protons is dominant.
  • a water fat separation unit is further provided, and the mask calculation unit uses a predetermined value of the fat content for each pixel calculated using the water image and the fat image as a threshold, and water is controlled as the plurality of region masks.
  • a water mask corresponding to a general region and a fat mask corresponding to a fat-dominated region are created.
  • the MRI apparatus when calculating the magnetic susceptibility distribution for each of the fat image and the water image, has the restriction that the susceptibility unevenness does not occur in the water region in the fat region, that is, the magnetic susceptibility.
  • the magnetic susceptibility By calculating the magnetic susceptibility under the constraint that becomes 0, the accuracy loss due to signal loss is reduced.
  • the background magnetic field is removed with the air region and the fat region as the background to calculate the magnetic susceptibility, so that the magnetic field change due to the water fat magnetic susceptibility difference is removed, and unevenness does not occur. That is, the magnetic susceptibility image can be calculated with high accuracy in the fat region and the water region.
  • any imaging section such as a transverse section, a coronal section, a sagittal section, and an oblique section, and the same effect can be obtained.
  • the present invention is not limited to this. At least one of these units may be constructed on an information processing apparatus independent of the MRI apparatus capable of transmitting and receiving data to and from the computer 209 of the MRI apparatus, for example.
  • the magnetic susceptibility calculation unit 540 calculates the magnetic susceptibility for each of the fat region and the water region.
  • the susceptibility of the fat region is determined in advance. It can also be calculated as a specific value.
  • the susceptibility calculation unit 540 of FIG. 3 the fat susceptibility calculation unit 542 can be omitted, and the background magnetic field removal processing S1501 and the fat susceptibility calculation processing S1502 are not performed in the flow shown in FIG. It is changed as follows.
  • the magnetic susceptibility calculation unit 540 first performs background removal in which the region defined by the low signal region mask is regarded as the background region (S1503), and then calculates the magnetic susceptibility image ⁇ using the following equation (S1504). .
  • M is a high signal area mask
  • Mf is a fat mask
  • ⁇ f is an arbitrary constant
  • e is a unit vector.
  • ⁇ f a predetermined difference in magnetic susceptibility between water and fat is used. For example, it is 0.61 ppm.
  • the second term represents a constraint for smoothing the magnetic susceptibility.
  • the third term represents the constraint that the susceptibility in the fat region is a specific value ⁇ f.
  • ⁇ and ⁇ fat are constants representing the strengths of the constraints of the second term and the third term, respectively.
  • Equation (21) can be minimized by iterative calculation using a conjugate gradient method or the like.
  • the initial vector is arbitrary, but for example, a zero vector is used.
  • the end condition of the repetitive calculation is arbitrary. For example, there is a method of ending the repetitive calculation when the number of repetitions exceeds a predetermined arbitrary threshold.
  • the formula for calculating the magnetic susceptibility of fat as a predetermined specific value is not limited to the formula (21).
  • a constraint based on the L1 norm may be used for the second term.
  • the magnetic susceptibility of the water region can be calculated with high accuracy by adding the restriction of the third term.
  • the calculation of the integrated magnetic susceptibility image ⁇ using the equation (20) for the calculated magnetic susceptibility of the water region and the magnetic susceptibility (specific value) of fat is the same as in the first embodiment.
  • Second Embodiment In the first embodiment, two high signal area masks, a fat mask and a water mask, are calculated by thresholding the fat content image, and the magnetic susceptibility is calculated by a different method in each area, and finally integrated.
  • the magnetic susceptibility image was calculated by
  • two high-signal area masks are calculated by thresholding not a fat content image but a previously calculated magnetic susceptibility image. Then, the magnetic susceptibility image is calculated by performing the same processing as in the first embodiment.
  • the present embodiment will be described focusing on differences from the first embodiment.
  • the MRI apparatus of this embodiment basically has the same configuration (FIG. 2) as that of the first embodiment.
  • the computer 209 of this embodiment instructs the sequencer 204 to perform echo measurement according to the measurement parameter and pulse sequence set (or set in advance) via the input device 212.
  • the echoes are measured, the obtained echoes are arranged in the k space, a calculation is performed on the echoes arranged in the k space, a magnetic susceptibility image is calculated, and the obtained image is displayed on the display device 210.
  • water fat separation is assumed, and two or more echoes having different echo times are measured.
  • the measurement may be performed with only one echo time.
  • the computer 209 of the present embodiment also includes the measurement unit 300, the image reconstruction unit 400, the image conversion unit 500, and the display processing unit 600, as in the first embodiment, but the configuration of the image conversion unit 500 is the same. Different.
  • FIG. 11 is a functional block diagram of the total image conversion unit 500 in the present embodiment. In FIG. 11, elements having the same functions as those in the first embodiment are denoted by the same reference numerals as those in FIG.
  • the image conversion unit 500 of this embodiment includes a magnetic field calculation unit 520, a mask calculation susceptibility calculation unit 550, a mask calculation unit 530, and a susceptibility calculation unit 540.
  • the image conversion unit 500 of the present embodiment may include the water / fat separation processing unit 510 as in the first embodiment.
  • the mask calculation magnetic susceptibility calculation unit 550 calculates a magnetic susceptibility image necessary for the processing in the mask calculation unit 530. Unlike the final magnetic susceptibility calculated by the magnetic susceptibility calculator 540, this mask-calculating magnetic susceptibility image is a provisional magnetic susceptibility image for making a mask.
  • the mask calculation unit 530 includes a low signal region mask calculation unit 810, a high magnetic susceptibility mask calculation unit 820, and a low magnetic susceptibility mask calculation unit 830.
  • the magnetic susceptibility calculation unit 540 includes background magnetic field removal units 841 and 843, a high magnetic susceptibility (first magnetic susceptibility image) calculation unit 842, a low magnetic susceptibility (second magnetic susceptibility image) calculation unit 844, and an integrated magnetic susceptibility calculation. Part 845.
  • the image conversion unit 500 of the present embodiment performs a magnetic field calculation process S1201, a mask calculation susceptibility image calculation process S1210, a mask calculation process S1203, and a susceptibility calculation process S1204 as shown in FIG.
  • Magnetic field calculation processing (S1202)
  • a frequency image is calculated by fitting, or a frequency image is calculated using the above formula (9) or (10).
  • a magnetic field image is calculated using a relational expression with the magnetic field (formula (11)).
  • the number of echoes having different echo times may be two or more, and water fat separation may be performed together with the calculation of the frequency image.
  • a provisional magnetic susceptibility image is calculated using the complex image calculated by the image reconstruction unit 400 (processing of the mask calculation magnetic susceptibility calculation unit 550).
  • the magnetic susceptibility image for mask calculation can calculate the magnetic susceptibility by, for example, a “new magnetic susceptibility calculation method” (method described in International Publication No. 2014076076).
  • the constraint condition based on the relational expression between the magnetic field and the magnetic susceptibility is defined in a region where the high signal region mask is 1.
  • step S1203 using the complex image calculated by the image reconstruction unit 400, the magnetic field image calculated by the magnetic field calculation unit 520, and the mask calculation susceptibility image calculated by the mask calculation susceptibility calculation unit 550, three masks, That is, a low signal region mask representing a low signal region such as air, a high magnetic susceptibility mask representing a relatively high magnetic susceptibility region, and a low magnetic susceptibility mask representing a relatively low magnetic susceptibility region are calculated.
  • FIG. 13 shows an example of a flowchart of the mask calculation process.
  • a low signal region mask and a high signal region mask are calculated by the same method as in the first embodiment (step S1601).
  • a high magnetic susceptibility mask is calculated from the high signal region mask and the mask calculating magnetic susceptibility image (step S1602).
  • a low magnetic susceptibility mask is calculated from the high signal region mask and the high magnetic susceptibility mask (step S1603). Details of each process will be described below.
  • the low signal region mask calculation method in step S1601 is basically the same as in the first embodiment. That is, the mask absolute value image is calculated from the complex image, and the low signal region mask and the high signal region mask are calculated from the mask absolute value image.
  • step S1602 the mask magnetic susceptibility image calculated in step S1210 is used as a discrimination image for dividing the high signal region, and a mask is calculated by performing threshold processing using an arbitrary threshold. That is, an area where the calculated magnetic susceptibility is higher than the threshold in the mask magnetic susceptibility image is 1, and an area where the calculated magnetic susceptibility is lower than the threshold is 0.
  • the threshold value is arbitrary, in this embodiment, it is set to 0.5 ppm, for example.
  • a high susceptibility mask is calculated by multiplying the calculated mask and the high signal region mask.
  • a process for removing isolated points or a process for setting a high magnetic susceptibility region whose area is equal to or smaller than an arbitrary threshold to 0 may be added to form a high magnetic susceptibility mask.
  • a low magnetic susceptibility mask is calculated (step S1603).
  • the low magnetic susceptibility mask is calculated by subtracting the high magnetic susceptibility mask from the high signal region mask.
  • the pixel values of the low signal region mask, the high magnetic susceptibility mask, and the low magnetic susceptibility mask are set to 0 or 1, but may be arbitrary values.
  • filter processing using a moving average filter or the like may be added to each mask.
  • the magnetic susceptibility calculation process (step S1204) is performed.
  • the first signal is calculated from the low signal region mask, the high magnetic susceptibility mask, and the low magnetic susceptibility mask calculated by the mask calculation unit 530 and the magnetic field image calculated by the magnetic field calculation unit 520.
  • the magnetic susceptibility image is calculated by the same method as in the embodiment. However, at this time, unlike the first embodiment, a high magnetic susceptibility mask is used instead of the fat mask, and a low magnetic susceptibility mask is used instead of the water mask.
  • the background removal processing unit 841 performs background removal processing on the magnetic field image, regarding the region defined by the low signal region mask as the background region.
  • a magnetic susceptibility calculation process is performed from the background-removed image, and the magnetic susceptibility image of the region defined by the high magnetic susceptibility mask (first magnetization) under the constraint that the region defined by the low magnetic susceptibility mask becomes zero Ratio image) (processing of high magnetic susceptibility calculation unit 842).
  • the background removal processing unit 843 performs background removal processing on the magnetic field image, assuming that the region defined by the low signal region mask and the high magnetic susceptibility mask is the background region.
  • a magnetic susceptibility calculation process is performed from the image from which the background has been removed by the background removal processing unit 843 to calculate a magnetic susceptibility image (second magnetic susceptibility image) in a region defined by the low magnetic susceptibility mask (low magnetic susceptibility calculation Section 844).
  • an integrated susceptibility image is calculated by integrating the first susceptibility image and the second susceptibility image (processing of the integrated susceptibility calculation unit 845).
  • the mask calculation unit uses the discrimination magnetic susceptibility image calculated using the magnetic field image, the high magnetic susceptibility region mask corresponding to the high magnetic susceptibility region, and A region mask corresponding to a region other than the high magnetic susceptibility region is created.
  • the mask calculation unit uses the discrimination magnetic susceptibility image calculated using the magnetic field image, the high magnetic susceptibility region mask corresponding to the high magnetic susceptibility region, and A region mask corresponding to a region other than the high magnetic susceptibility region is created.
  • the water image and the fat image are not used in the mask calculation unit 530, the water / fat separation process is not necessarily performed.
  • Third embodiment In the first and second embodiments, two high-signal area masks are calculated, magnetic susceptibility images are calculated by different methods in the respective areas, and finally integrated susceptibility images are calculated. In the third embodiment, three high signal area masks are calculated, magnetic susceptibility images are calculated by different methods in the respective areas, and finally integrated susceptibility images are calculated by integration.
  • a description will be given focusing on the configuration different from the first embodiment and the second embodiment.
  • the areas defined by the three high signal area masks are not limited.
  • an area having a high magnetic susceptibility such as bleeding, a fat area, and other areas can be assumed. In the following embodiment, these will be described as examples.
  • the MRI apparatus of this embodiment basically has the same configuration as that of the first embodiment.
  • the computer 209 of this embodiment instructs the sequencer 204 to perform echo measurement according to the measurement parameter and pulse sequence set (or set in advance) via the input device 212.
  • the echoes are measured, the obtained echoes are arranged in the k space, a calculation is performed on the echoes arranged in the k space, a magnetic susceptibility image is calculated, and the obtained image is displayed on the display device 210.
  • the computer 209 of this embodiment includes a measurement unit 300, an image reconstruction unit 400, an image conversion unit 500, and a display processing unit 600, as in the first embodiment.
  • FIG. 14 is a functional block diagram of the image conversion unit 500 in the present embodiment.
  • the image conversion unit 500 of the present embodiment includes a water / fat separation processing unit 910, a magnetic field calculation unit 920, a mask calculation unit 930, and a magnetic susceptibility calculation unit 940.
  • the mask calculation unit 930 is a low signal region mask calculation unit 931, a first high signal region mask calculation unit 932, a second high signal region mask calculation unit 933, and a third high signal region.
  • a mask calculation unit 934 is provided.
  • the magnetic susceptibility calculation unit 940 includes a background magnetic field removal unit 941, a first high signal magnetic susceptibility calculation unit 942, a background magnetic field removal unit 943, a second high signal magnetic susceptibility calculation unit 944, and a background.
  • a magnetic field removal unit 945, a third high signal susceptibility calculation unit 946, and an integrated susceptibility calculation unit 947 are provided.
  • the image conversion unit 500 of the present embodiment includes a water / fat separation process S1201 for calculating a water image / fat image / frequency image and a fat content image, and a magnetic field calculation process for calculating a magnetic field image from the frequency image.
  • S1202 a process for calculating a temporary magnetic susceptibility image for use in mask calculation from a magnetic field image (masking magnetic susceptibility image calculation process) S1210, a process for calculating a plurality of area masks (mask calculation process) S1213, a magnetic field image
  • the mask calculation unit 930 uses the complex image calculated by the image reconstruction unit 400, the fat content ratio image calculated by the water fat separation processing unit 910, the magnetic field image calculated by the magnetic field calculation unit 920, air, and the like.
  • a low signal region mask that represents a low signal region a first high signal region mask that represents a relatively high magnetic susceptibility region, a second high signal region mask that represents a fat region, and a second high signal region mask that represents the other high signal region.
  • Three high signal region masks are calculated.
  • FIG. 16 is a flowchart of the mask calculation process of this embodiment
  • FIG. 17 is a diagram showing the relationship between the image used for mask calculation and the calculated mask.
  • the mask calculating process first calculates the low-signal region mask M L and a high signal area mask M H (step S1701). Then calculated first high signal area mask M 1 from the high signal area mask M H and mask susceptibility image (step S1702). Next, the second high signal region mask M 2 is calculated from the high signal region mask M H , the first high signal region mask M 1 , and the fat content image (step S1703). Finally, the high signal area mask M H, a first high signal area mask M 1, the second high signal area mask M 2, to calculate a third high signal area mask M 3 (step S1704). Details of each process will be described below.
  • the low signal area mask calculation method in step S1701 is basically the same as that in the first embodiment.
  • the absolute value image for mask is calculated from the complex image, the absolute value image for mask is thresholded, and the low signal area mask is calculated. calculating the M L and a high signal area mask M H (processing of the low-signal area mask calculating section 931).
  • step S1702 the mask susceptibility image calculated in step S1210 is used as a discrimination image in the same manner as the high susceptibility mask calculation method in the second embodiment, and threshold processing is performed using an arbitrary threshold value.
  • a mask (first high signal area mask) M 1 is calculated (processing of the first high signal mask calculation unit 932).
  • step S1703 calculates a second high signal area mask M 2. Therefore, first, the mask Md is calculated by performing threshold processing on the fat content image calculated in the water / fat separation processing S1201 using an arbitrary threshold. For example, a region where the fat content is larger than the threshold is 1, and a region where the fat content is smaller than the threshold is 0. In this embodiment, the threshold is set to 0.8. Then, a mask Md calculated by threshold processing, by multiplying the high signal area mask M H masked by subtracting the first high signal area mask M 1 from M (H-1), a second high signal area mask to calculate the M 2.
  • the calculated four masks take only the values 0 and 1. However, similar to the first embodiment and the second embodiment, the calculated four masks may have arbitrary values. Further, as in the first embodiment and the second embodiment, additional processing such as isolated point removal processing may be added to each calculated mask.
  • the magnetic susceptibility calculation unit 940 includes a low signal region mask ML, a first high signal region mask M1, a second high signal region mask M2, a third high signal region mask M3, and a magnetic field calculation unit calculated by the mask calculation unit 930.
  • a magnetic susceptibility image is calculated from the magnetic field image calculated in 920.
  • FIG. 18 is a flowchart of the magnetic susceptibility calculation process of the present embodiment.
  • a background removal process is first performed on the magnetic field image, assuming that the region defined by the low signal region mask is the background region (step S1801).
  • a magnetic susceptibility calculation process is performed from the magnetic field image with the background removed in step S1801, and the magnetic susceptibility of the region defined by the sum of the second high signal region mask M2 and the third high signal region mask M3 is set to zero.
  • a first high signal susceptibility image is calculated (step S1802).
  • a background removal process is performed on the magnetic field image in which the region defined by the low signal region mask ML and the first high signal region mask M1 is regarded as the background (step S1803).
  • a second high signal susceptibility image is calculated from the magnetic field image with the background removed in step S1803 under the constraint that the third high signal region mask M3 is 0 (step S1804).
  • a background removal process is performed in which the area defined by the low signal area mask ML, the first high signal area mask M1, and the second high signal area mask M2 is regarded as the background (step S1805).
  • a magnetic susceptibility image is calculated from the image whose background has been removed in step S1805 (step S1806).
  • the background magnetic field removal method in steps S1801, S1803, and S1805 is basically the same as the method described in the first embodiment except for the difference in the background region.
  • the method for calculating the magnetic susceptibility image under the restriction that the magnetic susceptibility of some regions in step S1802 and step S1804 is basically the first except for the difference between the region of interest and the region to be restricted. This is the same as the method described in the embodiment.
  • the calculation method of the magnetic susceptibility image in step S1806 is basically the same as the method in step S1504 in the first embodiment except for the difference in the region of interest.
  • the first high signal susceptibility image, the second high signal susceptibility image, and the third high signal susceptibility image are integrated by Expression (22) to calculate an integrated susceptibility image ⁇ (step S1807). ).
  • M 1 , M 2 , and M 3 are the first high signal region mask, the second high signal region mask, and the third high signal region mask, and ⁇ 1 , ⁇ 2 , and ⁇ 3 are the first A high signal susceptibility image, a second high signal susceptibility image, and a third high signal susceptibility image are shown.
  • the image conversion unit uses the image created by the image reconstruction unit, and the water image, fat, or the like, in which signals from water protons are dominant among the components included in the subject.
  • the apparatus further includes a water / fat separation unit that separates a fat image in which signals from protons are dominant, and the mask calculation unit calculates the fat content and the magnetic field image calculated for each pixel using the water image and the fat image.
  • the discrimination magnetic susceptibility image calculated by using the high magnetic susceptibility region mask corresponding to the high magnetic susceptibility region, the fat mask corresponding to the fat-dominated region, and the high magnetic susceptibility region mask and the fat
  • An area mask corresponding to an area other than the area defined by the mask is created.
  • a high accuracy magnetic susceptibility image is calculated by reducing the accuracy degradation due to signal loss. can do.
  • three high signal area masks may be calculated by thresholding the mask susceptibility image using two threshold values.
  • Example of the first embodiment As shown in Table 1 below, a model is assumed in which the susceptibility in fat is 0.61 and the susceptibility outside fat is 0.
  • the susceptibility calculation method (invention method) of the first embodiment and the conventional method are used.
  • the magnetic susceptibility image was obtained by simulation.
  • the conventional method employs a method of calculating the magnetic susceptibility under the constraint of smoothing the magnetic susceptibility of the high signal region after removing the background of the low signal region.

Landscapes

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

Abstract

信号欠損による精度低下を低減して磁化率を高精度に算出する。このため、MRIを用いて空間的磁場不均一を反映したエコーを少なくとも一つ取得し、取得したエコーから複素画像を算出し、算出した複素画像から低信号領域をあらわす低信号領域マスク、高信号かつ脂肪含有率の大きい領域をあらわす第一の高信号領域マスク、高信号かつ脂肪含有率の小さい領域をあらわす第二の高信号領域マスクの三つのマスクを算出する。複素画像から作成した周波数画像或いは磁場画像から磁化率画像を作成する際、低信号領域マスクで定義される領域を背景、第二の高信号領域マスクで定義される領域の磁化率を予め定めた特定値とする制約のもと、磁化率画像を算出する。

Description

磁気共鳴イメージング装置、磁気共鳴イメージング方法及び磁化率算出プログラム
 本発明は、磁気共鳴イメージング(MRI)技術に関する。特に、再構成後の画像の画像処理技術に関する。
 磁気共鳴イメージング装置(以下、MRI装置ともいう) は、静磁場内に置かれた被検体に高周波磁場、傾斜磁場を印加し、核磁気共鳴により被検体から発生する信号を計測し、画像化する医用画像診断装置である。
 MRI装置では、プロトン(水素原子核)の密度や緩和時間(T1、T2)を反映した画像(絶対値を画素値とする画像)の他に、静磁場不均一や生体組織間の磁化率差などに起因した磁場変化を反映した画像(位相を画素値とする画像:位相画像)が得られる。
 近年、位相画像が組織間の磁化率差を反映することを利用し、位相画像から生体内の磁化率分布を推定する定量的磁化率マッピング(QSM:Quantitatively Susceptibility Mapping)法が提案されており、さまざまな手法が存在する(例えば、非特許文献1、2)。QSM法により推定された磁化率分布(以後、磁化率画像とよぶ)から大脳基底核の鉄濃度や脳静脈の酸素摂取率を算出することができる。そのためQSM法は、脳疾患である神経変性疾患や脳血管疾患の早期診断に役立つことが期待されている。一方、QSM法は、体幹部への適用も検討されている。代表的な例として、乳房の石灰化検出や肝臓の鉄濃度推定への応用などがある(非特許文献3)。また、前立腺においては、低酸素領域や石灰化検出への応用が期待される。
Sun H他、Quantitative susceptibility mapping using a superposed dipole inversion method: Application to intracranial hemorrhage、Magnetic Resonance in Medicine Sharma S他、Quantitative susceptibility mapping in the abdomen as an imaging biomarker of hepatic iron overload、Magnetic Resonance in Medicine、74巻、3号、673-683頁(2015年) Rochefort L他、Quantitative susceptibility map reconstruction from MR phase data using baysian regularization: Validation and application to brain imaging、Magnetic Resonance in Medicine、63巻、1号、194-206頁(2010年)
 前立腺は、下腹部に存在する男性特有の器官であり、臀部などの皮下脂肪に隣接している。この皮下脂肪の存在により、前立腺磁化率画像の算出精度が低下する可能性がある。QSM法において算出精度が低下する要因の一つに、対象組織の周辺領域におけるMRI信号の欠損がある。前立腺近辺の皮下脂肪では、体外や骨など信号が得られない領域や低信号領域が周辺にあることから、信号欠損による精度低下が生じやすい。また、磁化率の算出精度低下は隣接する領域にも伝搬するため、皮下脂肪に隣接した前立腺などの水領域で算出精度が低下する。
 非特許文献1には、脳内における出血領域とそれ以外の領域の磁化率を別々に算出し、統合する手法が開示されている。また非特許文献2には、脂肪領域を強く平滑化して肝臓の磁化率を算出する手法が開示されている。いずれの手法でも、磁化率画像の画素値を平滑化する制約のもと磁化率を算出している。しかしながら、平滑化によりノイズやアーチファクトを低減することはできるが、信号欠損による精度低下を低減することは困難であると考えられる。また非特許文献3には、脳外の領域を背景とみなして背景磁場除去を行った後に、脳外の領域における磁化率を0とする制約のもと脳内の磁化率を算出する手法が開示されている。本手法では、脳外の領域における算出磁化率に制約を加えることにより、高精度に磁化率画像を算出することができる。しかしながら、本手法を皮下脂肪の算出に適用した場合、高信号領域である水領域も背景とみなして背景磁場除去を行うため、水領域の磁場情報を用いて磁化率を算出できず精度が低下する。
 本発明は、上記事情に鑑みてなされたもので、信号欠損による精度低下を低減して磁化率を高精度に算出することを目的とする。
 本発明は、空間的磁場不均一を反映した核磁気共鳴信号(エコー)を計測し、エコーから複素画像を算出し、磁場画像を算出する。また複素画像から算出した絶対値画像から3つ以上のマスク、低信号領域をあらわす低信号領域マスク、高信号領域であって、判別用画像の所定の閾値で分割された複数のマスク(例えば第一の高信号領域マスク、第二の高信号領域マスク)を算出する。そして、低信号領域マスクを含む複数のマスクを用いて、高信号領域であって磁化率が異なる複数の領域について、それぞれ、磁化率を算出する。その際、一つの領域の磁化率について、他の領域の磁化率を予め定めた特定値とする制約(特定値制約)のもと、磁化率画像を算出する。
 本発明のMRI装置は、具体的には、静磁場磁石と、前記静磁場磁石が形成する空間に傾斜磁場及び高周波磁場を印加する磁場発生部と、前記空間に置かれた被検体から発生する核磁気共鳴信号を受信する受信部と、前記核磁気共鳴信号に対し演算を行う計算機と、を備え、前記計算機は、前記磁場発生部及び前記受信部の動作を制御し、磁場不均一の影響を含む核磁気共鳴信号を計測する計測部と、前記核磁気共鳴信号を用いて画像再構成を行い、複素画像を作成する画像再構成部と、前記複素画像を用いて前記被検体の磁化率画像を算出する画像変換部と、を備える。
 前記画像変換部は、前記複素画像から磁場反映画像を算出する磁場算出部と、前記画像再構成部が作成した画像を用いて、低信号領域に対応する低信号領域マスク及び高信号領域内の複数の領域に対応する複数の領域マスクを算出するマスク算出部と、前記磁場算出部が算出した磁場反映画像と前記マスク算出部が算出した前記低信号領域マスク及び複数の領域マスクとを用いて、前記複数の領域のそれぞれについて磁化率を算出する磁化率算出部と、を備え、前記磁化率算出部は、前記低信号領域を背景とする制約、及び、前記複数の領域のうち特定の領域の磁化率を特定値とする制約のもと、前記特定の領域と異なる領域の磁化率を算出する。なお磁場反映画像とは、周波数或いは磁場を画素値とする画像であって空間的磁場不均一を反映した画像である。
 信号欠損による精度低下を低減して磁化率を高精度に算出することができる。
(a)は垂直磁場方式の磁気共鳴イメージング装置の外観図、(b)水平磁場方式の磁気共鳴イメージング装置の外観図、(c)開放感を高めた磁気共鳴イメージング装置の外観図である。 第一実施形態におけるMRI装置の概略構成を示すブロック図である。 第一実施形態における計算機の機能ブロック図である。 第一実施形態における撮像処理のフローチャートである。 第一実施形態におけるグラディエントエコー型パルスシーケンスの一例を説明するための説明図である。 第一実施形態における画像変換処理のフローチャートである。 第一実施形態におけるマスク算出及び磁化率算出に用いる画像とマスクとの関係を示す図である。 第一実施形態における水脂肪分離処理のフローチャートである。 第一実施形態におけるマスク算出処理のフローチャートである。 第一実施形態における磁化率算出処理のフローチャートである。 第二実施形態の計算機における画像変換部の機能ブロック図である。 第二実施形態の画像変換処理のフローチャートである。 第二実施形態におけるマスク算出処理のフローチャートである。 第三実施形態の計算機における画像変換部の機能ブロック図である。 第三実施形態の画像変換処理のフローチャートである。 第三実施形態におけるマスク算出処理のフローチャートである。 第三実施形態におけるマスク算出及び磁化率算出に用いる画像とマスクとの関係を示す図である。 第三実施形態における磁化率算出処理のフローチャートである。 (a)~(e)は、実施例の結果を示す図である。
 以下、本発明を適用する実施形態について、図面を参照し説明する。以下、本発明の実施形態を説明するための全図において、同一機能を有するものは同一符号を付し、その繰り返しの説明は省略する。また、以下の記述により本発明が限定されるものではない。
 まず、本発明が適用されるMRI装置の概要について説明する。図1は、MRI装置の外観図である。図1(a)は、開放感を高めるために磁石を上下に分離したハンバーガー型(オープン型)の垂直磁場方式のMRI装置(垂直磁場MRI装置)100である。図1(b)は、ソレノイドコイルで静磁場を生成するトンネル型磁石を用いた水平磁場方式のMRI装置(水平磁場MRI装置)101である。また、図1(c)は、図1(b)と同じトンネル型磁石を用い、磁石の奥行を短くし且つ斜めに傾けることによって、開放感を高めたMRI装置102である。なお、図1に示す各MRI装置の形態は、それぞれ、垂直磁場方式、水平磁場方式の一例であり、本発明はこれらに限定されるものではない。以下、水平磁場MRI装置101を例にして説明する。
 図2は、MRI装置101の概略構成を示すブロック図である。このMRI装置101は、被検体に平行な方向に静磁場を発生するマグネット201と、傾斜磁場を発生する傾斜磁場コイル202と、シーケンサ204と、傾斜磁場電源205と、高周波磁場発生器206と、高周波磁場を照射するとともに核磁気共鳴信号(エコー)を検出するプローブ207と、受信器208と、計算機209と、表示装置210と、記憶装置211と、を備える。被検体(例えば、生体)203は寝台(テーブル)等に載置され、マグネット201の発生する静磁場空間内に配される。
 シーケンサ204は、傾斜磁場電源205と高周波磁場発生器206とに命令を送り、それぞれ傾斜磁場および高周波磁場を発生させる。発生された高周波磁場は、プローブ207を通じて被検体203に印加される。被検体203から発生したエコーはプローブ207によって受波され、受信器208で検波が行われる。
 検波の基準とする核磁気共鳴周波数(検波基準周波数f0)は、シーケンサ204によりセットされる。検波された信号は、計算機209に送られ、ここで画像再構成などの信号処理が行われる。その結果は、表示装置210に表示される。必要に応じて、記憶装置211に検波された信号や測定条件、信号処理後の画像情報などを記憶させてもよい。シーケンサ204は、予めプログラムされたタイミング、強度で各部が動作するように制御を行う。プログラムのうち、特に、高周波磁場、傾斜磁場、信号受信のタイミングや強度を記述したものはパルスシーケンスと呼ばれる。
 パルスシーケンスは、目的に応じて種々のものが知られているが、本実施形態では位相画像から得られる空間的磁場変化に基いて磁化率を算出することが目的であるため、磁場強度の空間分布の不均一性を反映したエコー、即ち、空間的磁場変化により各スピンの位相がずれるようなエコーを少なくとも一つ取得するパルスシーケンスを用いる。このようなパルスシーケンスとして、例えば、GrE(Gradient Echo)系のパルスシーケンスを用いることができる。GrE系のパルスシーケンスには、例えば、RSSG(RF-spoiled-Steady-state Acquisition with Rewound Gradient-Echo)シーケンスがある。
 本実施形態の計算機209は、入力装置212を介して設定された(または、あらかじめ設定された)計測パラメータとパルスシーケンスに従いエコーの計測をシーケンサ204に指示し、エコーを計測し、得られたエコーをk空間に配置し、k空間に配置されたエコーに対して演算を行い磁化率画像を算出し、得られた画像を表示装置210に表示する。また、必要に応じて、画像上にROIを設定し、ROI内の画素の統計値を算出する。
 計算機209のこれらの各機能は、記憶装置211に格納されたプログラムを、計算機209のCPUがメモリにロードして実行することにより実現される。なお計算機209の機能の一部或いは全部は、ASIC(Application Specific Integrated Circuit)やFPGA(Field Programmable Gate Array)などのハードウェアによって実現される場合もある。
 本実施形態のプログラムは、具体的には、所定のパルスシーケンスを実行する計測プログラム、計測したエコーを用いて複素画像や絶対値画像を作成する画像再構成プログラム、絶対値画像を用いて3以上の異なる領域マスクを算出するマスク算出プログラム、複素画像から周波数画像或いは磁場画像を算出する磁場算出プログラム、周波数画像或いは磁場画像と3以上の異なる領域マスクを用いて磁化率画像を算出する磁化率算出プログラムを含んでいる。
 これらプログラムと、それを実行することで実現される機能については、以下の実施形態において説明する。
<<第一実施形態>>
 本実施形態は、水領域(水プロトンの信号が支配的である領域)と脂肪領域(脂肪プロトンの信号が支配的である領域)とのそれぞれについて磁化率分布を算出し、統合する。磁化率分布の算出において、平滑化の制約に加え、特定領域について磁化率を特定値とする制約を与える。特定領域は例えば水領域とし、特定値はゼロとする。
 上記処理を行う本実施形態の計算機209の機能ブロック図を図3に示す。図示するように、計算機209は、主たる要素として、計測部300、画像再構成部400、画像変換部500、表示処理部600を備える。画像変換部500は、水脂肪分離処理部510、磁場算出部520、マスク算出部530及び磁化率算出部540を備える。マスク算出部530及び磁化率算出部540の詳細は後述する。
 以下、これら各部による撮像処理の詳細を、図4に示す処理フローに沿って説明する。
[計測(S1101)]
 各種の計測パラメータが設定され、撮像開始の指示を受け付けると、計測部300は、計測を行う(S1101)。ここでは、計測部300は、予め定められたパルスシーケンスに従って、シーケンサ204に指示を行い、エコー信号を取得し、k空間に配置する。シーケンサ204は、指示に従い、上述のように、傾斜磁場電源205と高周波磁場発生器206とに命令を送り、それぞれ傾斜磁場および高周波磁場を発生させる。そして、プローブ207によって受波され、受信器208で検波が行われたエコーを複素信号として受信する。
 このとき、本実施形態では、空間的磁場不均一に対応してスピン毎の位相ずれが生じているエコーを少なくとも一つ計測するパルスシーケンス、例えば、GrE系のパルスシーケンスを用いる。計測部300がステップS1101の計測において実行する計測シーケンスの一例を図5に示す。図5において、RFはRFパルスの、Gsはスライス選択傾斜磁場の、Gpは位相エンコード傾斜磁場の、Grは読み出し傾斜磁場の、それぞれ印加タイミングをそれぞれ示す。またエコーは、エコー信号の取得タイミングを示す。
 計測シーケンス710では、1回の繰り返し時間TR内に以下の手順でエコー信号の計測を行う。本実施形態では、水と脂肪の周波数差を利用して水画像と脂肪画像を算出するため、一つのTR内で2つ以上のエコーを取得する。ここでは、4つの異なるエコー時間でエコー信号を取得する場合を例示する。なお、最初のエコー時間をt、その後のエコー時間の間隔(エコー間隔)をΔtとする。
 まず、RFパルス711を照射し、被検体203の水素原子核スピンを励起する。この際、被検体203の特定のスライスを選択するためにスライス選択傾斜磁場(Gs)712をRFパルス711と同時に印加する。続いてエコー信号に位相エンコードするための位相エンコード傾斜磁場(Gp)713を印加する。その後、最初のRFパルス711照射から時間t後に、読み出し傾斜磁場(Gr)721を印加してエコー信号(第一エコー信号)731を計測する。更に、第一エコー信号731の計測から時間Δt後の時刻tに、極性の反転した読み出し傾斜磁場(Gr)722を印加してエコー信号(第二エコー信号)732を計測する。同様に、第二エコー信号732の計測から時間Δt後の時刻tに、極性の反転した読み出し傾斜磁場(Gr)723を印加してエコー信号(第三エコー信号)733を計測する。さらに、第三エコー信号733の計測から時間Δt後の時刻tに、極性の反転した読み出し傾斜磁場(Gr)724を印加してエコー信号(第四エコー信号)734を計測する。
 第一エコー信号731、第二エコー信号732、第三エコー信号733、第四エコー信号734を計測する際のエコー時間t、t、t、tの少なくとも一つは、水と脂肪の位相差が0にならないような時刻となるように、エコー時間tおよびエコー間隔Δtを設定する。なお、水と脂肪の周波数差をfwfとしたとき、水と脂肪が同位相になる時間をtInとすると、tInはm/fwfである。なお、mは整数である。
 計測シーケンス710では、上述の条件を満たすエコー時間t、t、t、t、またはエコー時間tおよびエコー間隔Δtが選択される。本実施形態では、エコー時間、エコー間隔、エコー取得回数は、ユーザにより入力装置212を介して設定される。あるいは、予め設定される。
 計測部300は、計測シーケンス710を、位相エンコード傾斜磁場713の強度を変化させながら、被検体203の予め定めた撮像領域へのRFパルス711の照射、および同領域からのエコー信号731、732、733、734の計測を、所定回数繰り返す。繰り返し回数は、例えば128回、256回等である。これにより、当該撮像領域の画像再構成に必要な数のエコー信号を繰り返し取得する。繰り返し回数分の第一エコー信号731により、1つの原画像(第一原画像)が形成され、繰り返し回数分の第二エコー信号732、第三エコー信号733、第四エコー信号734により、それぞれ、第二原画像、第三原画像、第四原画像が形成される。これらは、後述する水画像および脂肪画像を算出するための演算用の原画像として記憶装置等に保存され、使用される。
 なお、異なるエコー時間の数、すなわち原画像の数は4つに限らず、任意である。また、k空間において回転状にデータを取得するラジアルスキャンなど、ノンカーテシアン撮像を用いてもよい。
[画像再構成(S1102)]
 計測を終えると、画像再構成部400は、計測した各エコー時間t、t、t、tのエコー信号から画像を再構成する画像再構成処理を行う(ステップS1102)。ここでは各エコー信号を、k空間上に各々配置し、フーリエ変換する。これにより、各エコー時間t、t、t、tの原画像I(第一原画像)、原画像I(第二原画像)、原画像I(第三原画像)、原画像I(第四原画像)をそれぞれ算出する。ここで、各原画像は各画素値が複素数となる複素画像である。
[画像変換(S1103)]
 画像変換部500は、得られた複素画像に対して後述する種々の画像変換処理を行う(ステップS1103)。本実施形態では、画像変換部500は、画像再構成部400で得られた複素画像を水画像・脂肪画像・周波数画像に変換する水脂肪分離処理部510、水脂肪分離処理部510で算出した周波数画像を磁場画像に変換する磁場算出部520、複素画像および水・脂肪画像から低信号領域マスク、脂肪マスク、水マスクの三つのマスクを算出するマスク算出部530、磁場画像および三つのマスクから磁化率画像を算出する磁化率算出部540を有する。本実施形態の画像変換処理の詳細は後述する。
[表示(S1104)]
 表示処理部600は、得られた磁化率画像を濃淡画像として表示装置210に表示する(ステップS1104)。なお、磁化率画像は、最大値投影処理や最小値投影処理などの方法を用いて複数の空間的に連続する画像情報を統合させて表示してもよい。また、磁化率画像に画像処理を行い、磁化率画像と異なるコントラストの画像を作成し、表示装置210に表示させてもよい。例えば、磁化率画像から組織間の磁化率差を強調した強調マスクを作成し、それを絶対値画像にかけ合わせた磁化率差強調画像を表示してもよい。磁化率差強調画像は、磁化率差を強調する処理によって組織の磁化率の情報は失われるが、磁化率の高い組織とそれ以外の組織とのコントラストが増加する。そのため、磁化率の高い組織が明瞭に描出される。
 次に、本実施形態における画像変換部500の処理(S1103)の詳細を説明する。本実施形態の処理(S1103)の説明に先立って、まず一般的な磁化率算出方法の概要を説明する。
 一般的な磁化率算出法では、まず体内外の磁化率差などに起因する大域的な磁場変化を除き、生体組織間の磁化率差などに起因する局所磁場を算出する背景磁場除去処理を行う。その後、磁場変化と磁化率分布の関係式に基き磁化率を算出する。
 背景磁場除去法の代表的な例としてSHARP(sophisticated harmonic artifact reduction for phase data)法やPDF(projection onto dipole field)法などがある。SHARP法やPDF法では、関心領域(ここでは磁化率を算出すべき領域)と背景領域の二つに領域を分割し、背景除去処理を行う。
 背景除去処理後に、以下の式で表わされる相対的な磁場変化δと磁化率分布χの関係式(式(1))を用いて磁化率を算出する。
Figure JPOXMLDOC01-appb-M000001
式中、αは、ベクトル(r’-r)と静磁場方向とのなす角度、d(r)は、点ダイポール磁場をそれぞれ表す。rは画素の座標に対応する位置ベクトルを表す。
 式(1)に示すように、磁場分布δ(r)は、磁化率分布χ(r)と点ダイポール磁場d(r)との畳み込み積分で表される。したがって、式(1)の両辺をフーリエ変換することにより、式(1)は以下の式(2)に変換される。
Figure JPOXMLDOC01-appb-M000002
ここで、k=(k、k、k)は、k空間上の位置ベクトル、Δ(k)、X(k)、D(k)、は、磁場分布δ(r)、磁化率分布χ(r)、点ダイポール磁場d(r)のフーリエ成分をそれぞれ表す。
 式(2)に示すように、磁化率分布のフーリエ成分X(k)は、磁場分布のフーリエ成分Δ(k)を点ダイポール磁場のフーリエ成分D(k)で除算することによって算出できる。しかしながら、式(2)は、D(k)=0近傍の領域において、その逆数が発散してしまうため、直接的にX(k)を算出することができない。
 このD(k)=0となる領域はマジックアングルと呼ばれ、磁場方向に対しておよそ54.7°の2倍の頂角を持つ逆双円錐領域となる。マジックアングルの存在により磁場分布から磁化率分布を推定するQSM法は、不良条件逆問題(ill-conditioned inverse problem)に帰着され、いくつかの解法が提案されている。
 例えば、不良条件逆問題に対しては、正則化とよばれる制約項を用いた方法がよく用いられる(Bertero M and Boccacci P、Introduction to inverse problems in imaging、IOP Publishing、1998)。
 正規化制約項を用いた方法では、たとえば以下の式(3)に示す誤差関数e(χ)を用い、これを最小化する磁化率画像を求める。
Figure JPOXMLDOC01-appb-M000003
ここで、δは磁場画像の列ベクトル、χは磁化率画像候補の列ベクトル、Cはχに対する畳み込み演算に相当する行列、Wは重み、λは任意の定数、Gは微分演算子をあらわす。式(3)において、第二項が正則化項であり、磁化率値を平滑化する制約をあらわし、ノイズやアーチファクトを低減する。
 以上が一般的な磁化率算出法である。
 ところで、前立腺近辺の皮下脂肪では、体外や大腿骨など信号が得られない領域が周辺にあることから、信号欠損による精度低下が顕著になることが予想される。また、磁化率の算出精度低下は隣接する領域にも伝搬するため、皮下脂肪周辺の水領域で算出精度が低下するおそれがある。このとき、算出精度の低下は皮下脂肪周辺の水領域における磁化率値のムラとして生じる。たとえば式(3)のように平滑化の制約のもと磁化率を算出した場合、ノイズやアーチファクトは低減するが、信号欠損による脂肪領域の精度低下や水領域のムラを低減することは困難である。
 本実施形態では、前立腺磁化率マップにおける皮下脂肪領域およびその周辺領域における精度低下を低減するため、脂肪の算出磁化率低下が周囲に伝搬しないよう脂肪領域と水領域の磁化率を別々に算出する。具体的には、水脂肪分離処理により脂肪含有率画像を算出し、脂肪含有率画像を用いて脂肪領域と水領域に領域を分割し、脂肪領域と水領域の磁化率を別々に算出し、最後にそれらを統合する。水領域では、空気領域と脂肪領域を背景とみなした背景磁場除去を行い、磁化率を算出する。また、脂肪領域では、算出精度を改善するため、水領域に磁化率のムラが生じない制約のもとで磁化率を算出する。
 本実施形態の画像変換処理S1103を図6及び図7を参照して詳述する。図6は処理フローを示す図、図7は作成される画像とマスクの関係を示す図である。
 図6に示されるように、本実施形態における画像変換処理は、水脂肪分離処理S1201、磁場算出処理S1203、マスク画像算出処理S1203、磁化率画像算出処理S1204の4ステップからなる。
[水脂肪分離(S1201)]
 まず、本実施形態における水脂肪分離処理S1201の流れについて説明する。本実施形態の水脂肪分離処理部510は、水、脂肪、R (水と脂肪で共通化した見かけの横磁化緩和速度)、及びオフセット周波数分布であらわされる信号モデルに、計測により得られた信号値をフィッティングすることにより、水脂肪分離処理により、水画像・脂肪画像・周波数画像を算出する(図7の点線四角内の処理)。オフセット周波数は、静磁場不均一などによって空間的に変化する共鳴周波数のオフセット量であり、周波数画像はこのオフセット周波数を画素値とする画像である。また「水画像」及び「脂肪画像」は、それぞれ、水プロトンからの信号が支配的である画像、脂肪プロトンからの信号が支配的である画像を意味する。
 図8は、本実施形態の水脂肪分離処理のフローチャートである。まず、エコー時間による信号値の変化を表す信号モデルを設定する(ステップS1301)。次に、後述するフィッティング処理のための各種の初期値を設定する(ステップS1302)。そして、N個(Nは3以上の整数)の異なるエコー時間の原画像Iと初期値と信号モデルとを用いて、非線形最小二乗法により計測信号を信号モデルにフィッティングさせるフィッティング処理を行い(ステップS1303)、水・脂肪・周波数画像を算出する。最後に、得られた水画像と脂肪画像から脂肪含有率画像を算出する(ステップS1304)。以下、各処理の詳細を説明する。
 まず、ステップS1301の信号モデルの設定について説明する。本実施形態では、水の信号強度、脂肪の信号強度、オフセット周波数、見かけの緩和速度R をフィッティング変数とする信号モデルを設定する。ここで、n番目のエコー時間tで得たエコー信号から再構成した画像Iの、任意の画素における信号強度s(n=1、…、N)を表す信号モデルs’を、式(4)に示す。
Figure JPOXMLDOC01-appb-M000004
ここで、tは、n番目のエコー時間、Ψはオフセット周波数、ρおよびρは、水と脂肪の複素信号強度、Kは、時刻tにおける脂肪信号の位相変化量(複素数)、R は、水と脂肪で共通化した見かけの横磁化緩和速度をそれぞれ表す。なお、ステップS1301における信号モデルは、式(4)の形に限らない。たとえば、得られるエコーの数が少ない場合などには、水の信号強度、脂肪の信号強度、オフセット周波数のみをフィッティング変数とする信号モデルを設定してもよい。
 なお、脂肪信号は、その分子構造により、複数のスペクトルピークを持つことが知られている。したがって、P個(Pは1以上の整数)のピークを持つ脂肪信号を考えた場合、脂肪信号の位相変化量Kは、以下の式(5)で表される。
Figure JPOXMLDOC01-appb-M000005
およびfは、p番目(pは、1≦p≦Pを満たす整数)の脂肪ピークの相対信号強度および水との周波数差をそれぞれ表す。なお、aは、以下の式(6)の条件を満たす。
Figure JPOXMLDOC01-appb-M000006
 以下、本実施形態では、脂肪が6つのピークを持つ(P=6)信号モデルを用いる。
 次に、ステップS1302の、初期値の設定方法について説明する。設定する初期値は、各画素の水の複素信号強度および脂肪の複素信号強度と、オフセット周波数分布と、見かけの横磁化緩和速度R とである。水の複素信号強度ρおよび脂肪の複素信号強度ρの各画素の初期値は、実際に計測して得た原画像の信号値sの絶対値|s|を時間方向に最大値投影した値|smaxを用いる。オフセット周波数分布の初期値は、全画素で0とする。見かけの横磁化緩和速度R の初期値は、全画素で1とする。
 なお、水と脂肪の複素信号強度ρおよびρの各画素の初期値、オフセット周波数分布の初期値、見かけの横磁化緩和速度R の初期値は必ずしも上記の値でなくてもよく、非線形最小二乗法による演算結果が発散または振動しない値であればよい。
 次に、ステップS1303のフィッティング処理について説明する。本実施形態では、フィッティングにより真値を推定する変数を、各画素の、水の信号強度をρ、脂肪の信号強度をρ、見かけの横磁化緩和速度をR 、オフセット周波数をΨとする。そして、各変数の推定値をそれぞれ、ρ’、ρ’、R ’、Ψ’とし、真値と推定値との差分をそれぞれ、Δρ、Δρ、ΔR 、ΔΨとする。このとき、実際の計測で得た信号値をs、推定値ρ’、ρ’、R ’、Ψ’を式(4)で表される上記信号モデルに代入して得た推定信号をs’とすると、計測信号sと推定信号s’との差分Δsは、行列表記により、以下の式(7)で表される。
Figure JPOXMLDOC01-appb-M000007
従って、ベクトルxは、行列Aの擬似逆行列を用いて、以下の式(8)により算出できる。
Figure JPOXMLDOC01-appb-M000008
なお、Aは、Aの複素共役転置行列を表す。
 本実施形態におけるフィッティング処理S1303では、式(8)で算出した差分ベクトルxの各要素であるΔρ、Δρ、ΔR 、ΔΨを、推定値ρ’、ρ’、R ’、Ψ’にそれぞれに加えて、推定信号s’を再計算した後、再び式(8)を用いて差分ベクトルxを算出する。この手順を繰り返すことにより、差分ベクトルxを最小化し、推定値を真値へ近づけていく。本実施形態におけるフィッティング処理では、任意の閾値を設定し、差分ベクトルxのノルムがその閾値以下となるまで、上記手順を繰り返す。そして、最終的に得られた、各画素の水の信号強度ρ’脂肪の信号強度ρ’、見かけの横磁化緩和速度R ’、オフセット周波数Ψ’を、それぞれ、水画像、脂肪画像、R 画像、周波数画像とする。なお、非線形最小二乗法には、レーベンバーグ・マーカート法など任意の公知の方法を用いることもできる。
 次に、フィッティング処理S1303で得た水画像と脂肪画像から脂肪含有率画像を算出する(ステップS1304)。脂肪含有率画像は、各画素の画素値が脂肪含有率となる画像である。本実施形態における脂肪含有率の各画素値は、脂肪画像の画素値を、水画像と脂肪画像の画素値の和で除算することにより算出される。脂肪含有率画像は、後述するマスク算出処理において、脂肪マスクを算出する際の判定用画像として用いられる。
 なお、上述した水脂肪分離処理S1201では、エコー時間の異なるエコーを4つ用いて、4つの未知数である水画像、脂肪画像、周波数画像及びR 画像を得たが、R 画像を得ることは必須ではない。従って、本処理はエコー時間の異なるエコーの数Neが3以上の場合に適用できる。
 さらに周波数画像については、エコーの数が1つの場合でも、原画像Iの各画素値I(r)とエコー時間tを用いて以下の式(9)により周波数画像の各画素値Ψ’を算出できる。
Figure JPOXMLDOC01-appb-M000009
ここで、arg(c)は、複素数cの偏角をあらわす。
また、エコーの数が2つの場合は、第一原画像Iおよび第二原画像Iのそれぞれの画素値(I(r)およびI(r))とエコー時間(tおよびt)を用いて以下の式(10)により周波数画像を算出できる。
Figure JPOXMLDOC01-appb-M000010
 従って、上述したフィッティングによる水脂肪分離処理とは別に、周波数画像を取得する処理を行ってもよい。
[磁場算出処理(S1202)]
 次に、算出した周波数画像Ψ’に対し、以下に述べる磁場算出処理を行い、磁場画像δを算出する(ステップS1202)。磁場画像は、静磁場強度で規格化した相対的な磁場変化をあらわす画像である。磁場算出部520は、周波数画像Ψ’(r)を以下の式(11)により変換することにより磁場画像δ(r)を算出する。
Figure JPOXMLDOC01-appb-M000011
ここで、γはプロトンの磁気回転比、Bは静磁場強度である。
[マスク算出処理(S1203)]
 次に、本実施形態におけるマスク算出処理S1203の流れについて説明する。本実施形態のマスク算出部530は、画像再構成部400で算出した複素画像および水脂肪分離処理部510で算出した水画像・脂肪画像から、空気などの低信号領域をあらわす低信号領域マスク、ボクセル内で水が支配的となる領域をあらわす水マスク、ボクセル内で脂肪が支配的となる領域をあらわす脂肪マスクを算出する。このためマスク算出部530は、図3に示すように、低信号領域マスク算出部531、脂肪マスク算出部532、水マスク算出部533を備える。
 図9に、本実施形態のマスク算出部530の処理フローチャートを示す。本実施形態では、まず、複素画像からマスク用絶対値画像を算出し、マスク用絶対値画像から低信号領域マスクおよび高信号領域マスクを算出する(ステップS1401)。次に、高信号領域マスクおよび上記ステップS1304で算出した脂肪含有率画像から脂肪マスクを算出する(ステップS1402)。最後に、高信号領域マスクおよび脂肪マスクから水マスクを算出する(ステップS1403)。以下、各処理の詳細を説明する。
 まず、ステップS1401の低信号領域マスク算出方法(図7の一点鎖線の四角内の処理)について説明する。まずマスク用絶対値画像を算出する。マスク用絶対値画像は、各画素が0以上の実数であらわされる画像である。本実施形態におけるマスク用絶対値画像は、全エコーの絶対値画像から画素ごとに時間方向の二乗和平方根を計算した1枚の絶対値画像とする。なお、マスク用絶対値画像は、任意の方法で算出される。たとえば、1エコー目の絶対値画像としてもよい。
 次に、マスク用絶対値画像から低信号領域マスクを算出する。本実施形態では、マスク用絶対値画像から閾値処理により低信号領域マスクMaを算出する。低信号領域マスクは、予め定めた閾値を用い、マスク用絶対値画像の画素値が当該閾値より小さい領域を1、大きい領域を0とする。なお、閾値は、絶対値画像の全画素の画素値分布から判別分析法等の方法を用いて求めてもよい。また、閾値処理を行った後に、孤立点除去処理などの画像処理を行ってもよい。なお、低信号領域マスク算出法には様々な方法がある。例えば、体内外の境界部分を検出し、検出結果をもとに体内領域を0、体外領域を1としてもよい。ノイズマスク処理では、こうした他の方法を用いてもよい。
 また、画素ごとに1から低信号領域マスクの画素値を減算したものを高信号領域マスクMhとして算出する。
 次に、ステップS1402の脂肪マスク算出方法について説明する。本実施形態では、高信号領域マスクMhを脂肪マスクMfと水マスクMwに分割するために、両者を判別するための画像(判別用画像)を用意し、判別用マスクを作成する。ここでは、判別用画像として、水脂肪分離処理S1201(S1304)で算出した脂肪含有率画像を用い、任意の閾値を用いて閾値処理することにより判別用マスクMdを算出する。例えば、脂肪含有率画像の、脂肪含有率が当該閾値より大きい領域を1、脂肪含有率が当該閾値より小さい領域を0とする。本実施形態では、閾値は0.8に設定する。このように算出した判別用マスクMdと高信号領域マスクMhをかけあわせることにより、脂肪マスクMfを算出する。
 なお、脂肪マスクの算出方法は任意であり、たとえば脂肪画像を閾値処理する等の方法で算出してもよい。この場合には、水脂肪分離処理(図6、S1201)における脂肪含有率画像算出処理(図8、S1304)は省略することができる。また、閾値処理後に、孤立点を除去する処理や、任意の閾値以下の脂肪領域を0とする処理を加えて脂肪マスクとしてもよい。この処理を加えることより、ノイズにより脂肪領域と判別された領域を除去することができる。
 最後に、水マスクMwを算出する(ステップS1403)。本実施形態では、高信号領域マスクMhから脂肪マスクMfを減算することにより水マスクMwを算出する。なお、水マスクの算出方法は任意であり、たとえば水画像を閾値処理する等の方法で算出してもよい。
 また本実施形態では、低信号領域マスクMa、脂肪マスクMf、水マスクMwの画素値は0または1としたが、任意の値とすることもできる。たとえば、脂肪マスクは脂肪含有率画像に高信号領域マスクをかけあわせ、水マスクは画素ごとに1から脂肪含有率を減算した画像に高信号領域マスクをかけあわせることにより算出してもよい。
[磁化率算出処理S1204]
 次に、本実施形態における磁化率算出処理(図7の二点鎖線の四角内の処理)について説明する。本実施形態の磁化率算出部540は、図3に示すように、マスク算出部530で算出した低信号領域マスク、水マスク、脂肪マスクと、磁場算出部520で算出した磁場画像とから、磁化率画像を算出する。このため磁化率算出部540は、図3に示すように、背景磁場除去部541、脂肪磁化率算出部542、背景磁場除去部543、水磁化率算出部544、統合磁化率算出部545を備える。
 本実施形態の磁化率算出部540の処理S1204の詳細を図10に示すフローチャートを参照して説明する。図10の各ステップS1501~S1505は、それぞれ、磁化率算出部540の各部541~545の処理に対応する。本実施形態では、まず、磁場画像に対し、低信号領域マスクで定義される領域を背景領域とみなした背景除去処理を行う(ステップS1501)。次に、S1501で背景除去した画像から磁化率算出処理を行い、水領域が0となる制約のもと脂肪磁化率画像を算出する(ステップS1502)。次に、磁場画像に対し、低信号領域マスクと脂肪マスクで定義される領域を背景領域とみなした背景除去処理を行う(ステップS1503)。次に、S1503で背景除去した画像から磁化率算出処理を行い、水磁化率画像を算出する(ステップS1504)。最後に、脂肪磁化率画像と水磁化率画像を統合することにより統合磁化率画像を算出する(ステップS1505)。以下、各処理の詳細を説明する。
 まず、ステップS1501で、磁場算出処理S1202で算出した磁場画像δに対し、低信号領域マスクで定義される領域を背景領域とみなした背景磁場除去を行う。背景磁場除去は任意の方法を用いて行うことができるが、本実施形態ではSHARP法を用いる。
 具体的には、背景領域から生じる磁場が関心領域内(マスク算出処理S1203で算出した高信号領域マスクで定義される領域内)で調和関数の性質を満たすことを利用し、関心領域内における下記の式を打ち切り特異値分解法や最小二乗法などにより解くことで局所磁場を算出する。
[数12]    MHδlocal = MHδtotal  (12)
 ここで、Mは関心領域を1とする二値マスク(高信号マスク)を対角成分にもつ直交行列、Hは畳込み演算をあらわす行列、δlocalは局所磁場画像(背景除去後の磁場δ’)の列ベクトル、δtotalは背景除去前の磁場画像(磁場算出処理S1202で算出した磁場画像δ)の列ベクトルである。ここでは、関心領域はMの対角成分で定義される領域であり、背景領域は、Iを単位行列としたとき、I-Mの対角成分で定義される領域である。
 次に、ステップS1502では、S1501で背景磁場除去した磁場画像から、磁化率を平滑化する制約をおく平滑制約と特定の領域の磁化率を特定値とする制約をおく特定値制約のもと、脂肪磁化率画像χfを算出する。特定の領域とは、ここでは水マスクで定義される領域であり、特定値を0とする。
 以下、平滑制約と特定値制約を含む、磁化率画像χfの算出手法の具体例を説明する。
 具体例の一つは、式(13)で表されるような、磁化率を平滑化する制約をおく平滑制約項と特定値領域の磁化率を特定値とする制約をおく特定値制約項とを有する目的関数を最小化することにより脂肪磁化率画像χfを算出する。
Figure JPOXMLDOC01-appb-M000012
またCは畳込み演算をあらわす演算子、Gは微分演算子を表す。また、Mは高信号領域マスク、δ’はS1501で背景除去した磁場、Mwは水マスク、Mfは脂肪マスクを表す。式(13)の第一項はマスクM内を関心領域として磁化率と周波数の関係式の誤差を最小化する項、第二項は水領域を0とする制約をおく項、第三項は脂肪内を平滑化する項である。λwおよびλfは、各制約の大きさを規定するもので、任意の定数である。以下、λwを特定値パラメータ、λfを平滑パラメータと呼ぶ。
 特定値制約の大きさを規定する特定値パラメータλwと平滑制約の大きさを規定する平滑パラメータλfは、第二項の精度改善効果と第三項のノイズ低減効果を最大にするという観点から、脂肪磁化率画像の脂肪マスク内における平均値と標準偏差の比が最大になる値で設定することが好ましい。なお、これらのパラメータは任意の方法で決めることができる。例えば、予め定めた任意の値としてもよい。
 ただし、本実施形態では、λwは所定の値(λwth)以下の値に設定する。非特許文献3には、脳外の領域を背景領域とみなした背景磁場除去後に、背景領域である脳外の磁化率を0とする制約のもと脳内の磁化率を算出する手法が開示されている。この手法では、背景磁場除去処理により磁化率の寄与が除かれた背景領域に0とおく制約を加えるため、制約に用いる正則化項の正則化パラメータが1020と大きい。一方、式(13)においては、制約を加える領域である水領域は背景領域ではない。そのため、大きすぎる特定値パラメータにより水領域における算出磁化率の精度低下が生じ、その結果周辺の脂肪領域における磁化率の精度低下が生じる。λwの値を所定の値(λwth)以下とすることにより、精度低下を防ぎ、パラメータ探索のために要する計算時間を低減することができる。なお、λwthは第一項に用いるMやWの値に応じて異なるが、本実施形態では、λwthは例えば10に設定する。
 式(13)は共役勾配法などを用いた繰り返し演算により最小化することができる。その際、初期ベクトルは任意であるが、たとえば0ベクトルを用いる。また、繰り返し演算の終了条件は任意であるが、たとえば繰り返し回数が、あらかじめ定めた任意の閾値をこえたときに繰り返し演算を終了するなどの方法がある。
 式(13)を用いた算出方法では、第二項を用いて水領域にムラが生じない制約のもと磁化率を算出することにより、脂肪磁化率の算出精度を改善することができる。
 なお、式(13)における第一項目では、マスクMのかわりに重みWを用いてもよい。重みWは、通常、背景領域を除いた関心領域内で定義され、たとえば位相画像に基づいて算出される(本発明者らによる特開2015-62637号公報記載の方法)。具体的には、位相画像Pから、各画素の位相ばらつきを算出し、算出した位相ばらつきの大きさに基づいて重みWを算出する。例えば、位相ばらつきが大きいほど、重みWが小さくなるよう算出する。位相ばらつきには、位相ばらつきを算出する対象画素の周辺領域の、複数の画素の画素値(位相値)を用いて算出した、例えば、標準偏差等の統計指標を用いる。
 脂肪磁化率を算出する式としては、式(13)の他に、例えば次式(14)を用いてもよい。
Figure JPOXMLDOC01-appb-M000013
 式(14)では、第一項にMfを用いることにより、水領域の磁化率を0にする特定値制約をおくことができる。第二項は、式(13)の第三項と同様の平滑制約である。なお、第一項のMfのかわりに、水領域が任意の定数c(例えば0.2)、脂肪領域が任意の定数c(例えば0.8)となる重み関数Mf’を用いてもよい(ただしc<c)。
 さらに別の磁化率算出手法として、予め定めた初期磁化率分布χf0と磁場分布δとにより定義した誤差関数に基き解を更新する処理と平滑化処理を交互に行う方法(本発明者らによる特願2014-228843号:国際公開公報2016076076号に記載の方法:以下、「新磁化率算出法」と呼ぶ)に基き、磁化率画像を算出してもよい。
 具体的には、磁場分布δと初期磁化率分布χf0とから、更新磁化率分布χUPDを算出し、更新磁化率分布χUPDを平滑化してフィルタ処理後の磁化率分布(処理後磁化率分布)χSMTを算出し、更新磁化率分布χUPDと処理後磁化率分布χSMTとを合成して合成磁化率分布χADDを算出し、繰り返し演算の終了判定を行うとともに、初期磁化率分布χf0を合成磁化率分布χADDに更新する処理を繰り返し行う。
 更新磁化率分布χUPDは、下式に示す誤差関数e(χf0)を用いて初期磁化率分布χf0から算出する。
Figure JPOXMLDOC01-appb-M000014
ここでは、式(15)から誤差関数e(χf0)の勾配▽e(χf0)を算出し、算出した勾配▽e(χf0)を用いて、以下の式(16)により更新磁化率分布χUPDを算出する。
[数16]
    χUPD=χf0+▽e(χf0) ・・・(16)
 なお、最初に設定する初期磁化率分布χf0は、任意のものを設定できる。例えば、全ての値を0としたものを初期磁化率分布χf0としてもよい。以上の処理では、式(15)においてMfを用いることにより、水領域の磁化率を0にする制約をおくことができる。なお、式(15)におけるMfのかわりに、水領域が任意の定数c(例えば0.2)、脂肪領域が任意の定数c(例えば0.8)となる重み関数Mf’を用いてもよい(ただしc<c)。
 処理後磁化率分布χSMTは、更新磁化率分布χUPDに、エッジ保存平滑化フィルタ(国際公開公報2016076076号記載のフィルタ)を適用して平滑化することにより算出する(平滑制約に相当する処理)。このエッジ保存平滑化フィルタは、局所分散が小さい領域では平滑化が強く働き、局所分散が大きい領域では平滑化が弱く働くという性質を持った平滑化フィルタであり、磁化率の定量値を維持できるとともに微細区像を含むエッジ部分を保存する。なお、平滑化処理は、脂肪領域のみに行ってもよい。
 合成磁化率分布χADDは、更新磁化率分布χUPDと処理後磁化率分布χSMTとを合成することにより算出する。具体的には、まず、更新磁化率分布XUPDおよび処理後磁化率分布XSMTそれぞれを、フーリエ変換によりk空間データに変換する。そして、変換後の各k空間データ(フーリエ成分)に重み付けを行い、重み付けを行った各k空間データを加算する。そして、加算後、逆フーリエ変換により合成磁化率分布χADDを得る。
 k空間における重み付けは、重み画像G(k)(kは、k空間のデータ位置を示す。以下、同様)を作成し、重み画像Gをk空間データに掛け合わせることにより行う。このとき、更新磁化率分布XUPDから得たk空間データに乗算する重み画像GUPDでは、例えば、マジックアングル領域近傍のデータに付与する重みは、他の領域のデータに付与する重みより小さく設定される。また、処理後磁化率分布XSMTから得たk空間データに乗算する重み画像GSMTでは、逆に、マジックアングル領域近傍のデータに付与する重みは、他の領域のデータに付与する重みより大きく設定される。これらの条件を満たす重み画像GUPD(k)、GSMT(k)として、点ダイポール磁場のフーリエ成分D(k)とマジックアングル領域の大きさを調整する調整パラメータbとを用いて、例えば、以下の式(17)より定義する重み画像を用いることができる。
Figure JPOXMLDOC01-appb-M000015
上記式(17)では、更新磁化率分布χUPDにおいて、フーリエ成分の値が所定の値(b)より小さい領域を、マジックアングル領域とする。
 繰り返し演算により誤差関数e(χf0)を最小化する最小化処理において、繰り返し演算の終了判定を行う。処理を継続する、すなわち、繰り返すと判定した場合は、その時点で得られた最新の合成磁化率分布χADDを、初期磁化率分布χf0として設定することにより、初期磁化率分布χf0を更新する。なお、終了と判定した際は、その時点の最新の合成磁化率分布χADDを磁化率分布(ここでは脂肪磁化率分布)として出力する。例えば、以下の式(18)に示す評価関数f(χUPD)を用い、その値が予め定めた閾値εより小さい場合、終了と判定する。
Figure JPOXMLDOC01-appb-M000016
 以上、ステップS1502における脂肪内磁化率算出処理の具体例を説明したが、ステップS1502では、磁化率を平滑化する制約をおく平滑制約と特定の領域の磁化率を特定値とする制約をおく特定値制約のもと、脂肪磁化率画像χfを算出する方法であれば、これら具体例に限定することなく採用することが可能である。
 なお、脂肪磁化率算出S1502では、補間などの方法によりδ’、M、Mf、Mwなどの入力画像(背景除去後の磁場画像、高信号マスク画像、脂肪マスク画像及び水マスク画像)の解像度をおとしてから磁化率画像χfを算出し、算出後に磁化率画像χfの解像度を補間などの方法により元の大きさに戻してもよい。例えば、各入力画像の各方向における画像サイズを半分にしてからχfを算出し、算出後に各方向における画像サイズを元の大きさに戻してもよい。診断領域が脂肪領域になく高い解像度が必要とされない場合に、本処理を行うことにより、χfの算出に要する計算時間を短縮することができる。
 次に、ステップS1503において、磁場算出処理S1202で算出した磁場画像δに対し、低信号領域マスクと脂肪マスクの和で定義される領域を背景領域とみなした背景磁場除去を行う(ステップS1503)。本処理により、水と脂肪の磁化率差による周波数変化(磁場変化)を除去することができる。ステップS1503の背景磁場除去についても、SHARP法やPDF法など任意の方法により背景磁場除去を行うことができる。例えば、SHARP法を用いた場合、本ステップS1503における背景磁場除去は、背景領域として、低信号領域マスクと脂肪マスクの和で定義される領域を用いること以外は、ステップS1501で説明した背景磁場除去と同様である。
 次に、ステップS1504において、S1503で背景磁場除去した磁場画像から水マスクで定義される領域を関心領域とした磁化率算出処理を行い、水磁化率画像χwを算出する。ステップS1504でも脂肪磁化率算出処理のステップS1502と同様に任意の方法により磁化率算出を行うことができるが、本実施形態では新磁化率算出法(国際公開公報2016076076号記載の方法)により磁化率を算出する。このとき、磁場と磁化率との関係式に基づく制約条件は水マスクが1となる領域で定義される。
すなわち、脂肪磁化率算出に用いた上述した式(15)に代えて、次式(19)を用いる。
Figure JPOXMLDOC01-appb-M000017
 その他の算出方法は、脂肪磁化率の算出方法で説明した方法と同様であり、予め定めた初期磁化率分布χw0を更新する処理と平滑化処理とを交互に行う処理を繰り返し、所定回数の繰り返し演算後に得た合成磁化率分布χADDを水磁化率とする。
 最後に、ステップSステップ1505において、S1502で得た脂肪磁化率(脂肪内磁化率画像)とステップS1504で得た水磁化率(脂肪外磁化率画像)とを式(20)により統合し、統合磁化率画像χを算出する(ステップS1505)。
Figure JPOXMLDOC01-appb-M000018
 以上の処理(S1501~S1505)で、磁化率算出処理(図6:S1204)が完了し、図4の画像変換処理ステップS1103が完了する。統合磁化率画像に対し、必要に応じて所定の画像処理を行い、磁化率画像(絶対値画像)、或いは磁化率強調画像を作成し、表示部600により表示すること(ステップS1104)は前述のとおりである。
 以上、説明したように、本実施形態のMRI装置は、静磁場磁石と、前記静磁場磁石が形成する空間に傾斜磁場及び高周波磁場を印加する磁場発生部と、前記空間に置かれた被検体から発生する核磁気共鳴信号を受信する受信部と、前記核磁気共鳴信号に対し演算を行う計算機と、を備え、前記計算機は、前記磁場発生部及び前記受信部の動作を制御し、磁場不均一の影響を含む核磁気共鳴信号を計測する計測部と、前記核磁気共鳴信号を用いて画像再構成を行い、複素画像を作成する画像再構成部と、前記複素画像を用いて前記被検体の磁化率画像を算出する画像変換部と、を備える。画像変換部は、前記複素画像から磁場反映画像を算出する磁場算出部と、前記画像再構成部が作成した画像を用いて、低信号領域に対応する低信号領域マスク及び高信号領域内の複数の領域に対応する複数の領域マスクを算出するマスク算出部と、前記磁場算出部が算出した磁場反映画像と前記マスク算出部が算出した前記低信号領域マスク及び複数の領域マスクとを用いて、前記複数の領域のそれぞれについて磁化率を算出する磁化率算出部と、を備える。磁化率算出部は、前記低信号領域を背景とする制約、及び、前記複数の領域のうち特定の領域の磁化率を特定値とする制約のもと、前記特定の領域と異なる領域の磁化率を算出する。
 また、磁化率算出部は、前記低信号領域及び前記特定の領域と異なる領域を背景とする制約のもと、前記特定の領域の磁化率を算出し、前記特定の領域と異なる領域の磁化率と統合する。
 本実施形態の画像変換部は、前記画像再構成部が作成した画像を用いて、水プロトンからの信号が支配的である水画像、脂肪プロトンからの信号が支配的である脂肪画像を分離する水脂肪分離部をさらに備え、マスク算出部は、前記水画像及び前記脂肪画像を用いて算出した画素毎の脂肪含有率の所定の値を閾値として用い、前記複数の領域マスクとして、水が支配的な領域に対応する水マスク及び脂肪が支配的な領域に対応する脂肪マスクを作成する。
 以上、説明したとおり、本実施形態のMRI装置は、脂肪画像と水画像のそれぞれについて、磁化率分布を算出する際、脂肪領域では、水領域に磁化率のムラが生じない制約、すなわち磁化率が0となる制約のもとで磁化率を算出することにより、信号欠損による精度低下が低減する。水領域では、空気領域と脂肪領域を背景とみなした背景磁場除去を行い磁化率を算出するため、水脂肪の磁化率差による磁場変化が除去され、ムラが生じない。すなわち、脂肪領域および水領域において高精度に磁化率画像を算出することができる。
 なお、本実施形態では水平磁場MRIについて説明したが、垂直磁場MRIやその他の装置を用いても、同様の処理が適用でき、同様の効果が得られる。また、撮像断面も、横断面、冠状断面、矢状断面、オブリーク断面など任意の撮像断面で同様の処理が適用でき、同様の効果が得られる。
 また、本実施形態では、画像再構成部、画像変換部、表示処理部の各部の機能を、MRI装置が備える計算機内で実現する場合を例にあげて説明したが、これに限られない。これらの各部の少なくとも1つは、例えば、MRI装置の計算機209とデータの送受信が可能なMRI装置とは独立した情報処理装置上に構築されていてもよい。
<変更例>
 上記第一実施形態では、磁化率算出部540が、脂肪領域と水領域のそれぞれについて磁化率を算出した場合を説明したが、脂肪の磁化率が既知の場合、脂肪領域の磁化率を予め定めた特定値として算出することもできる。この場合、図3の磁化率算出部540において、脂肪磁化率算出部542は省略でき、また図10に示すフローにおいて背景磁場除去処理S1501及び脂肪磁化率算出処理S1502は行わず、処理内容は次のように変更される。
 磁化率算出部540は、まず、低信号領域マスクで定義される領域を背景領域とみなした背景除去を行ったあと(S1503)、下記の式を用いて磁化率画像χを算出する(S1504)。
Figure JPOXMLDOC01-appb-M000019
ここで、Mは高信号領域マスク、Mfは脂肪マスク、Δχfは任意の定数、eは単位ベクトルをあらわす。Δχfには、予め定めた水と脂肪の磁化率差などを用いる。たとえば、0.61ppmとする。第二項は磁化率を平滑化する制約をあらわす。また、第三項は脂肪領域における磁化率を特定値Δχfとする制約を表す。λ及びλfatは、それぞれ第二項及び第三項の制約の強さをあらわす定数である。
 式(21)は共役勾配法などを用いた繰り返し演算により最小化することができる。その際、初期ベクトルは任意であるが、たとえば0ベクトルを用いる。また、繰り返し演算の終了条件は任意であるが、たとえば繰り返し回数が、あらかじめ定めた任意の閾値をこえたときに繰り返し演算を終了するなどの方法がある。
 なお、脂肪の磁化率を予め定めた特定値として算出する式は、式(21)の形に限られない。たとえば、第二項には、L1ノルムによる制約を用いてもよい。
 本変更例によれば、脂肪の磁化率が既知の場合、第三項の制約を加えることにより、水領域の磁化率を高精度に算出することができる。算出した水領域の磁化率と脂肪の磁化率(特定値)とを式(20)を用いて統合磁化率画像χを算出することは第一実施形態と同じである。
<<第二実施形態>>
 第一実施形態では、脂肪含有率画像を閾値処理することにより脂肪マスクと水マスクの二つの高信号領域マスクを算出し、それぞれの領域で異なる方法で磁化率を算出し、最後に統合することにより磁化率画像を算出した。第二実施形態では、脂肪含有率画像ではなく事前に算出した磁化率画像を閾値処理することにより二つの高信号領域マスクを算出する。そして、第一実施形態と同様の処理を行うことにより磁化率画像を算出する。以下、第一実施形態と異なる点を中心に本実施形態を説明する。
 本実施形態のMRI装置は、基本的に第一実施形態と同様の構成(図2)を有する。
 本実施形態の計算機209は、第一実施形態と同じように、入力装置212を介して設定された(または、あらかじめ設定された)計測パラメータとパルスシーケンスに従いエコーの計測をシーケンサ204に指示し、エコーを計測し、得られたエコーをk空間に配置し、k空間に配置されたエコーに対して演算を行い磁化率画像を算出し、得られた画像を表示装置210に表示する。なお第一実施形態では、水脂肪分離を前提としており、エコー時間の異なる2以上のエコーを計測したが、本実施形態は一つのエコー時間のみで計測を行ってもよい。
 また本実施形態の計算機209も、計測部300、画像再構成部400、画像変換部500、表示処理部600を備えることは、第一実施形態と同様であるが、画像変換部500の構成が異なる。図11に、本実施形態における計画像変換部500の機能ブロック図を示す。図11において、第一実施形態と同様の機能を有する要素は、図3と同じ符号で示し、重複する説明は省略する。
 図11に示すように、本実施形態の画像変換部500は、磁場算出部520、マスク算出用磁化率算出部550、マスク算出部530、磁化率算出部540を備える。なお図11では図示を省略しているが、本実施形態の画像変換部500は、第一実施形態と同様に、画像変換部500が水脂肪分離処理部510を備えていてもよい。マスク算出用磁化率算出部550は、マスク算出部530における処理に必要な磁化率画像を算出する。このマスク算出用磁化率画像は、磁化率算出部540が算出する最終的な磁化率とは異なり、マスクを作るための暫定的な磁化率画像である。マスク算出部530は、低信号領域マスク算出部810、高磁化率マスク算出部820、低磁化率マスク算出部830を備える。磁化率算出部540は、背景磁場除去部841、843、高磁化率(第一の磁化率画像)算出部842、低磁化率(第二の磁化率画像)算出部844、及び統合磁化率算出部845を備える。
 以上の構成を踏まえ、本実施形態における画像変換部500の処理の流れについて、図12及び図13を参照して説明する。図12において、図6と同じ処理は同じ符号で示し詳細な説明は省略する。
 本実施形態の画像変換部500は、図12に示すように、磁場算出処理S1201、マスク算出用磁化率画像算出処理S1210、マスク算出処理S1203、磁化率算出処理S1204を行う。
[磁場算出処理(S1202)]
 磁場算出処理S1201では、第一実施形態と同様に、フィッティイングにより周波数画像を算出し、或いは前述の式(9)、又は(10)を用いて周波数画像を算出し、周波数画像から、周波数と磁場との関係式(式(11))を用いて磁場画像を算出する。なおフィッティングによる周波数画像を算出する場合には、エコー時間の異なるエコーの数は2以上とし、周波数画像の算出とともに水脂肪分離を行ってもよい。
[暫定磁化率算出処理(S1210)]
 次にステップS1210において、画像再構成部400で算出した複素画像を用いて暫定的な磁化率画像を算出する(マスク算出用磁化率算出部550の処理)。マスク算出用磁化率画像は、例えば「新磁化率算出法」(国際公開公報2016076076号記載の方法)により磁化率を算出できる。このとき、磁場と磁化率との関係式に基づく制約条件は高信号領域マスクが1となる領域で定義される。
[マスク算出処理(S1203)]
 ステップS1203では、画像再構成部400で算出した複素画像、磁場算出部520で算出した磁場画像及びマスク算出用磁化率算出部550で算出したマスク算出用磁化率画像を用いて、3つのマスク、即ち、空気などの低信号領域をあらわす低信号領域マスク、相対的に高磁化率な領域をあらわす高磁化率マスク、相対的に低磁化率な領域をあらわす低磁化率マスクを算出する。
 図13は、マスク算出処理のフローチャートの一例を示す。この例では、まず、第一の実施形態と同様の方法で低信号領域マスクと高信号領域マスクを算出する(ステップS1601)。次に、高信号領域マスクおよびマスク算出用磁化率画像から高磁化率マスクを算出する(ステップS1602)。最後に、高信号領域マスクおよび高磁化率マスクから低磁化率マスクを算出する(ステップS1603)。以下、各処理の詳細を説明する。
 ステップS1601の低信号領域マスク算出方法は、基本的に第一の実施形態と同様である。すなわち、複素画像からマスク用絶対値画像を算出し、マスク用絶対値画像から低信号領域マスクおよび高信号領域マスクを算出する。
 次に、ステップS1602では、ステップS1210で算出したマスク用磁化率画像を、高信号領域を分割するための判別用画像とし、任意の閾値を用いて閾値処理することによりマスクを算出する。すなわち、マスク用磁化率画像における算出磁化率が当該閾値より大きい領域を1、算出磁化率が当該閾値より小さい領域を0とする。閾値は任意であるが、本実施形態では、たとえば0.5ppmに設定する。そして、算出したマスクと高信号領域マスクをかけあわせることにより、高磁化率マスクを算出する。
 この際、第一の実施形態と同様、閾値処理後に、孤立点を除去する処理や、面積が任意の閾値以下となる高磁化率領域を0とする処理を加えて高磁化率マスクとしてもよい。
 最後に、低磁化率マスクを算出する(ステップS1603)。本実施形態では、高信号領域マスクから高磁化率マスクを減算することにより低磁化率マスクを算出する。
 なお本実施形態では、低信号領域マスク、高磁化率マスク、低磁化率マスクの画素値は0または1としたが、任意の値とすることもできる。例えば、たとえば、それぞれのマスクに移動平均フィルタなどを用いたフィルタ処理を加えてもよい。
[磁化率算出処理(S1204)]
 上述のようにして3つのマスクを作成した後、磁化率算出処理(ステップS1204)を行う。磁化率算出部540における磁化率算出処理では、マスク算出部530で算出した低信号領域マスク、高磁化率マスク、及び低磁化率マスクと、磁場算出部520で算出した磁場画像とから、第一の実施形態と同様の方法で磁化率画像を算出する。ただしこのとき、第一の実施形態と異なり、脂肪マスクの代わりに高磁化率マスク、水マスクの代わりに低磁化率マスクを用いる。
 すなわち、まず、背景除去処理部841において、磁場画像に対し、低信号領域マスクで定義される領域を背景領域とみなした背景除去処理を行う。次に、背景除去した画像から磁化率算出処理を行い、低磁化率マスクで定義される領域が0となる制約のもと高磁化率マスクで定義される領域の磁化率画像(第一の磁化率画像)を算出する(高磁化率算出部842の処理)。次に、背景除去処理部843において、磁場画像に対し、低信号領域マスクと高磁化率マスクで定義される領域を背景領域とみなした背景除去処理を行う。次に、背景除去処理部843で背景除去した画像から磁化率算出処理を行い、低磁化率マスクで定義される領域の磁化率画像(第二の磁化率画像)を算出する(低磁化率算出部844の処理)。最後に、第一の磁化率画像と第二の磁化率画像を統合することにより統合磁化率画像を算出する(統合磁化率算出部845の処理)。
 以上、説明したように本実施形態のMRI装置は、マスク算出部は、磁場画像を用いて算出した判別用磁化率画像を用いて、高磁化率領域に対応する高磁化率領域マスク、及び、前記高信号領域であって前記高磁化率領域以外の領域に対応する領域マスクを作成する。
 本実施形態により、出血を含む脳画像を算出する場合など、磁化率が大きく異なる二つの領域を含み、かつそれらを水・脂肪画像を用いて分割できない場合に、信号欠損による精度低下を低減して高精度な磁化率画像を算出することができる。
 なお本実施形態では、第一の実施形態と異なり、マスク算出部530で水画像と脂肪画像を用いないため、水脂肪分離処理は必ずしも行わなくてもよい。
<<第三実施形態>>
 第一および第二実施形態では、二つの高信号領域マスクを算出し、それぞれの領域において異なる方法で磁化率画像を算出し、最後に統合することにより統合磁化率画像を算出した。第三実施形態では、三つの高信号領域マスクを算出し、それぞれの領域で異なる方法で磁化率画像を算出し、最後に統合することにより統合磁化率画像を算出する。以下、第一実施形態および第二実施形態と異なる構成に主眼をおいて説明する。
 なお三つの高信号領域マスクでそれぞれ定義される領域は、限定されるものではないが、例えば、出血などの高磁化率となる領域、脂肪領域、それ以外の領域などを想定することができる。以下の実施形態ではこれらを例に説明する。
 本実施形態のMRI装置は、基本的に第一実施形態と同様の構成を有する。
 本実施形態の計算機209は、第一実施形態と同じように、入力装置212を介して設定された(または、あらかじめ設定された)計測パラメータとパルスシーケンスに従いエコーの計測をシーケンサ204に指示し、エコーを計測し、得られたエコーをk空間に配置し、k空間に配置されたエコーに対して演算を行い磁化率画像を算出し、得られた画像を表示装置210に表示する。
 本実施形態の計算機209は、第一実施形態と同じように、計測部300、画像再構成部400、画像変換部500、表示処理部600を備える。図14に、本実施形態における画像変換部500の機能ブロック図を示す。本実施形態の画像変換部500は、第一の実施形態と同様に、水脂肪分離処理部910、磁場算出部920、マスク算出部930、磁化率算出部940を備える。マスク算出部930は、第一実施形態とは異なり、低信号領域マスク算出部931、第一の高信号領域マスク算出部932、第二の高信号領域マスク算出部933、第三の高信号領域マスク算出部934を備える。磁化率算出部940は、第一実施形態とは異なり、背景磁場除去部941、第一の高信号磁化率算出部942、背景磁場除去部943、第二の高信号磁化率算出部944、背景磁場除去部945、第三の高信号磁化率算出部946、統合磁化率算出部947を備える。
 本実施形態の画像変換部500は、図15に示すように、水画像・脂肪画像・周波数画像及び脂肪含有率画像を算出する水脂肪分離処理S1201、周波数画像から磁場画像を算出する磁場算出処理S1202、磁場画像からマスク算出に用いるための暫定的な磁化率画像を算出する処理(マスク算出用磁化率画像算出処理)S1210、複数の領域マスクを算出する処理(マスク算出処理)S1213、磁場画像と複数の領域マスクを用いて磁化率を算出する処理(磁化率算出処理)S1214を行う。図15において、図6(第一実施形態)及び図12(第二実施形態)と同じ符号で示す処理は同様であり、説明を省略し、以下、マスク算出処理S1213及び磁化率算出処理S1214について説明する。
[マスク算出処理(S1213)]
 マスク算出処理S1213では、マスク算出部930が、画像再構成部400で算出した複素画像、水脂肪分離処理部910で算出した脂肪含有率画像、磁場算出部920で算出した磁場画像から、空気などの低信号領域をあらわす低信号領域マスク、相対的に高磁化率な領域をあらわす第一の高信号領域マスク、脂肪領域をあらわす第二の高信号領域マスク、それ以外の高信号領域をあらわす第三の高信号領域マスクを算出する。
 図16は、本実施形態のマスク算出処理のフローチャート、図17はマスク算出に用いる画像と算出されるマスクとの関係を示す図である。マスク算出処理では、まず、低信号領域マスクMと高信号領域マスクMを算出する(ステップS1701)。次に、高信号領域マスクMおよびマスク用磁化率画像から第一の高信号領域マスクMを算出する(ステップS1702)。次に、高信号領域マスクM、第一の高信号領域マスクM、脂肪含有率画像から、第二の高信号領域マスクMを算出する(ステップS1703)。最後に、高信号領域マスクM、第一の高信号領域マスクM、第二の高信号領域マスクMから、第三の高信号領域マスクMを算出する(ステップS1704)。以下、各処理の詳細を説明する。
 ステップS1701の低信号領域マスク算出方法は、基本的に第一の実施形態と同様であり、複素画像からマスク用絶対値画像を算出し、マスク用絶対値画像を閾値処理し、低信号領域マスクMおよび高信号領域マスクMを算出する(低信号領域マスク算出部931の処理)。
 次に、ステップS1702において、第二の実施形態における高磁化率マスクの算出方法と同様の方法で、ステップS1210で算出したマスク用磁化率画像を判別用画像とし、任意の閾値を用いて閾値処理することによりマスク(第一の高信号領域マスク)Mを算出する(第一の高信号マスク算出部932の処理)。
 次に、ステップS1703において、第二の高信号領域マスクMを算出する。このため、まず、水脂肪分離処理S1201で算出した脂肪含有率画像を任意の閾値を用いて閾値処理することによりマスクMdを算出する。例えば、脂肪含有率が当該閾値より大きい領域を1、脂肪含有率が当該閾値より小さい領域を0とする。本実施形態では、閾値は0.8に設定する。そして、閾値処理により算出したマスクMdと、高信号領域マスクMから第一の高信号領域マスクMを減算したマスクM(H-1)をかけ合わせることにより、第二の高信号領域マスクMを算出する。
 最後に、高信号領域マスクMから第一の高信号領域マスクMおよび第二の高信号領域マスクMを減算することにより、第三の高信号領域マスクMを算出する(ステップS1704)。
 なお本実施形態では、算出した四つのマスクは0と1の値のみをとるが、第一の実施形態および第二の実施形態と同様、任意の値とすることもできる。また、第一の実施形態および第二の実施形態と同様、算出したそれぞれのマスクに対して、孤立点除去処理などの追加の処理を加えてもよい。
[磁化率算出処理(S1214)]
 次に、本実施形態における磁化率算出処理S1214の流れについて説明する。磁化率算出部940は、マスク算出部930で算出した低信号領域マスクML、第一の高信号領域マスクM1、第二の高信号領域マスクM2、第三の高信号領域マスクM3と磁場算出部920で算出した磁場画像から、磁化率画像を算出する。
 図18は、本実施形態の磁化率算出処理のフローチャートである。本実施形態では、磁場画像について、まず低信号領域マスクで定義される領域を背景領域とみなした背景除去処理を行う(ステップS1801)。次に、ステップS1801で背景除去した磁場画像から磁化率算出処理を行い、第二の高信号領域マスクM2と第三の高信号領域マスクM3の和で定義される領域の磁化率を0とする制約のもと、第一の高信号磁化率画像を算出する(ステップS1802)。次に、磁場画像について、低信号領域マスクMLと第一の高信号領域マスクM1で定義される領域を背景とみなした背景除去処理を行う(ステップS1803)。次に、ステップS1803で背景除去した磁場画像から、第三の高信号領域マスクM3を0とする制約のもと、第二の高信号磁化率画像を算出する(ステップS1804)。次に、低信号領域マスクMLと第一の高信号領域マスクM1と第二の高信号領域マスクM2で定義される領域を背景とみなした背景除去処理を行う(ステップS1805)。次に、ステップS1805で背景除去した画像から磁化率画像を算出する(ステップS1806)。
 なお、ステップS1801、ステップS1803、ステップS1805における背景磁場除去法は、背景領域の違いを除き、基本的には第一の実施形態に記した方法と同様である。またステップS1802、ステップS1804における、一部の領域の磁化率を0とする制約のもと磁化率画像を算出する方法は、関心領域および制約をおく領域の違いを除き、基本的には第一の実施形態に記した方法と同様である。ステップS1806における磁化率画像の算出方法は、関心領域の違いを除き、基本的には第一の実施形態におけるステップS1504の方法と同様である。
 最後に、第一の高信号磁化率画像、第二の高信号磁化率画像、第三の高信号磁化率画像を、式(22)により統合し、統合磁化率画像χを算出する(ステップS1807)。
Figure JPOXMLDOC01-appb-M000020
ここで、M、M、Mはそれぞれ第一の高信号領域マスク、第二の高信号領域マスク、第三の高信号領域マスク、χ、χ、χはそれぞれ第一の高信号磁化率画像、第二の高信号磁化率画像、第三の高信号磁化率画像をあらわす。
 本実施形態のMRI装置は、画像変換部が、前記画像再構成部が作成した画像を用いて、前記被検体に含まれる成分の内、水プロトンからの信号が支配的である水画像、脂肪プロトンからの信号が支配的である脂肪画像を分離する水脂肪分離部をさらに備え、前記マスク算出部は、前記水画像及び前記脂肪画像を用いて画素毎に算出した脂肪含有率と前記磁場画像を用いて算出した判別用磁化率画像とを用いて、高磁化率領域に対応する高磁化率領域マスク、脂肪が支配的な領域に対応する脂肪マスク、及び、前記高磁化率領域マスク及び脂肪マスクで定義される領域以外の領域に対応する領域マスクを作成する。
 本実施形態により、出血領域、脂肪領域、それ以外の水領域など、磁化率が大きく異なる三つの領域を含む画像の場合に、信号欠損による精度低下を低減して高精度な磁化率画像を算出することができる。
 なお、三つの場合に限らず四つ以上の領域に対しても本処理を同様に適用することが可能である。また、二つの閾値を用いてマスク用磁化率画像を閾値処理することにより、三つの高信号領域マスクを算出してもよい。
<第一実施形態の実施例>
 下表1に示すように、脂肪内磁化率を0.61、脂肪外磁化率を0とするモデルを想定し、このモデルについて、第一実施形態の磁化率算出法(発明方法)及び従来法により磁化率画像をシミュレーションにより得た。従来法は、低信号領域の背景除去後に高信号領域の磁化率を平滑化する制約のもと磁化率を算出する手法を採用した。
 シミュレーションの結果を表1及び図19に示す。図19の(a)、(b)、(c)は、それぞれ、モデル、従来法、発明方法の磁化率画像、(d)、(e)はモデルの磁化率画像との差で、(d)は従来法、(e)は発明方法を示す。
Figure JPOXMLDOC01-appb-T000021
 シミュレーション結果から、発明方法では誤差が大幅に低減し磁化率算出の精度が向上することが確認された。
 100:垂直磁場MRI装置、101:水平磁場MRI装置、102:MRI装置、201:マグネット、202:傾斜磁場コイル、203:被検体、204:シーケンサ、205:傾斜磁場電源、206:高周波磁場発生器、207:プローブ、208:受信器、209:計算機、210:表示装置、211:記憶装置、212:入力装置、300:計測部、400:画像再構成部、500:画像変換部、510:水脂肪分離処理部、520:磁場算出部、530:マスク算出部、531:低信号領域マスク算出部、532:脂肪マスク算出部、533:水マスク算出部、540:磁化率算出部、541:背景磁場除去部、542:脂肪磁化率算出部、543:背景磁場除去部、544:水磁化率算出部、545:統合磁化率算出部、550:マスク用磁化率算出部、600:表示処理部、710:GrE型パルスシーケンス、711:RFパルス、712:スライス選択傾斜磁場、713:位相エンコード傾斜磁場、721;読み出し傾斜磁場、722;読み出し傾斜磁場、723;読み出し傾斜磁場、724;読み出し傾斜磁場、731:第一エコー信号、732:第二エコー信号、733:第三エコー信号、734:第四エコー信号、810:低信号領域マスク算出部、820:高磁化率マスク算出部、830:低磁化率マスク算出部、910:水脂肪分離処理部、920:磁場算出部、930:マスク算出部、940:磁化率算出部

Claims (15)

  1.  静磁場磁石と、前記静磁場磁石が形成する空間に傾斜磁場及び高周波磁場を印加する磁場発生部と、前記空間に置かれた被検体から発生する核磁気共鳴信号を受信する受信部と、前記核磁気共鳴信号に対し演算を行う計算機と、を備えた磁気共鳴イメージング装置であって、
     前記計算機は、前記磁場発生部及び前記受信部の動作を制御し、磁場不均一の影響を含む核磁気共鳴信号を計測する計測部と、前記核磁気共鳴信号を用いて画像再構成を行い、複素画像を作成する画像再構成部と、前記複素画像を用いて前記被検体の磁化率画像を算出する画像変換部と、を備え、
     前記画像変換部は、前記複素画像から磁場反映画像を算出する磁場算出部と、前記画像再構成部が作成した画像を用いて、低信号領域に対応する低信号領域マスク及び高信号領域内の複数の領域に対応する複数の領域マスクを算出するマスク算出部と、前記磁場算出部が算出した磁場反映画像と前記マスク算出部が算出した前記低信号領域マスク及び複数の領域マスクとを用いて、前記複数の領域のそれぞれについて磁化率を算出する磁化率算出部と、を備え、
     前記磁化率算出部は、前記低信号領域を背景とする制約、及び、前記複数の領域のうち特定の領域の磁化率を特定値とする制約のもと、前記特定の領域とは異なる領域の磁化率を算出することを特徴とする磁気共鳴イメージング装置。
  2.  請求項1に記載の磁気共鳴イメージング装置において、
     前記磁化率算出部は、前記低信号領域及び前記特定の領域とは異なる領域を背景とする制約のもと、前記特定の領域の磁化率を算出し、前記特定の領域とは異なる領域の磁化率と統合することを特徴とする磁気共鳴イメージング装置。
  3.  請求項1又は2に記載の磁気共鳴イメージング装置であって、
     前記画像変換部は、前記画像再構成部が作成した画像を用いて、水プロトンからの信号が支配的である水画像、脂肪プロトンからの信号が支配的である脂肪画像を分離する水脂肪分離部をさらに備え、
     前記マスク算出部は、前記水画像及び前記脂肪画像を用いて算出した画素毎の脂肪含有率の所定の値を閾値として用い、前記複数の領域マスクとして、水が支配的な領域に対応する水マスク及び脂肪が支配的な領域に対応する脂肪マスクを作成することを特徴とする磁気共鳴イメージング装置。
  4.  請求項3に記載の磁気共鳴イメージング装置であって、
     前記磁化率算出部による算出において、前記特定の領域は、水が支配的な領域であり、前記特定値はゼロであることを特徴とする磁気共鳴イメージング装置。
  5.  請求項3に記載の磁気共鳴イメージング装置であって、
     前記磁化率算出部による算出において、前記特定の領域は、脂肪が支配的な領域であり、前記特定値は予め定めた正の定数であることを特徴とする磁気共鳴イメージング装置。
  6.  請求項3に記載の磁気共鳴イメージング装置であって、
     前記計測部は、エコー時間が異なる3以上の核磁気共鳴信号を計測し、
     前記水脂肪分離部は、前記エコー時間が異なる複数の核磁気共鳴信号の信号モデルへのフィッティングによって、前記水画像、前記脂肪画像及び前記被検体の周波数分布を表す周波数画像を算出することを特徴とする磁気共鳴イメージング装置。
  7.  請求項6に記載の磁気共鳴イメージング装置であって、
     前記磁場算出部は、前記水脂肪分離部が算出した周波数画像から前記磁場反映画像を算出することを特徴とする磁気共鳴イメージング装置。
  8.  請求項1に記載の磁気共鳴イメージング装置であって、
     前記マスク算出部は、前記磁場反映画像を用いて算出した判別用磁化率画像を用いて、高磁化率領域に対応する高磁化率領域マスク、及び、前記高信号領域であって前記高磁化率領域以外の領域に対応する領域マスクを作成することを特徴とする磁気共鳴イメージング装置。
  9.  請求項2に記載の磁気共鳴イメージング装置であって、
     前記画像変換部は、前記画像再構成部が作成した画像を用いて、前記被検体に含まれる成分の内、水プロトンからの信号が支配的である水画像、脂肪プロトンからの信号が支配的である脂肪画像を分離する水脂肪分離部をさらに備え、
     前記マスク算出部は、前記水画像及び前記脂肪画像を用いて画素毎に算出した脂肪含有率と前記磁場反映画像を用いて算出した判別用磁化率画像とを用いて、高磁化率領域に対応する高磁化率領域マスク、脂肪が支配的な領域に対応する脂肪マスク、及び、前記高磁化率領域マスク及び脂肪マスクで定義される領域以外の領域に対応する領域マスクを作成することを特徴とする磁気共鳴イメージング装置。
  10.  磁場不均一の影響を含む核磁気共鳴信号を計測し、
     前記核磁気共鳴信号を用いて複素画像を作成し、
     前記複素画像を用いて磁場反映画像を作成し、前記磁場反映画像を用いて磁化率画像を作成する磁気共鳴イメージング方法であって、
     前記複素画像から作成した絶対値画像を用いて、低信号領域マスク、高信号マスクを作成し、前記高信号マスクから、任意の判別用画像の画素値が所定の閾値より大きい第一の領域を表す第一の領域マスク及び前記判別用画像の画素値が所定の閾値より小さい第二の領域を表す第二の領域マスクを作成し、
     前記磁場反映画像と前記低信号領域マスク、第一及び第二の領域マスクを用いて、前記信号マスクで定義される領域を背景とする制約及び前記第二の領域の磁化率を特定値とする制約のもと前記第一の領域の磁化率を算出し、前記信号マスクで定義される領域及び前記第一の領域を背景とする制約のもと前記第二の領域の磁化率画像を算出することを特徴とする磁気共鳴イメージング方法。
  11.  請求項10に記載の磁気共鳴イメージング方法であって、
     前記第一の領域の磁化率の算出は、さらに、磁化率を平滑化する制約を含むことを特徴とする磁気共鳴イメージング方法。
  12.  請求項11に記載の磁気共鳴イメージング方法であって、
     前記磁化率の算出において、前記特定値とする制約の強さを表すパラメータの値と、前記平滑化する制約の強さを表すパラメータの値とを、前記第一の領域マスクで定義される領域における磁化率の平均値と標準偏差が最大になるときの値に設定することを特徴とする磁気共鳴イメージング方法。
  13.  請求項10に記載の磁気共鳴イメージング方法であって、
     前記複数の領域マスクは、水プロトンからの信号が支配的である領域に対応する水マスクと、脂肪プロトンからの信号が支配的である領域に対応する脂肪マスクとを含み、
     磁化率の算出において、前記水マスクで定義される領域の特定値をゼロとして磁化率を算出することを特徴とする磁気共鳴イメージング方法。
  14.  請求項13に記載の磁気共鳴イメージング方法であって、さらに
     前記複素画像を用いて、水が支配的である領域の水画像、及び、脂肪が支配的である領域の脂肪画像を作成し、
     前記水画像と前記脂肪画像を用いて、画素毎の脂肪含有率を算出し、
     前記画素毎の脂肪含有率の所定の値を閾値として、前記高信号マスクを前記第一及び第二の領域マスクに分割することを特徴とする磁気共鳴イメージング方法。
  15.  コンピュータに、
     磁気共鳴イメージング装置で計測された、磁場不均一の影響を含む核磁気共鳴信号を用いて複素画像を作成するステップと、
     前記複素画像を用いて磁場反映画像を作成するステップと、
     前記複素画像から絶対値画像を作成し、当該絶対値画像を用いて、低信号領域マスク、高信号マスクを作成し、前記高信号マスクから複数の領域に対応する複数の領域マスクを作成するステップと、
     前記磁場反映画像と前記低信号領域マスク及び前記複数の領域マスクとを用いて、前記低信号領域マスクで定義される領域を背景とする制約、及び、前記複数の領域のうち特定の領域の磁化率を特定値とする制約のもとで、前記特定の領域以外の領域の第一の磁化率を算出するステップと、
     前記低信号領域マスクで定義される領域及び前記特定の領域以外に対応する領域マスクで定義される領域を背景とする制約のもとで、前記特定の領域の第二の磁化率を算出するステップと、
     前記第一の磁化率と前記第二の磁化率を統合するステップと、
     を実行させるための磁化率算出プログラム。
PCT/JP2017/019381 2016-05-25 2017-05-24 磁気共鳴イメージング装置、磁気共鳴イメージング方法及び磁化率算出プログラム WO2017204251A1 (ja)

Priority Applications (1)

Application Number Priority Date Filing Date Title
US16/097,284 US11141078B2 (en) 2016-05-25 2017-05-24 Magnetic resonance imaging device, magnetic resonance imaging method and susceptibility calculation program

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
JP2016-104358 2016-05-25
JP2016104358A JP6608764B2 (ja) 2016-05-25 2016-05-25 磁気共鳴イメージング装置、磁気共鳴イメージング方法及び磁化率算出プログラム

Publications (1)

Publication Number Publication Date
WO2017204251A1 true WO2017204251A1 (ja) 2017-11-30

Family

ID=60412402

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2017/019381 WO2017204251A1 (ja) 2016-05-25 2017-05-24 磁気共鳴イメージング装置、磁気共鳴イメージング方法及び磁化率算出プログラム

Country Status (3)

Country Link
US (1) US11141078B2 (ja)
JP (1) JP6608764B2 (ja)
WO (1) WO2017204251A1 (ja)

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP6776217B2 (ja) * 2017-11-24 2020-10-28 株式会社日立製作所 画像処理装置、画像処理方法、画像処理プログラム及び磁気共鳴イメージング装置
CN108829639B (zh) * 2018-07-06 2023-10-27 上海联影医疗科技股份有限公司 一种磁共振成像方法和设备
US10955511B2 (en) * 2019-01-25 2021-03-23 Siemens Healthcare Gmbh Magnetic resonance imaging coil normalization by using a reference image
US11475558B2 (en) * 2019-11-13 2022-10-18 Raytheon Company Organ isolation in scan data

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20140142417A1 (en) * 2012-11-21 2014-05-22 Wisconsin Alumni Research Foundation System and Method for Assessing Susceptibility of Tissue Using Magnetic Resonance Imaging
US20150002148A1 (en) * 2013-06-26 2015-01-01 Medimagemetric LLC Joint estimation of chemical shift and quantitative susceptibility map using mri signal

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP4251763B2 (ja) * 2000-08-11 2009-04-08 株式会社日立メディコ 磁気共鳴イメージング装置
WO2014057716A1 (ja) * 2012-10-10 2014-04-17 株式会社日立製作所 磁気共鳴イメージング装置

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20140142417A1 (en) * 2012-11-21 2014-05-22 Wisconsin Alumni Research Foundation System and Method for Assessing Susceptibility of Tissue Using Magnetic Resonance Imaging
US20150002148A1 (en) * 2013-06-26 2015-01-01 Medimagemetric LLC Joint estimation of chemical shift and quantitative susceptibility map using mri signal

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
KHALIDOV,I. ET AL.: "Susceptibility mapping: computation of the field map using water-fat separation at 7T", PROCEEDINGS OF INTERNATIONAL SOCIETY FOR MAGNETIC RESONANCE IN MEDICINE, vol. 19, 2011, pages 4475 *
WEI,HONGJIANG ET AL.: "IMAGING MAGNETIC SUSCEPTIBILITY OF THE HUMAN KNEE JOINT AT 3 AND 7 TESLA", PROCEEDINGS OF INTERNATIONAL SOCIETY FOR MAGNETIC RESONANCE IN MEDICINE, vol. 23, 5 June 2015 (2015-06-05), pages 288 *

Also Published As

Publication number Publication date
US20190076049A1 (en) 2019-03-14
JP6608764B2 (ja) 2019-11-20
US11141078B2 (en) 2021-10-12
JP2017209288A (ja) 2017-11-30

Similar Documents

Publication Publication Date Title
JP5843876B2 (ja) 磁気共鳴イメージング装置および磁化率強調画像生成方法
CN107072592B (zh) 磁共振成像装置以及定量性磁化率匹配方法
JP5982494B2 (ja) 磁気共鳴イメージング装置
JP6085545B2 (ja) 磁気共鳴イメージング装置、画像処理装置および磁化率画像算出方法
WO2014076808A1 (ja) 磁気共鳴イメージング装置および定量的磁化率マッピング法
JP6595393B2 (ja) 磁気共鳴イメージング装置、及び、画像処理方法
WO2017204251A1 (ja) 磁気共鳴イメージング装置、磁気共鳴イメージング方法及び磁化率算出プログラム
US20150355298A1 (en) Method and device for accurate quantification of t2 relaxation times based on fast spin-echo nmr sequences
US11241163B2 (en) Measuring blood vessel characteristics with MRI
US10617343B2 (en) Methods and systems for quantitative brain assessment
JP7230149B2 (ja) 磁気共鳴イメージングで得た画像の処理方法、画像処理プロブラム、及び、計算機
JP5636058B2 (ja) 磁気共鳴撮影装置
US10114100B2 (en) Method for magnetic resonance imaging
Wang et al. High-fidelity direct contrast synthesis from magnetic resonance fingerprinting
JP6113012B2 (ja) 磁気共鳴イメージング装置及び補正用b1マップを計算する方法
JP6776217B2 (ja) 画像処理装置、画像処理方法、画像処理プログラム及び磁気共鳴イメージング装置
Xu et al. Improving the lesion appearance on FLAIR images synthetized from quantitative MRI: a fast, hybrid approach
Milles et al. Computation of Transmitted and Received $ rm B1 $ Fields in Magnetic Resonance Imaging
Banerjee High resolution high field quantitative parallel magnetic resonance imaging for osteoporosis and other clinical applications

Legal Events

Date Code Title Description
NENP Non-entry into the national phase

Ref country code: DE

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

Ref document number: 17802837

Country of ref document: EP

Kind code of ref document: A1

122 Ep: pct application non-entry in european phase

Ref document number: 17802837

Country of ref document: EP

Kind code of ref document: A1