WO2010137516A1 - 磁気共鳴撮影装置 - Google Patents

磁気共鳴撮影装置 Download PDF

Info

Publication number
WO2010137516A1
WO2010137516A1 PCT/JP2010/058564 JP2010058564W WO2010137516A1 WO 2010137516 A1 WO2010137516 A1 WO 2010137516A1 JP 2010058564 W JP2010058564 W JP 2010058564W WO 2010137516 A1 WO2010137516 A1 WO 2010137516A1
Authority
WO
WIPO (PCT)
Prior art keywords
phase
correction
magnetic resonance
resonance imaging
imaging apparatus
Prior art date
Application number
PCT/JP2010/058564
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 JP2011515997A priority Critical patent/JP5165791B2/ja
Publication of WO2010137516A1 publication Critical patent/WO2010137516A1/ja

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01R—MEASURING ELECTRIC VARIABLES; MEASURING MAGNETIC VARIABLES
    • G01R33/00—Arrangements or instruments for measuring magnetic variables
    • G01R33/20—Arrangements or instruments for measuring magnetic variables involving magnetic resonance
    • G01R33/44—Arrangements or instruments for measuring magnetic variables involving magnetic resonance using nuclear magnetic resonance [NMR]
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01R—MEASURING ELECTRIC VARIABLES; MEASURING MAGNETIC VARIABLES
    • G01R33/00—Arrangements or instruments for measuring magnetic variables
    • G01R33/20—Arrangements or instruments for measuring magnetic variables involving magnetic resonance
    • G01R33/44—Arrangements or instruments for measuring magnetic variables involving magnetic resonance using nuclear magnetic resonance [NMR]
    • G01R33/48—NMR imaging systems
    • G01R33/483—NMR imaging systems with selection of signals or spectra from particular regions of the volume, e.g. in vivo spectroscopy
    • G01R33/4833—NMR imaging systems with selection of signals or spectra from particular regions of the volume, e.g. in vivo spectroscopy using spatially selective excitation of the volume of interest, e.g. selecting non-orthogonal or inclined slices

Definitions

  • the present invention relates to spectroscopic imaging, and more particularly to a technique for improving the accuracy of eddy current correction for correcting spectral distortion caused by eddy currents.
  • MRS Magnetic Resonance Spectroscopy
  • CSI Chemical Shift
  • Major metabolites of the human body that can be detected by MRS and MRSI include choline (Cho), creatine (Cr), N-acetylaspartic acid (NAA), lactic acid (Lac), and the like. From the amount of these metabolites, it is possible to determine the degree of progression and early diagnosis of metabolic disorders such as cancer. It is also considered possible to perform non-invasive diagnosis of tumor malignancy.
  • an eddy current is generated by a gradient magnetic field applied during measurement. Eddy currents cause static magnetic field inhomogeneities spatially and temporally and distort the shape of the spectrum obtained by measurement.
  • This spectral distortion is generally corrected using the phase value of the reference signal.
  • a water signal having a signal intensity greater than that of the metabolite is used as the reference signal.
  • a phase value is acquired from this water signal and used for correction (for example, see Non-Patent Document 1).
  • spatial and temporal phase values are calculated from this water signal, phase correction is performed on metabolite image data, and spectral distortion due to eddy current is corrected.
  • CSI or MRSI has a very small number of matrices (voxels) to be measured, about 8 ⁇ 8 to 32 ⁇ 32, from the viewpoint of measurement time and SNR (signal to noise ratio). For this reason, truncation occurs due to the Fourier transform performed in the image reconstruction, and a signal of a distant voxel is mixed. As a result, when the static magnetic field inhomogeneity exists, a water signal having a frequency different from that of the water signal in the target voxel is mixed.
  • phase jump region occurs in which the variation in the phase value per unit time is prominently larger than other locations in the time variation of the phase value used for eddy current correction. To do.
  • the magnitude of the phase change amount in this phase jump region is proportional to the concentration of the mixed water signal.
  • Fig. 13 (a) shows the spectrum of the water signal by computer simulation.
  • Water 1 and Water 2 in FIG. 13A are water signal spectra having frequencies of 2 Hz and 5 Hz, respectively, and a concentration ratio of 1.0 to 0.9.
  • FIG. 13B shows the time change of the phase value in the time domain when the two signals are mixed. As shown in FIGS. 13A and 13B, a phase jump proportional to the concentration ratio occurs at a time interval of 1 / ⁇ f with respect to the frequency difference ⁇ f.
  • static magnetic field fluctuations due to eddy currents are added to cause a gradual phase change.
  • FIG. 13C shows a state before eddy current correction
  • FIG. 13D shows a state after eddy current correction.
  • FIG. 13D when correction is performed using phase data having a phase jump, ringing artifacts are generated by the phase jump, and the spectrum is deteriorated conversely by the eddy current correction process.
  • Non-Patent Document 2 there is a method of reducing ringing artifacts by applying a low-pass filter to the spectrum of the water signal (for example, see Non-Patent Document 2).
  • a method of correcting a phase jump appearing in the phase data and reducing ringing artifacts see, for example, Non-Patent Document 3.
  • the timing at which the absolute value intensity of the water signal in the time domain takes an extreme value is defined as the occurrence of a phase jump.
  • a first derivative with respect to the time t is calculated with respect to the phase data of the water signal in the time domain, and the periphery of the phase jump occurrence point is fitted with a model function.
  • the fitted region is set as a phase skip correction region, and the phase skip is corrected by linearly interpolating the phase value.
  • Non-Patent Document 2 cannot completely eliminate the ringing artifact because the phase skipping low frequency component remains. Further, since the high-frequency component of the eddy current is cut, a sufficient eddy current correction effect cannot be obtained.
  • Non-Patent Document 3 when a phase jump occurs even if there is no extreme value in the absolute value intensity of the water signal in the time domain, the phase jump is not corrected. In addition, when there are a plurality of locations where the phase change is steep, it is difficult to identify and extract the location where the phase jumps. Furthermore, since the fitting area is not specified, the fitting accuracy decreases when there are a plurality of phase jumps having different phase change amounts. That is, in order to increase the fitting accuracy, it is necessary to perform fitting for each phase jump occurrence location, and the process becomes complicated and it is difficult to automatically correct the phase jump.
  • the present invention has been made in view of the above circumstances, and an object of the present invention is to prevent artifacts due to correction processing and improve the accuracy of correction results when performing processing for correcting the influence of eddy currents in an MRI apparatus.
  • the phase data includes the phase change amount in the time domain.
  • the region where the phase skip occurs is specified, and the phase skip correction is performed only for the specified location.
  • the magnitude of the amount of change is determined in each area divided according to a predetermined area division.
  • a magnetic resonance imaging apparatus comprising: receiving means for receiving a nuclear magnetic resonance signal from the subject; and eddy current correction means for performing eddy current correction on the nuclear magnetic resonance signal, wherein the eddy current correction means comprises: Phase data calculation means for calculating phase data of a nuclear magnetic resonance signal of a substance having a signal intensity greater than that of a metabolite to be measured; phase data dividing means for dividing the phase data into a plurality of small regions in the time direction; Phase jump generation area extracting means for extracting a phase jump occurrence area in which phase jump exists from the small area, and phase jump correction means for correcting the phase jump in the phase jump generation area are provided.
  • the eddy current correction means comprises: Phase data calculation means for calculating phase data of a nuclear magnetic resonance signal of a substance having a signal intensity greater than that of a metabolite to be measured; phase data dividing means for dividing the phase data into a plurality of small regions in the time direction; Phase jump generation area extracting means for extracting a phase jump occurrence
  • (A)-(c) is an external view of the magnetic resonance imaging apparatus of embodiment of this invention.
  • (A) is a functional block diagram of the magnetic resonance imaging apparatus of embodiment of this invention
  • (b) is a functional block diagram of the computer with which the magnetic resonance imaging apparatus of embodiment of this invention is provided. It is a figure which shows an example of the MRSI pulse sequence of embodiment of this invention.
  • (A)-(c) is a figure for demonstrating the area
  • (A) is a flowchart of the area dividing process of the embodiment of the present invention
  • (b) and (c) are diagrams for explaining the area dividing process of the embodiment of the present invention. It is a figure for demonstrating the phase variation calculation process of embodiment of this invention. It is a figure for demonstrating the phase jump correction area
  • (A) is the phase data before performing the phase skip correction processing of the embodiment of the present invention and the phase data after correcting one phase skip region.
  • (B) is a figure which shows the phase data before implementing the phase jump correction process of embodiment of this invention, and the phase data after correct
  • (A)-(c) is a figure for demonstrating the effect by embodiment of this invention using a computer simulation result.
  • (A)-(c) is a figure for demonstrating the effect by embodiment of this invention using a phantom experiment result.
  • (A)-(d) is a figure for demonstrating the cause which a ringing artifact arises when eddy current correction is performed using the phase data of a water signal.
  • FIG. 1 is an external view of the MRI apparatus of this embodiment.
  • FIG. 1A shows a horizontal magnetic field type MRI apparatus 100 using a tunnel magnet that generates a static magnetic field with a solenoid coil.
  • FIG. 1B shows a hamburger type (open type) vertical magnetic field type MRI apparatus 110 in which magnets are separated into upper and lower sides in order to enhance the feeling of opening.
  • FIG. 1C shows an MRI apparatus 120 that uses the same tunnel-type magnet as in FIG. 1A and has a feeling of openness by shortening the depth of the magnet and tilting it obliquely. In the present embodiment, any of these MRI apparatuses having these appearances can be used.
  • the MRI apparatus of the present embodiment is not limited to these forms.
  • various known MRI apparatuses can be used regardless of the form and type of the apparatus.
  • the MRI apparatus 100 is representative.
  • FIG. 2A is a functional configuration diagram of the MRI apparatus 100 of the present embodiment.
  • the MRI apparatus 100 of the present embodiment generates a static magnetic field coil 2 that generates a static magnetic field and a gradient magnetic field in each of the x, y, and z directions in the space where the subject 1 is placed.
  • the shim coil 4 that adjusts the static magnetic field distribution
  • the measurement high-frequency coil 5 that irradiates the measurement region of the subject 1 with a high-frequency magnetic field
  • a receiving high-frequency coil 6 that receives a nuclear magnetic resonance signal (hereinafter simply referred to as a receiving coil).
  • the gradient magnetic field coil 3 and the shim coil 4 are driven by a gradient magnetic field power supply unit 12 and a shim power supply unit 13, respectively.
  • a case where separate transmission coils 5 and reception coils 6 are used will be described as an example.
  • the transmission coil 5 and the reception coil 6 are configured as a single coil. May be.
  • the high frequency magnetic field irradiated by the transmission coil 5 is generated by the transmitter 7.
  • the nuclear magnetic resonance signal detected by the receiving coil 6 is sent to the computer 9 through the receiver 8.
  • the sequence control device 14 controls the operations of the gradient magnetic field power supply unit 12 that is a drive power supply for the gradient magnetic field generating coil 3, the shim power supply unit 13 that is the drive power supply for the shim coil 4, the transmitter 7, and the receiver 8. Control the timing of application of gradient magnetic field and high-frequency magnetic field and reception of nuclear magnetic resonance signals.
  • the control time chart is called a pulse sequence, is preset according to measurement, and is stored in a storage device or the like provided in the computer 9 described later.
  • the computer 9 performs various arithmetic processes on the received nuclear magnetic resonance signal, generates image information and spectrum information, and controls the overall operation of the MRI apparatus 100.
  • the computer 9 is an information processing device including a CPU, a memory, a storage device, and the like, and a display 10, an external storage device 11, an input device 15, and the like are connected to the computer 9.
  • the display 10 is an interface that displays results obtained by the arithmetic processing to the operator.
  • the input device 15 is an interface for an operator to input conditions, parameters, and the like necessary for arithmetic processing performed in the present embodiment.
  • the external storage device 11 holds data used for various types of arithmetic processing executed by the computer 9, data obtained by the arithmetic processing, input conditions, parameters, and the like together with the storage device included in the computer 9.
  • FIG. 2B is a functional block diagram of the computer 9 of this embodiment.
  • the computer 9 of this embodiment includes a measurement control unit 910, a display information generation unit 920, and an eddy current correction unit 930.
  • the measurement control unit 910 operates the sequence control device 14 according to the pulse sequence and controls each unit to perform measurement.
  • the display information generation unit 920 performs various arithmetic processes on the nuclear magnetic resonance signal obtained by measurement to generate image information and spectrum information.
  • the eddy current correction unit 930 performs eddy current correction processing for correcting spectral distortion due to eddy current.
  • spectral distortion due to eddy current is corrected by phase data of a water signal having a signal intensity greater than that of a metabolite.
  • the eddy current correction unit 930 corrects all phase jumps in the phase data by repeating calculation with a phase data calculation unit 940 that calculates the phase data of the water signal from the obtained nuclear magnetic resonance signal.
  • a phase data correction unit 950 Furthermore, the phase data correction unit 950 of the present embodiment extracts an area where phase jump occurs from the divided area 953 that divides the phase data into a plurality of small areas in the time direction.
  • a phase jump generation area extraction unit 951 and a phase jump correction unit 952 that corrects the phase jump in the small area extracted by the phase jump generation area extraction unit 951.
  • the various functions realized by the computer 9 are realized by the CPU loading a program stored in the storage device into the memory and executing it.
  • at least one function is realized by an information processing apparatus that is independent of the MRI apparatus 100 and capable of transmitting and receiving data to and from the MRI apparatus 100. It may be.
  • a water signal is used as the eddy current correction signal. Therefore, first, the measurement control unit 910 controls the sequence control device 14 according to a predetermined pulse sequence, performs non-water suppression measurement, and obtains a water signal F (t).
  • t is a sampling time and is a variable taking a discrete value.
  • the phase data calculation unit 940 calculates phase data ⁇ (t).
  • the phase data correction unit 950 obtains phase data ⁇ c (t) in which all phase skips existing in the phase data ⁇ (t) are corrected.
  • the measurement control unit 910 controls the sequence control device 14 according to a predetermined pulse sequence, performs water suppression measurement, and obtains a metabolite signal S (t).
  • the eddy current correction unit 930 corrects the obtained metabolite signal S (t) with the phase data ⁇ c (t) to obtain the metabolite signal S ecc (t) after the eddy current correction.
  • the display information generation unit 920 performs a Fourier transform on the metabolite signal S ecc (t) to obtain a spectrum or distribution image of the metabolite.
  • an example of a pulse sequence used by the measurement control unit 910 for the measurement will be described.
  • a description will be given by taking as an example a region selective MRSI pulse sequence (hereinafter referred to as an MRSI pulse sequence) for imaging a metabolite.
  • FIG. 3 is an example of an MRSI pulse sequence 300.
  • RF indicates the application timing of the high frequency magnetic field pulse.
  • Gx, Gy, and Gz indicate application timings of gradient magnetic field pulses in the x, y, and z directions, respectively.
  • a / D indicates a signal measurement period Tp1.
  • the MRSI pulse sequence 300 shown in FIG. 3 is the same as a known MRSI pulse sequence. That is, a predetermined time region of interest is selectively excited using one excitation pulse RF1 and two inversion pulses RF2 and RF3, and an FID signal (free induction decay) FID1 is obtained from this region of interest. Repeat at TR interval.
  • FIG. 4A and 4B are scout images for positioning obtained by measurement performed prior to the main measurement.
  • FIG. 4A is a trans image
  • FIG. 4B is a sagittal image
  • FIG. 4C is a coronal image. It is.
  • the relationship between the operation of each part and the excited region will be described with reference to FIGS.
  • a high-frequency magnetic field RF1 and gradient magnetic field pulses Gs1 and Gs1 'in the z direction are applied to excite the cross section 401 in the z direction.
  • a high frequency magnetic field RF2 and a gradient magnetic field pulse Gs2 in the y direction are applied after TE / 4 (where TE is an echo time).
  • TE is an echo time
  • the high frequency magnetic field RF3 and the gradient magnetic field pulse Gs3 in the x direction are applied after TE / 2 from the application of the high frequency magnetic field RF2.
  • the phase of nuclear magnetization in the region of interest 404 where the cross section 401 in the z direction, the cross section 402 in the y direction, and the cross section 403 in the x direction intersect is inverted, and a free induction decay signal FID1 is generated therefrom.
  • This free induction decay signal FID1 is measured.
  • the gradient magnetic field pulses Gd1 to Gd3 and Gd1 ′ to Gd3 ′ in each direction do not disturb the phase of the nuclear magnetization excited by the high-frequency magnetic field RF1, but dephase the phase of the nuclear magnetization excited by the RF2 and RF3. It is a gradient magnetic field.
  • the phase encode gradient magnetic fields Gp1 and Gp2 are applied after the high-frequency magnetic field RF3.
  • the nuclear magnetic resonance signal of the region of interest 404 is obtained.
  • the eddy current correction unit 930 uses the phase signal ⁇ c (t) of the water signal after the phase jump correction from the metabolite signal S (t), and the metabolite signal S ecc (t) after the eddy current correction, It calculates according to the following formula
  • S ecc (t) S (t) ⁇ exp ( ⁇ i ⁇ ⁇ c (t)) (1)
  • i is an imaginary unit.
  • exp ( ⁇ i ⁇ ⁇ c (t)) will be referred to as a phase correction amount.
  • the phase data calculation unit 940 calculates the phase data ⁇ (t) of the water signal from the measured water signal F (t) according to the following equation (2).
  • ⁇ (t) tan ⁇ 1 (Im (F (t))) / (Re (F (t))) (2)
  • tan ⁇ 1 represents an arctangent function
  • Im (F (t)) represents an imaginary part of the complex number F (t)
  • Re (F (t)) represents a real part of the complex number F (t).
  • phase data correction unit 950 corrects all phase skips existing in the phase data ⁇ (t) by repeated calculation.
  • the phase data correction unit 950 performs a phase return connection process on the phase data ⁇ (t) obtained by the above equation (2) prior to the execution of the phase data correction process.
  • the phase loopback connection process will be described with reference to FIG.
  • the horizontal axis represents the time (ms) from the start of water signal measurement
  • the vertical axis represents the phase value (rad).
  • the broken line indicates the phase data ⁇ (t) calculated from the measurement result
  • solid line indicates the phase data ⁇ z (t) after the phase return connection process.
  • the phase data ⁇ (t) is calculated as a value between ⁇ and + ⁇ .
  • the phase of the water signal acquired by the MRSI pulse sequence 300 has a value exceeding the range of ⁇ to + ⁇ .
  • it is folded back to obtain a value between ⁇ and + ⁇ .
  • phase loop connection processing for removing such a temporal change in the discontinuous phase value is performed to obtain phase data ⁇ z (t) that is the original phase change state.
  • the phase return connection process can use various existing phase return connection processes. At this time, smoothing may be appropriately performed in order to prevent phase variations due to noise.
  • phase data correction process for the phase data ⁇ z (t) after the phase loop back connection process by the phase data correction unit 950 will be described.
  • the region where the phase skip occurs includes a point where the primary time differential value of the phase data is convex upward or downward and the secondary time differential value is zero.
  • the change amount of the primary time differential value of the phase data is simply referred to as “phase change amount”, and the point at which the secondary time differential value is 0 is simply referred to as “inflection point”.
  • the phase data correction unit 950 identifies a region in the phase data ⁇ z (t) where a phase jump occurs based on the magnitude of the phase change amount, and is identified as a phase jump occurring. Phase jump correction is performed only within the specified area.
  • the phase data correction unit 950 divides the phase data ⁇ z (t) into a plurality of small areas in the time direction, and calculates a phase change amount within the small areas. Then, the phase jump is corrected in order from the region where the amount of change is large. At this time, every time correction is performed, the influence of the phase correction amount on the value obtained by Fourier transform is evaluated. When the influence on this value is almost eliminated, it is determined that the correction of the portion where the phase jump has occurred is completed.
  • the absolute value of the difference between the maximum value and the minimum value of the primary differential ⁇ z ′ (t) of the phase data ⁇ z (t) in the divided small region is calculated as the phase change amount.
  • FIG. 6 is a processing flow of phase jump correction processing of the phase data ⁇ z (t) of the present embodiment by the phase data correction unit 950.
  • the area dividing unit 953 performs area dividing processing for dividing the phase data ⁇ z (t) into a plurality of small areas in the time direction (step S601).
  • phase jump generation area extraction unit 951 calculates a phase change amount for each small area divided in step S601 (step S602). Then, the phase jump generation area extraction unit 951 associates each small area with the calculated phase change amount, and holds it in a storage device or the like (step S603).
  • phase jump generation area extraction unit 951 extracts, as a phase jump generation area, a small area held in association with the largest phase change quantity from among the phase change quantities held in the storage device (step). S604).
  • phase jump correction unit 952 determines a phase jump correction area in the phase jump generation area extracted in step S604 (step S605). Thereafter, the phase skip correction unit 952 corrects the phase skip in the phase skip correction area (step S606).
  • the phase jump generation area extraction unit 951 calculates an evaluation value according to a predetermined evaluation function (step S607).
  • the phase jump generation area extraction unit 951 determines whether or not the calculated evaluation value satisfies a predetermined end condition (step S608). If the end condition is not satisfied, it is assumed that the phase jump generation region still remains, and the phase change amount held in the storage device corresponding to the small region corrected in step S606 is set to 0 (step S609).
  • the process returns to S604.
  • the phase data correction unit 950 ends the process, and finally obtains the phase data ⁇ c (t) after the phase jump correction.
  • the process in step S609 is not limited to the above, and any process may be used as long as the information stored in step S603 is updated so that the processed small area can be excluded from the subsequent processes.
  • the area dividing unit 953 calculates an inflection point of the phase data ⁇ z (t), divides the divided area so that one inflection point is included, and obtains each small area.
  • FIG. 7 is a diagram for explaining the region division processing of the present embodiment by the region dividing unit 953, (a) is a processing flow, and (b) and (c) are explanatory diagrams of the division method. .
  • the region dividing unit 953 calculates the second derivative ⁇ z ′′ (t) of the phase data ⁇ z (t) of the water signal, and the point where the value becomes 0 is shown in FIG. 7 is extracted as an inflection point 810 of the phase data ⁇ z (t) of the signal (step S801)
  • the horizontal axis represents time (ms) from the start of measurement of the water signal, and the vertical axis.
  • tau a second-order differential value
  • the region dividing unit 953 identifies the time t IPn (n is any one of 1, 2, 3, 4, 5) of each inflection point 810 (step S802), and sequentially sets the boundary points 811 in the time direction. Determine (step S803).
  • the region dividing unit 953 picks up each inflection point 810 in order in the time direction, and the first inflection point 810 and the last inflection point 810 are respectively the nearest data end point (t 0 or t F ).
  • an intermediate time t IPM1, t IPM6 between, in the case of other inflection point 810, inflection point 810 (t IPn and t IPn-1) the intermediate time t IPMN between, each boundary point 811 Time.
  • the region dividing unit 953 changes the time region of the phase data ⁇ z (t) of the water signal by the boundary point 811 as shown in FIG. 7C, for example, small regions 821, 822, 823, 824, 825. , 826, 827.
  • a horizontal axis is time (ms) from the measurement start of a water signal
  • shaft is a phase value (rad) of phase data.
  • the area dividing unit 953 associates each area divided by the above procedure, for example, each of the small areas 821, 822, 823, 824, 825, 826, and 827 with each end.
  • time boundary point 811 (t 0, t IPM1) , (t IPM1, t IPM2), (t IPM2, t IPM3), (t IPM3, t IPM4), (t IPM4, t IPM5), (t IPM5, t IPM6 ), (t IPM6 , t F ) are stored.
  • FIG. 8 is a graph 831 of the first derivative ⁇ z ′ (t) of the phase data ⁇ z (t), the horizontal axis is the time t (ms) from the start of measurement of the water signal, and the vertical axis is the first derivative. (Au).
  • the phase jump generation area extraction unit 951 first obtains the first derivative ⁇ z ′ (t) of the phase data ⁇ z (t), and for each small area obtained by the division in step S601, the corresponding area The absolute value of the difference between the maximum value and the minimum value of the primary differential ⁇ z ′ (t) is calculated. Then, the phase jump generation area extraction unit 951 sets the calculated absolute value as the phase change amount.
  • phase change amount calculation processing for calculating the phase change amount is not limited to this.
  • the first derivative ⁇ z ′ (t) of the phase data ⁇ z (t) is obtained, and the first derivative ⁇ z ′ (t) is fitted with a model function such as Gaussian for each small region obtained by the division in step S601.
  • the height of the fitting function may be the phase change amount.
  • FIG. 9 is a first-order differential graph 831 in the vicinity of the small region determined as the phase jump generation region, where the horizontal axis represents time (ms) and the vertical axis represents the first-order differential value (au).
  • the small region 822 is determined as the phase jump generation region 841.
  • the phase jump correction unit 952 calculates the peak half-value width of the first-order differential ⁇ z ′ (t) in the phase jump generation area 841, and sets the double width as the phase jump correction area 842.
  • the method of determining the phase skip correction area 1002 is not limited to this. Any region may be used as long as the phase change amount in the small region is minimized. For example, when one point is excluded from both sides of the primary differential peak point corresponding to the inflection point, the phase change amount is minimized. You may decide to become. Further, in the phase change amount calculation process as described above, when the phase skip generation region 841 is determined by fitting with a model function such as Gaussian, the phase skip correction region 842 is set to twice the half-value width of the fitting function. It may be configured.
  • the phase skip correction unit 952 performs linear interpolation that connects both ends of the phase skip correction region 842 with a straight line on the graph of the primary differential ⁇ z ′ (t), for example, as shown in FIG. To determine the corrected primary differential ⁇ c ′ (t). Then, the phase jump correction unit 952 calculates the phase value ⁇ c (t) after the phase jump correction by integrating the corrected first derivative ⁇ c ′ (t).
  • the phase difference ⁇ c (t) after phase jump correction may be calculated from the corrected primary differential ⁇ c ′ (t), not from the integral, but from the sum of powers.
  • FIG. 10A shows an example of a graph of phase data before and after the phase jump correction process.
  • a broken line is a graph of the phase data ⁇ z (t) of the water signal before phase jump correction
  • a solid line is a graph of the water signal phase data ⁇ c (t) after correcting the phase jump of the phase jump generation region 841 (small region 822). is there.
  • the horizontal axis represents time (ms) and the vertical axis represents phase (rad).
  • ms time
  • rad phase
  • phase skip correction unit 952 interpolates the phase skip correction region r so that the primary differential ⁇ ′ (t) becomes a smooth curve on both sides of the phase skip correction region r, and the corrected primary differential ⁇ c ′ ( t) may be determined.
  • the phase jump generation area extraction unit 951 calculates the following equation (3) with respect to the phase data ⁇ n (t) after the n-th phase jump correction.
  • g n (t) exp (i ⁇ ⁇ n (t)) (3)
  • the phase jump generation region extraction unit 951 calculates G n (f) obtained by Fourier transform of g n (t) (formula (4)).
  • G n (f) F [gn (t)] (4)
  • F [] indicates Fourier transform.
  • Phase jump occurs region extraction unit 951 the absolute value of the calculated G n (f)
  • a function, and ⁇ P n obtained by calculating the evaluation function is an evaluation value.
  • the obtained evaluation value is less than a threshold value.
  • a threshold value a value close to 0 is set.
  • after phase jump correction approaches 0 means that the effect does not change even if a new area is extracted and phase jump correction is performed. Means. This is because the phase change other than the phase jump due to eddy current or noise is not so large as to significantly change the evaluation value. Therefore, in this embodiment, when this end condition is satisfied, it is determined that the phase skip correction has been completed for all the portions where the phase skip has occurred, and the process is ended.
  • FIG. 10B shows an example of a phase data graph before the phase skip correction processing and a phase data graph after all the phase skip regions have been corrected.
  • a broken line is a graph of the phase data ⁇ z (t) of the water signal before the phase jump correction process of the present embodiment.
  • the solid line is a graph of the phase data ⁇ c (t) of the water signal after correcting all the phase jumps.
  • the horizontal axis represents time (ms) and the vertical axis represents phase (rad).
  • the phase data of the eddy current correction signal is divided into a plurality of regions in the time direction, and the region having the largest phase change amount is extracted.
  • the phase skip of the region is corrected, and the evaluation value is calculated based on the evaluation function. This process is repeated, and when the evaluation value approaches 0, the process ends.
  • a phase jump having a plurality of phase changes is automatically determined and removed. Further, eddy current correction is performed using the phase data after removal.
  • the phase jump generation region is extracted until the evaluation value approaches 0, and the phase jump correction is repeated. Therefore, the phase jump is performed at a place where the phase data value is irregularly changed. Only the location can be extracted and corrected. Therefore, it is possible to remove the phase jump of the phase data efficiently and accurately. Thereby, in eddy current correction, phase data having no phase jump can be used, and a good spectrum without ringing artifacts can be acquired.
  • FIG. 11 shows the computer simulation results.
  • FIG. 11A is a metabolite spectrum when eddy current correction is not performed
  • FIG. 11B is a case where eddy current correction is performed using the phase data that is not subjected to phase jump correction of this embodiment.
  • (C) shows the metabolite spectrum when the eddy current correction is performed using the phase data subjected to the phase jump correction of the present embodiment.
  • FIG. 12A shows a metabolite spectrum when eddy current correction is not performed
  • FIG. 12B shows a metabolite when eddy current correction is performed using phase data that is not subjected to phase jump correction according to this embodiment.
  • Spectrum (c) shows a metabolite spectrum when eddy current correction is performed using the phase data subjected to the phase skip correction of the present embodiment.
  • phase data used for correcting the eddy current distortion includes a plurality of locations where the phase change amount varies greatly, such as a large amount of phase change. Since the phase skip correction can be performed by extracting only the phase skip portion that causes the ringing artifact, the phase skip of the phase data can be corrected efficiently. Therefore, in correcting the spectral distortion due to the eddy current, the generation of ringing artifacts can be effectively prevented, and the spectral distortion due to the eddy current can be corrected well.
  • the case where the phase data used for eddy current correction is obtained from the water signal is taken as an example, but the present invention is not limited to this.
  • the signal is not particularly limited to a water signal as long as it is a signal of a substance having a signal intensity greater than that of a metabolite to be measured.
  • the phase jump generation region extraction unit 951 calculates the phase change amount for all the small regions and holds it in the storage device. I can't.
  • the phase skip correction process may be configured to re-divide the region for each phase skip correction and calculate a small region having the largest phase change amount. That is, after step S602, the result is not retained, the phase jump generation region in step S604 is extracted, and the process proceeds with steps S605, S606, S607, and S608. If the condition is not satisfied, the corrected result is obtained in step S606. The process is repeated from step S601 for the phase data.
  • the evaluation function is used to evaluate the influence of the phase jump correction on the eddy current correction coefficient, and the phase jump generation region is repeatedly extracted.
  • a threshold value may be set in advance, a region where the phase change amount is equal to or greater than the threshold value may be extracted as a phase skip generation region, and the phase skip correction process may be performed on all the extracted regions.
  • the phase jump is a phase change caused by a vector sum of two vectors having different intensities and frequencies. Therefore, the closer the intensity of the two vectors, the more rapidly the phase changes in a short time width in which the two vectors face each other. On the contrary, as the intensity of the two vectors is different, the influence of the vector having a smaller intensity is reduced, so that the phase gradually changes in the opposing time width. Therefore, the magnitude of the phase change amount and the time width (change width) of the phase skip correction region are in an inversely proportional relationship.
  • phase skip correction is not performed in order from the region having the largest phase change amount, but the phase skip correction is performed in order from the region having the smaller change width until the influence of the phase jump correction is eliminated. Good.
  • a smoothing process is performed to reduce the influence of noise.
  • the magnitude of the phase change amount and the above-described change width may be used in combination. For example, when a threshold value is set in advance and the phase change amount is equal to or greater than a predetermined first threshold value and the change width is equal to or less than a predetermined second threshold value, it is determined as a phase jump occurrence region. The phase jump correction may be performed.

Landscapes

  • Physics & Mathematics (AREA)
  • High Energy & Nuclear Physics (AREA)
  • Condensed Matter Physics & Semiconductors (AREA)
  • General Physics & Mathematics (AREA)
  • Magnetic Resonance Imaging Apparatus (AREA)

Abstract

　ＭＲＩ装置において、例えば、ＭＲスペクトロスコピー、ＭＲスペクトロスコピックイメージングといった計測において、渦電流の影響を補正する処理を行うとき、補正処理によるアーチファクトを防ぎ、補正結果の精度を向上させることを目的とする。　計測対象の代謝物質より信号強度の大きい物質の信号の位相データを用いて渦電流の影響を補正する渦電流補正処理において、位相変化量の大きさにより、位相データ内の位相とびが発生している領域箇所を特定し、位相とびと特定された箇所のみ、位相とび補正を行う。変化量の変動の大きさは、予め定められた領域分割に従って分割された各領域内で決定する。

Description

磁気共鳴撮影装置
　本発明は、スペクトロスコピックイメージングに関し、特に渦電流などによるスペクトル歪みを補正する渦電流補正の精度を高める技術に関する。
　ＭＲＩで計測される核磁気共鳴信号には、分子構造の違いによって、共鳴周波数が僅かに異なる化学シフト現象が生じる。この現象を利用して、分子（代謝物質）ごとに核磁気共鳴信号を分離し、スペクトルを得るＭＲＳ（Ｍａｇｎｅｔｉｃ　Ｒｅｓｏｎａｎｃｅ　Ｓｐｅｃｔｒｏｓｃｏｐｙ）、代謝物質ごとの空間的な信号強度分布を画像化するＣＳＩ（Ｃｈｅｍｉｃａｌ　Ｓｈｉｆｔ　Ｉｍａｇｉｎｇ：化学シフトイメージング）、ＭＲＳＩ（Ｍａｇｎｅｔｉｃ　Ｒｅｓｏｎａｎｃｅ　Ｓｐｅｃｔｒｏｓｃｏｐｉｃ　Ｉｍａｇｉｎｇ：核磁気共鳴スペクトロスコピックイメージング）が知られている。
　ＭＲＳやＭＲＳＩで検出できる人体の主な代謝物質には、コリン（Ｃｈｏ）、クレアチン（Ｃｒ）、Ｎ－アセチルアスパラギン酸（ＮＡＡ）、乳酸（Ｌａｃ）等がある。これら代謝物質の量から、がん等の代謝異常疾患の進行度判定や早期診断が可能となる。また、腫瘍の悪性度診断を非侵襲的に行うことが可能になると考えられる。
　ＭＲＳおよびＭＲＳＩでは、計測時に印加する傾斜磁場によって渦電流が生じる。渦電流は、空間的、時間的に静磁場不均一を生じさせ、計測によって得られたスペクトルの形状を歪ませる。このスペクトル歪みは、一般に参照信号の位相値を用いて補正される。例えば、参照信号として代謝物質よりも信号強度が大きい水信号を用いる。この水信号から、位相値を取得し、補正に使用する（例えば、非特許文献１参照。）。非特許文献１に開示の手法では、この水信号から空間的、時間的な位相値を算出し、代謝物質画像データに対して位相補正を行い、渦電流によるスペクトル歪みを補正している。
　ＣＳＩあるいはＭＲＳＩは、計測時間およびＳＮＲ（信号対雑音比）の観点から、計測するマトリクス（ボクセル）数が８×８～３２×３２程度と非常に少ない。このため、画像再構成で実施するフーリエ変換によって、トランケーションが生じ、遠くのボクセルの信号が混入する。その結果、静磁場不均一が存在する場合、対象とするボクセル内の水信号と異なる周波数の水信号が混入する。
　異なる周波数の水信号が混入すると、渦電流補正に用いる位相値の時間変化に、単位時間当たりの位相値の変化量の変動が他の箇所に比べて突出して大きい、いわゆる、位相とび領域が発生する。この位相とび領域における位相変化量の大きさは、混入した水信号の濃度に比例する。ここで、図１３を用いて、この位相とびについて説明する。
　図１３（ａ）に、計算機シミュレーションによる、水信号のスペクトルを示す。図１３（ａ）におけるＷａｔｅｒ１およびＷａｔｅｒ２は、周波数がそれぞれ２Ｈｚと５Ｈｚ、濃度比が１．０対０．９の水信号のスペクトルである。この２つの信号が混在する場合の、時間領域における位相値の時間変化を図１３（ｂ）に示す。図１３（ａ）、（ｂ）に示すように、周波数差Δｆに対して、１／Δｆの時間間隔で、濃度比に比例した位相とびが生じる。なお、ここでは、渦電流による静磁場変動を付加し、緩やかな位相変化を生じさせている。
　ここで、時間方向に図１３（ｂ）に示す変化を示す位相データを用いて、代謝物質の渦電流によるスペクトル歪みを補正する場合の、補正前後の様子を、それぞれ図１３（ｃ）、（ｄ）に示す。図１３（ｃ）には、渦電流補正前の様子を、図１３（ｄ）は渦電流補正後の様子をそれぞれ示す。図１３（ｄ）に示すように、位相とびのある位相データを用いて補正すると、位相とびによって、リンギングアーチファクトが生じ、渦電流補正処理によって逆にスペクトルが劣化する。
　これに対し、例えば、水信号のスペクトルに対してローパスフィルタをかけることによって、リンギングアーチファクトを低減する手法がある（例えば、非特許文献２参照。）。また、位相データに現れる位相とびを補正し、リンギングアーチファクトを低減する手法がある（例えば、非特許文献３参照。）。ここでは、時間領域における水信号の絶対値強度が極値を取るタイミングを位相とびの発生箇所とする。そして、時間領域における水信号の位相データに関して時間ｔについての１次導関数を算出し、上記の位相とび発生箇所の周辺をモデル関数でフィッティングする。このとき、フィッティングした領域を位相とび補正領域とし、上記位相値を線形補間することにより位相とびを補正する。
Ｕｗｅ　Ｋｌｏｓｅ　"Ｉｎ　Ｖｉｖｏ　Ｐｒｏｔｏｎ　Ｓｐｅｃｔｒｏｓｃｏｐｙ　ｉｎ　Ｐｒｅｓｅｎｃｅ　ｏｆ　Ｅｄｄｙ　Ｃｕｒｒｅｎｔｓ"　Ｍａｇｎｅｔｉｃ　Ｒｅｓｏｎａｎｃｅ　Ｉｎ　Ｍｅｄｉｃｉｎｅ、１４巻、２６－３０頁（１９９０年） Ｊ．　Ｍ．　Ｗｉｌｄ　"Ａｒｔｉｆａｃｔｓ　Ｉｎｔｒｏｄｕｃｅｄ　ｂｙ　Ｚｅｒｏ　Ｏｒｄｅｒ　Ｐｈａｓｅ　Ｃｏｒｒｅｃｔｉｏｎ　ｉｎ　Ｐｒｏｔｏｎ　ＮＭＲ　Ｓｐｅｃｔｒｏｓｃｏｐｙ　ａｎｄ　ａ　Ｍｅｔｈｏｄ　ｏｆ　Ｅｌｉｍｉｎａｔｉｏｎ　ｂｙ　Ｐｈａｓｅ　Ｆｉｌｔｅｒｉｎｇ"　Ｊｏｕｒｎａｌ　ｏｆ　Ｍａｇｎｅｔｉｃ　Ｒｅｓｏｎａｎｃｅ、１３７巻、４３０－４３６頁（１９９９年） Ａ．　Ｗ．　Ｓｉｍｏｎｅｔｔｉ，ｅｔ　ａｌ．"Ａｕｔｏｍａｔｅｄ　ｃｏｒｒｅｃｔｉｏｎ　ｏｆ　ｕｎｗａｎｔｅｄ　ｐｈａｓｅ　ｊｕｍｐｓ　ｉｎ　ｒｅｆｅｒｅｎｃｅ　ｓｉｇｎａｌｓ　ｗｈｉｃｈ　ｃｏｒｒｕｐｔ　ＭＲＳＩ　ｓｐｅｃｔｒａ　ａｆｔｅｒ　ｅｄｄｙ　ｃｕｒｒｅｎｔ　ｃｏｒｒｅｃｔｉｏｎ"　Ｊｏｕｒｎａｌ　ｏｆ　Ｍａｇｎｅｔｉｃ　Ｒｅｓｏｎａｎｃｅ、１５９巻、１５１－１５７頁（２００２年）
　しかしながら、非特許文献２に記載の手法では、位相とびの低周波成分が残るため、リンギングアーチファクトを完全には除去できない。また、渦電流の高周波成分をカットしてしまうため、十分な渦電流の補正効果が得られない。
　非特許文献３に記載の手法では、時間領域における水信号の絶対値強度に極値が存在しなくても位相とびが発生するときには、位相とびが補正されない。また、位相変化が急峻な箇所が複数ある場合、その中から位相とびの箇所を特定し、抽出することが難しい。さらに、フィッティング領域を指定していないため、位相変化量が異なる複数の位相とびがある場合、フィッティング精度が低下する。すなわち、フィッティング精度を高めるためには、位相とび発生箇所毎にフィッティングを行う必要があり、処理が煩雑となるとともに、自動的に位相とびを補正することが難しい。
　本発明は、上記事情に鑑みてなされたものであり、ＭＲＩ装置における渦電流の影響を補正する処理を行うとき、補正処理によるアーチファクトを防ぎ、補正結果の精度を向上させることを目的とする。
　本発明は、計測対象の代謝物質より信号強度の大きい物質の信号の位相データを用いて渦電流の影響を補正する渦電流補正処理において、時間領域における位相変化量の大きさにより、位相データ内の位相とびが発生している領域箇所を特定し、位相とびと特定された箇所のみ、位相とび補正を行う。変化量の大きさは、予め定められた領域分割に従って分割された各領域内で決定する。
　具体的には、被検体に静磁場を印加する静磁場印加手段と、前記被検体に傾斜磁場を印加する傾斜磁場印加手段と、前記被検体に高周波磁場パルスを照射する高周波磁場パルス照射手段と、前記被検体から核磁気共鳴信号を受信する受信手段と、前記核磁気共鳴信号に渦電流補正を実施する渦電流補正手段と、を備える磁気共鳴撮影装置であって、前記渦電流補正手段は、計測対象の代謝物質より信号強度の大きい物質の核磁気共鳴信号の位相データを算出する位相データ算出手段と、前記位相データを、時間方向に複数の小領域に分割する位相データ分割手段と、前記小領域の中から位相とびが存在する位相とび発生領域を抽出する位相とび発生領域抽出手段と、前記位相とび発生領域内の位相とびを補正する位相とび補正手段と、を備えることを特徴とする磁気共鳴撮影装置を提供する。
　本発明によれば、ＭＲＩ装置における渦電流の影響を補正する処理を行うとき、補正処理によるアーチファクトを防ぎ、補正結果の精度を向上させることができる。
（ａ）～（ｃ）は、本発明の実施形態の磁気共鳴撮影装置の外観図である。 （ａ）は、本発明の実施形態の磁気共鳴撮影装置の機能構成図であり、（ｂ）は、本発明の実施形態の磁気共鳴撮影装置が備える計算機の機能ブロック図である。 本発明の実施形態のＭＲＳＩパルスシーケンスの一例を示す図である。 （ａ）～（ｃ）は、本発明の実施形態のＭＲＳＩパルスシーケンスにより励起される領域を説明するための図である。 本発明の実施形態の位相折り返し接続処理の一例を説明するための図である。 本発明の実施形態の位相とび補正処理のフローチャートである。 （ａ）は、本発明の実施形態の領域分割処理のフローチャートであり、（ｂ）および（ｃ）は、本発明の実施形態の領域分割処理を説明するための図である。 本発明の実施形態の位相変化量算出処理を説明するための図である。 本発明の実施形態の位相とび補正領域決定処理および位相とび補正処理を説明するための図である。 （ａ）は、本発明の実施形態の位相とび補正処理を実施する前の位相データと、位相とび領域を１箇所補正した後の位相データである。（ｂ）は、本発明の実施形態の位相とび補正処理を実施する前の位相データと、全ての位相とび領域を補正した後の位相データを示す図である。 （ａ）～（ｃ）は、計算機シミュレーション結果を用いて本発明の実施形態による効果を説明するための図である。 （ａ）～（ｃ）は、ファントム実験結果を用いて本発明の実施形態による効果を説明するための図である。 （ａ）～（ｄ）は、水信号の位相データを用いて渦電流補正を行う場合、リンギングアーチファクトが生じる原因を説明するための図である。
　以下、本発明を適用する実施形態について説明する。以下、実施形態を説明するための全図において、同一機能を有するものは同一符号を付し、その繰り返しの説明は省略する。
　まず、本実施形態の磁気共鳴撮影装置（ＭＲＩ装置）について説明する。図１は、本実施形態のＭＲＩ装置の外観図である。図１（ａ）は、ソレノイドコイルで静磁場を生成するトンネル型磁石を用いた水平磁場方式のＭＲＩ装置１００である。図１（ｂ）は、開放感を高めるために磁石を上下に分離したハンバーガー型（オープン型）の垂直磁場方式のＭＲＩ装置１１０である。また、図１（ｃ）は、図１（ａ）と同じトンネル型磁石を用い、磁石の奥行を短くし且つ斜めに傾けることによって、開放感を高めたＭＲＩ装置１２０である。本実施形態では、これらの外観を有するＭＲＩ装置のいずれを用いることもできる。なお、これらは一例であり、本実施形態のＭＲＩ装置はこれらの形態に限定されるものではない。本実施形態では、装置の形態やタイプを問わず、公知の各種のＭＲＩ装置を用いることができる。以下、特に区別する必要がない場合は、ＭＲＩ装置１００で代表する。
　図２（ａ）は、本実施形態のＭＲＩ装置１００の機能構成図である。本図に示すように、本実施形態のＭＲＩ装置１００は、被検体１が置かれる空間に、静磁場を生成する静磁場コイル２と、ｘ方向、ｙ方向、ｚ方向にそれぞれ傾斜磁場を発生させる傾斜磁場コイル３と、静磁場分布を調整するシムコイル４と、被検体１の計測領域に対し高周波磁場を照射する計測用高周波コイル５（以下、単に送信コイルという）と、被検体１から発生する核磁気共鳴信号を受信する受信用高周波コイル６（以下、単に受信コイルという）とを備える。
　静磁場コイル２は、図１に示した各ＭＲＩ装置１００、１１０、１２０の構造に応じて、種々の形態のものが採用される。傾斜磁場コイル３及びシムコイル４は、それぞれ傾斜磁場用電源部１２及びシム用電源部１３により駆動される。なお、本実施形態では、送信コイル５と受信コイル６とに別個のものを用いる場合を例にあげて説明するが、送信コイル５と受信コイル６との機能を兼用する１のコイルで構成してもよい。送信コイル５が照射する高周波磁場は、送信機７により生成される。受信コイル６が検出した核磁気共鳴信号は、受信機８を通して計算機９に送られる。
　シーケンス制御装置１４は、傾斜磁場発生コイル３の駆動用電源である傾斜磁場用電源部１２、シムコイル４の駆動用電源であるシム用電源部１３、送信機７及び受信機８の動作を制御し、傾斜磁場、高周波磁場の印加および核磁気共鳴信号の受信のタイミングを制御する。制御のタイムチャートはパルスシーケンスと呼ばれ、計測に応じて予め設定され、後述する計算機９が備える記憶装置等に格納される。
　計算機９は、受け取った核磁気共鳴信号に対して様々な演算処理を行い、画像情報やスペクトル情報を生成するとともに、ＭＲＩ装置１００全体の動作を制御する。計算機９は、ＣＰＵ、メモリ、記憶装置などを備える情報処理装置であり、計算機９にはディスプレイ１０、外部記憶装置１１、入力装置１５などが接続される。ディスプレイ１０は、演算処理で得られた結果等をオペレータに表示するインタフェースである。入力装置１５は、本実施形態で行われる演算処理に必要な条件、パラメータ等をオペレータが入力するためのインタフェースである。外部記憶装置１１は、計算機９が備える記憶装置とともに、計算機９が実行する各種の演算処理に用いられるデータ、演算処理により得られるデータ、入力された条件、パラメータ等を保持する。
　以下、本実施形態の計算機９が実現する機能について説明する。図２（ｂ）は、本実施形態の計算機９の機能ブロック図である。本実施形態の計算機９は、計測制御部９１０と、表示情報生成部９２０と、渦電流補正部９３０とを備える。計測制御部９１０は、パルスシーケンスに従ってシーケンス制御装置１４を動作させるとともに各部を制御し、計測を行う。表示情報生成部９２０は、計測により得られた核磁気共鳴信号に対して様々な演算処理を行い画像情報やスペクトル情報を生成する。
　渦電流補正部９３０は、渦電流によるスペクトル歪みを補正する渦電流補正処理を行う。本実施形態では、渦電流によるスペクトル歪みを、代謝物質よりも信号強度の大きい水信号の位相データにより補正する。これを実現するため、渦電流補正部９３０は、得られた核磁気共鳴信号から水信号の位相データを算出する位相データ算出部９４０と、位相データ内の全ての位相とびを繰り返し計算によって補正する位相データ補正部９５０とを備える。さらに、本実施形態の位相データ補正部９５０は、位相データを時間方向に複数の小領域に分割する領域分割部９５３と、分割された小領域の中から位相とびが発生している領域を抽出する位相とび発生領域抽出部９５１と、前記位相とび発生領域抽出部９５１が抽出した小領域内の位相とびを補正する位相とび補正部９５２とを備える。
　なお、計算機９が実現する各種の機能は、記憶装置が保持するプログラムを、ＣＰＵがメモリにロードして実行することにより実現される。また、計算機９が実現する各種の機能のうち、少なくとも１の機能は、ＭＲＩ装置１００とは独立した、情報処理装置であって、ＭＲＩ装置１００とデータの送受信が可能な情報処理装置により実現されていてもよい。
　以下、上記各機能による、本実施形態の計測全体の流れについて簡単に説明する。本実施形態では、渦電流補正用信号として水信号を用いる。従って、まず、計測制御部９１０は、予め定められたパルスシーケンスに従ってシーケンス制御装置１４を制御し、非水抑圧計測を行い、水信号Ｆ（ｔ）を得る。なお、ここで、ｔは、サンプリング時間であり、離散値をとる変数である。得られた水信号Ｆ（ｔ）から位相データ算出部９４０は、位相データφ（ｔ）を算出する。そして、位相データ補正部９５０は、位相データφ（ｔ）に存在する全ての位相とびが補正された位相データφｃ（ｔ）を得る。次に、計測制御部９１０は、予め定められたパルスシーケンスに従ってシーケンス制御装置１４を制御し、水抑圧計測を行い、代謝物質信号Ｓ（ｔ）を得る。渦電流補正部９３０は、得られた代謝物質信号Ｓ（ｔ）を位相データφｃ（ｔ）で補正し、渦電流補正後の代謝物質信号Ｓｅｃｃ（ｔ）を得る。表示情報生成部９２０は、代謝物質信号Ｓｅｃｃ（ｔ）をフーリエ変換し、代謝物質のスペクトルまたは分布画像を得る。
　ここで、計測制御部９１０が上記計測に用いるパルスシーケンスの一例について説明する。ここでは、代謝物質を画像化する領域選択型ＭＲＳＩのパルスシーケンス（以降、ＭＲＳＩパルスシーケンスと呼ぶ）を例にあげて説明する。
　図３は、ＭＲＳＩパルスシーケンス３００の一例である。図３において、ＲＦは高周波磁場パルスの印加タイミングを示す。Ｇｘ、Ｇｙ、Ｇｚは、それぞれ、ｘ、ｙ、ｚ方向の傾斜磁場パルスの印加タイミングを示す。Ａ／Ｄは、信号の計測期間Ｔｐ１を示す。図３に示すＭＲＳＩパルスシーケンス３００は、公知のＭＲＳＩパルスシーケンスと同じである。すなわち、１つの励起パルスＲＦ１と２つの反転パルスＲＦ２およびＲＦ３とを用いて、所定の関心領域を選択的に励起し、この関心領域からＦＩＤ信号（自由誘導減衰）ＦＩＤ１を得ることを、繰り返し時間ＴＲ間隔で繰り返す。
　このＭＲＳＩパルスシーケンス３００に従って、励起される領域を図４に示す。図４は、本計測に先立って行われる計測により得る位置決め用スカウト画像であって、それぞれ、図４（ａ）はトランス像、図４（ｂ）はサジタル像、図４（ｃ）はコロナル像である。以下、各部の動作と励起される領域との関係を図３および図４を用いて説明する。
　まず高周波磁場ＲＦ１とz方向の傾斜磁場パルスＧｓ１、Ｇｓ１’とを印加して、ｚ方向の断面４０１を励起する。ＴＥ／４（ここで、ＴＥはエコー時間）時間後に、高周波磁場ＲＦ２とｙ方向の傾斜磁場パルスＧｓ２とを印加する。その結果、ｚ方向の断面４０１とｙ方向の断面４０２とが交差する領域における核磁化の位相のみが反転する。続いて、高周波磁場ＲＦ２印加からＴＥ／２後に高周波磁場ＲＦ３とｘ方向の傾斜磁場パルスＧｓ３とを印加する。それによって、ｚ方向の断面４０１、ｙ方向の断面４０２、ｘ方向の断面４０３が交差する関心領域４０４における核磁化の位相のみが反転し、ここから自由誘導減衰信号ＦＩＤ１が生じる。この自由誘導減衰信号ＦＩＤ１を計測する。なお各方向の傾斜磁場パルスＧｄ１～Ｇｄ３およびＧｄ１’～Ｇｄ３’は、高周波磁場ＲＦ１で励起された核磁化の位相は乱さず、ＲＦ２、ＲＦ３で励起された核磁化の位相をディフェイズするための傾斜磁場である。また、高周波磁場ＲＦ３の後には、位相エンコード傾斜磁場Ｇｐ１、Ｇｐ２を印加する。以上により、関心領域４０４の核磁気共鳴信号を得る。
　次に、渦電流補正部９３０による渦電流補正を説明する。渦電流補正部９３０は、代謝物質信号をＳ（ｔ）から、位相とび補正後の水信号の位相データφｃ（ｔ）を用いて、渦電流補正後の代謝物質信号Ｓｅｃｃ（ｔ）を、以下の式（１）に従って算出する。
　Ｓｅｃｃ（ｔ）＝Ｓ（ｔ）・ｅｘｐ（－ｉ・φｃ（ｔ））　　　（１）
ここで、ｉは虚数単位である。
なお、今後、ｅｘｐ（－ｉ・φｃ（ｔ））を位相補正量と呼ぶ。
　まず、位相データ算出部９４０による位相データの算出処理を説明する。位相データ算出部９４０は、計測された水信号Ｆ（ｔ）から水信号の位相データφ（ｔ）を、以下の式（２）に従って算出する。
　φ（ｔ）＝ｔａｎ－１（Ｉｍ（Ｆ（ｔ）））／（Ｒｅ（Ｆ（ｔ）））　　　（２）
ここで、ｔａｎ－１は、アークタンジェント関数、Ｉｍ（Ｆ（ｔ））は複素数Ｆ（ｔ）の虚部、Ｒｅ（Ｆ（ｔ））は複素数Ｆ（ｔ）の実部を表す。
　次に、位相データ補正部９５０による、位相データφ（ｔ）に存在する全ての位相とびを繰り返し計算により補正する位相データ補正処理について説明する。なお、位相データ補正部９５０は、位相データ補正処理の実施に先立ち、上記式（２）で得られた位相データφ（ｔ）に対し、位相折り返し接続処理を行う。まず、この位相折返し接続処理について図５を用いて説明する。図５において、横軸は、水信号の計測開始からの時間（ｍｓ）、縦軸は位相値（ｒａｄ）である。また、破線は計測結果から算出された位相データφ（ｔ）、実線は位相折返し接続処理後の位相データφｚ（ｔ）を示す。
　位相データφ（ｔ）は、－πから＋πの間の値として算出される。しかし、上記ＭＲＳＩパルスシーケンス３００により取得する水信号の位相には、－πから＋πの範囲を超えた値が生じる。位相にこのような値が発生するタイミングでは、折り返されて－πから＋πの間の値が得られる。本図に示すように、折り返された部分では、値に不連続な状態が生じる。そこで、このような不連続な位相値の時間的な変化を取り除く位相折り返し接続処理を行い、本来の位相変化の状態である位相データφｚ（ｔ）を得る。なお、位相折り返し接続処理は、既存の各種の位相折り返し接続処理を用いることができる。このとき、ノイズによる位相ばらつきを防ぐため、適宜スムージングを行ってもよい。
　次に、位相データ補正部９５０による、位相折返し接続処理後の位相データφｚ（ｔ）に対する位相データ補正処理について説明する。一般に位相とびが発生している箇所は、その近傍で位相値に急激な変化が発生し、そこに変曲点が生じる。すなわち、位相とびが生じる領域は、位相データの１次時間微分値が上あるいは下に凸となり、２次時間微分値が０となる点が含まれる。以降、位相データの１次時間微分値の変化量を単に「位相変化量」、２次時間微分値が０となる点を単に「変曲点」と呼ぶ。ここで、位相値に不規則な変化が発生している箇所であっても、その位相変化量が所定以下の箇所は、位相とびによるものではなく、渦電流の影響など、他の理由による。本実施形態では、位相データ補正部９５０は、位相変化量の大きさにより、位相データφｚ（ｔ）内の位相とびが発生している領域を特定し、位相とびが発生していると特定された領域内でのみ、位相とび補正を行う。
　具体的には、位相データ補正部９５０は、位相データφｚ（ｔ）を時間方向に複数の小領域に分割し、小領域内で位相変化量を算出する。そして、この変化量の大きい領域から順に位相とびを補正する。このとき、補正する毎に、位相補正量をフーリエ変換して得る値に与える影響を評価する。この値に与える影響が略なくなったところで、位相とびが発生している箇所の補正を終えたものと判断する。
　なお、本実施形態では、分割した各小領域内で、当該領域内の位相データφｚ（ｔ）の１次微分φｚ’（ｔ）の最大値と最小値との差分の絶対値を位相変化量とする。
　図６は、位相データ補正部９５０による本実施形態の位相データφｚ（ｔ）の位相とび補正処理の処理フローである。まず、領域分割部９５３は、位相データφｚ（ｔ）を時間方向に複数の小領域に分割する領域分割処理を行う（ステップＳ６０１）。
　次に、位相とび発生領域抽出部９５１は、ステップＳ６０１で分割した小領域毎に、位相変化量を算出する（ステップＳ６０２）。そして、位相とび発生領域抽出部９５１は、各小領域と算出した位相変化量とを対応づけて、記憶装置等に保持する（ステップＳ６０３）。
　次に、位相とび発生領域抽出部９５１は、記憶装置に保持された位相変化量の中から、最も大きい位相変化量に対応づけて保持された小領域を、位相とび発生領域として抽出する（ステップＳ６０４）。
　次に、位相とび補正部９５２は、ステップＳ６０４で抽出した位相とび発生領域において、位相とび補正領域を決定する（ステップＳ６０５）。その後、位相とび補正部９５２は、位相とび補正領域において、位相とびを補正する（ステップＳ６０６）。
　次に、位相とび発生領域抽出部９５１は、予め定められた評価関数に従って評価値を算出する（ステップＳ６０７）。位相とび発生領域抽出部９５１は、算出した評価値が予め定められた終了条件を満たすか否かを判別する（ステップＳ６０８）。終了条件を満たさない場合、位相とび発生領域がまだ残っているものとして、ステップＳ６０６で補正を行った小領域に対応して記憶装置に保持された位相変化量を０とし（ステップＳ６０９）、ステップＳ６０４に戻る。一方、終了条件を満たす場合、位相データ補正部９５０は、処理を終了し、最終的に位相とび補正後の位相データφｃ（ｔ）を得る。なお、ステップＳ６０９の処理は、上記に限られず、処理済みの小領域を、以降の処理から除くことができるようステップＳ６０３で記憶した情報を更新するものであればよい。
　ここで、上記ステップＳ６０１の領域分割処理の詳細について説明する。上述したように、位相とびが生じる領域は変曲点が含まれる。したがって、本実施形態では、領域分割部９５３は、位相データφｚ（ｔ）の変曲点を算出し、分割領域に変曲点が１つ含まれるように分割し、各小領域を得る。図７は、領域分割部９５３による、本実施形態の領域分割処理を説明するための図であり、（ａ）は、処理フロー、（ｂ）および（ｃ）は、分割手法の説明図である。
　まず、領域分割部９５３は、水信号の位相データφｚ（ｔ）の２次微分φｚ”（ｔ）を計算し、その値が０となる点を、図７（ｂ）に示すように、水信号の位相データφｚ（ｔ）の変曲点８１０として抽出する（ステップＳ８０１）。なお、図７（ｂ）のグラフにおいて、横軸は、水信号の計測開始からの時間（ｍｓ）、縦軸は２次微分値（ａ．ｕ）である。ここでは、計測期間（ｔ０～ｔＦ）中に５個の変曲点８１０を抽出した場合を例示する。
　次に、領域分割部９５３は、各変曲点８１０の時間ｔＩＰｎ（ｎは１、２、３、４、５のいずれか）を特定し（ステップＳ８０２）、時間方向に順に境界点８１１を決定する（ステップＳ８０３）。ここでは、領域分割部９５３は、各変曲点８１０を時間方向に順に取り上げ、最初の変曲点８１０および最後の変曲点８１０については、それぞれ最寄りのデータ端点（ｔ０またはｔＦ）との間の中間の時間ｔＩＰＭ１、ｔＩＰＭ６を、他の変曲点８１０の場合は、変曲点８１０（ｔＩＰｎとｔＩＰｎ－１）間の中間の時間ｔＩＰＭｎを、それぞれ境界点８１１の時間とする。そして、領域分割部９５３は、境界点８１１により、水信号の位相データφｚ（ｔ）の時間領域を、図７（ｃ）に示すように、例えば、小領域８２１、８２２、８２３、８２４、８２５、８２６、８２７に分割する。ここで、図７（ｃ）のグラフにおいて、横軸は、水信号の計測開始からの時間（ｍｓ）、縦軸は位相データの位相値（ｒａｄ）である。
　本実施形態では、領域分割部９５３は、以上の手順で分割した各領域を特定する情報、例えば、各小領域８２１、８２２、８２３、８２４、８２５、８２６、８２７に対応づけて、それぞれの両端の境界点８１１の時間（ｔ０、ｔＩＰＭ１）、（ｔＩＰＭ１、ｔＩＰＭ２）、（ｔＩＰＭ２、ｔＩＰＭ３）、（ｔＩＰＭ３、ｔＩＰＭ４）、（ｔＩＰＭ４、ｔＩＰＭ５）、（ｔＩＰＭ５、ｔＩＰＭ６）、（ｔＩＰＭ６、ｔＦ）を記憶する。
　次に、上記ステップＳ６０２における、位相とび発生領域抽出部９５１による位相変化量を算出する位相変化量算出処理について、図８を用いて説明する。図８は、位相データφｚ（ｔ）の、１次微分φｚ’（ｔ）のグラフ８３１であり、横軸は、水信号の計測開始からの時間ｔ（ｍｓ）、縦軸は１次微分値（ａ．ｕ）である。位相変化量算出処理では、位相とび発生領域抽出部９５１は、まず、位相データφｚ（ｔ）の１次微分φｚ’（ｔ）を求め、ステップＳ６０１の分割で得た小領域毎に、当該領域内の１次微分φｚ’（ｔ）の最大値と最小値との差の絶対値を算出する。そして、位相とび発生領域抽出部９５１は、算出した絶対値を、位相変化量とする。
　なお、位相変化量を算出する位相変化量算出処理はこれに限られない。位相データφｚ（ｔ）の１次微分φｚ’（ｔ）を求め、ステップＳ６０１の分割で得た小領域毎に、ガウシアンなどのモデル関数で１次微分φｚ’（ｔ）をフィッティングし、各領域のフィッティング関数の高さを位相変化量としてもよい。
　次に、位相とび補正部９５２による、ステップＳ６０５の位相とび補正領域を決定する補正領域決定処理と、ステップＳ６０６の位相とび補正処理と、について、図９を用いて説明する。図９は、位相とび発生領域として決定された小領域近傍の１次微分グラフ８３１であり、横軸は時間（ｍｓ）、縦軸は１次微分値（ａ．ｕ．）である。ここでは、一例として小領域８２２が、位相とび発生領域８４１と決定されたものとして説明する。
　補正領域決定処理では、位相とび補正部９５２は、位相とび発生領域８４１において、一次微分φｚ’（ｔ）のピーク半値幅を算出し、その２倍の幅を位相とび補正領域８４２とする。なお、位相とび補正領域１００２の決定手法はこれに限られない。当該小領域における位相変化量が、最小となるような領域であればよく、例えば、変曲点に対応する１次微分ピーク点の両側を１点ずつ除外したときに、位相変化量が最小となるように決定してもよい。また、上記のように位相変化量算出処理において、ガウシアンなどのモデル関数でフィッティングし、位相とび発生領域８４１を決定した場合は、フィッティング関数の半値幅の２倍を位相とび補正領域８４２とするよう構成してもよい。
　また、位相とび補正処理では、位相とび補正部９５２は、例えば、図９に示すように、１次微分φｚ’（ｔ）のグラフ上で、位相とび補正領域８４２の両端を直線で結ぶ線形補間により補正後の１次微分φｃ’（ｔ）を決定する。そして、位相とび補正部９５２は、補正後の１次微分φｃ’（ｔ）を積分することで位相とび補正後の位相値φｃ（ｔ）を算出する。ここで、補正後の１次微分φｃ’（ｔ）を、積分ではなく累乗和から位相とび補正後の位相値φｃ（ｔ）を算出してもよい。
　ここで、上記位相とび補正処理の前後の位相データのグラフの一例を図１０（ａ）に示す。破線が位相とび補正前の水信号の位相データφｚ（ｔ）のグラフ、実線が、位相とび発生領域８４１（小領域８２２）の位相とびを補正後の水信号位相データφｃ（ｔ）のグラフである。本図において、横軸は、時間（ｍｓ）、縦軸は、位相（ｒａｄ）を示す。本図に示すように、上記位相とび補正処理により、水信号の位相とびが解消されていることがわかる。
　なお、位相とび補正の手法は、これに限られない。例えば、位相とび補正部９５２は、位相とび補正領域ｒの両側で一次微分φ’（ｔ）が滑らかな曲線になるよう、位相とび補正領域ｒを補間し、補正後の１次微分φｃ’（ｔ）を決定してもよい。
　次に、ステップＳ６０７およびＳ６０８における、位相とび発生領域抽出部９５１による評価処理について説明する。本実施形態では、位相とび発生領域抽出部９５１は、ｎ回目の位相とび補正後の位相データφｎ（ｔ）に関して、以下の式（３）を計算する。
　ｇｎ（ｔ）＝ｅｘｐ（ｉ・φｎ（ｔ））　　　（３）
次に、位相とび発生領域抽出部９５１は、ｇｎ（ｔ）をフーリエ変換したＧｎ（ｆ）を算出する（式（４））。
　Ｇｎ（ｆ）＝Ｆ[ｇｎ（ｔ）]　　　（４）
ここで、Ｆ[]は、フーリエ変換を示す。位相とび発生領域抽出部９５１は、算出したＧｎ（ｆ）の絶対値｜Ｇｎ（ｆ）｜の最大ピークＰｎを算出する。そして、位相とび発生領域抽出部９５１は、最大ピークＰｎと、ｎ－１回目の位相とび補正後の最大ピークＰｎ－１との差分の絶対値｜Ｐｎ－Ｐｎ－１｜を評価関数とし、この評価関数を計算することにより得られるΔＰｎを評価値とする。
　また、本実施形態の評価処理の終了条件は、得られた評価値が閾値未満となることとする。閾値としては、限りなく０に近い値を設定する。位相とび補正後の｜Ｇｎ（ｆ）｜の最大ピークＰｎの変化が０に漸近するということは、これ以上、新たな領域を抽出して位相とび補正を行っても効果が変わらないことを意味する。なぜなら、渦電流やノイズなどに起因する位相とび以外の位相変化は、上記評価値を著しく変化させるほど大きな変化ではないからである。従って、本実施形態では、本終了条件を満たす場合、位相とびが発生している全ての箇所の位相とび補正を終えたものと判断し、処理を終了させる。
　なお、評価関数として、｜Ｇｎ（ｆ）｜の最大ピーク値の位相とび補正毎の変化ではなく、｜Ｇｎ（ｆ）｜の標準偏差の位相とび補正毎の変化を用いてもよい。また、｜Ｇｎ（ｆ）｜が閾値Ｔ以上の値をとる信号値ｆの数Ｎｎの、位相とび補正ごとの変化を用いてもよい。
　図１０（ｂ）に、位相とび補正処理を実施する前の位相データのグラフと、全ての位相とび領域を補正した後の位相データのグラフの一例を示す。破線が本実施形態の位相とび補正処理前の水信号の位相データφｚ（ｔ）のグラフである。そして、実線が、全ての位相とびを補正した後の水信号の位相データφｃ（ｔ）のグラフである。本図において、横軸は、時間（ｍｓ）、縦軸は、位相（ｒａｄ）を示す。本図に示すように、上記繰り返し計算による位相とび補正処理により、水信号の全ての位相とびが解消されていることがわかる。
　以上説明したように、本実施形態によれば、渦電流補正用信号の位相データを時間方向に複数の領域に分割し、その中で位相変化量が最も大きい領域を抽出する。その領域の位相とびを補正し、評価関数を基に評価値を算出する。この処理を繰り返し、評価値が０に近づいたとき、処理を終了する。このように構成することで、複数の位相変化がある位相とびを自動で判別し除去する。また、除去後の位相データを用い、渦電流補正を実施する。
　本実施形態によれば、評価値が０に漸近するまで、位相とび発生領域を抽出し、位相とび補正を繰り返すため、位相データの値に不規則な変化が発生している箇所において、位相とび箇所のみ抽出して補正することができる。従って、効率よく、また、精度よく、位相データの位相とびを除去することができる。これにより、渦電流補正において、位相とびのない位相データを用いることができ、リンギングアーチファクトが生じない良好なスペクトルを取得することができる。
　ここで、本実施形態による渦電流補正の効果を計算機シミュレーションおよびファントム実験にて確認する。
　計算機シミュレーション結果を図１１に示す。ここで、図１１（ａ）は、渦電流補正を行わなかった場合の代謝物質スペクトル、（ｂ）は、本実施形態の位相とび補正を行わない位相データを用いて渦電流補正を行った場合の代謝物質スペクトル、（ｃ）は、本実施形態の位相とび補正を実施した位相データを用いて渦電流補正を行った場合の代謝物質スペクトルをそれぞれ示す。
　これらの図に示すように、本実施形態の位相とび補正後に渦電流補正を行うと、位相とびによるリンギングアーチファクトが除去でき、かつ、渦電流によるスペクトル歪みも補正できることがわかる。
　次に、ファントムを用い、実際に計測した結果を図１２に示す。ここで用いたファントムは、濃度１０ｍＭのアセチルアラニンファントムである。図１２（ａ）は、渦電流補正を行わなかった場合の代謝物質スペクトル、（ｂ）は、本実施形態の位相とび補正を行わない位相データを用いて渦電流補正を行った場合の代謝物質スペクトル、（ｃ）は、本実施形態の位相とび補正を実施した位相データを用いて渦電流補正を行った場合の代謝物質スペクトルをそれぞれ示す。
　これらの図に示すように、ファントムを用いて実際に計測した結果においても、本実施形態の位相とび補正後に渦電流補正を行うと、位相とびによるリンギングアーチファクトが除去でき、かつ、渦電流によるスペクトル歪みも補正できることがわかる。
　以上、説明したように、本実施形態によれば、渦電流歪みを補正するために用いる位相データに、位相変化量が大きく異なる等の不規則な変化をする箇所が複数ある場合であっても、リンギングアーチファクトの原因となる位相とびの箇所のみを抽出して位相とび補正を行うことができるため、位相データの位相とびを効率よく補正することができる。従って、渦電流によるスペクトル歪みの補正にあたり、リングングアーチファクトの発生を効果的に防ぐことができ、渦電流によるスペクトル歪みを良好に補正することができる。
　このため、本実施形態によれば、渦電流の影響を補正する処理を行うとき、補正処理によるアーチファクトを防ぎ、補正結果の精度を高めることができる。例えば、ＭＲＳＩにおいて、リンギングアーチファクトのない、良好な代謝物質スペクトルを得ることができる。
　なお、本実施形態では、渦電流の補正に用いる位相データを水信号から得る場合を例にあげているが、これに限られない。計測対象の代謝物質よりも信号強度の大きい物質の信号であれば、特に水信号に限られない。
　また、上記位相とび補正処理では、領域分割部９５３が領域分割後、位相とび発生領域抽出部９５１が全小領域について、位相変化量を算出し、記憶装置に保持しているが、これに限られない。例えば、位相とび補正処理を、位相とび補正毎に、領域を分割し直し、位相変化量が最も大きい小領域を算出するよう構成してもよい。すなわち、ステップＳ６０２の後に、その結果を保持せず、ステップＳ６０４の位相とび発生領域を抽出し、ステップＳ６０５、Ｓ６０６、Ｓ６０７、Ｓ６０８と処理を進め、条件を満たさない場合、ステップＳ６０６で補正後の位相データについて、ステップＳ６０１から処理を繰り返す。
　さらに、上記実施形態では、評価関数を用い、位相とび補正が渦電流補正係数に与える影響を評価し、位相とび発生領域を繰り返し抽出しているが、これに限られない。例えば、予め閾値を定めておき、位相変化量が閾値以上である領域を位相とび発生領域として抽出し、抽出した全領域に対し、位相とび補正処理を行うよう構成してもよい。
　また、一般に、局所的な時間範囲で考えた場合、位相とびは、強度および周波数の異なる２つのベクトルのベクトル和によって生じる位相変化である。そのため、２つのベクトルの強度が近ければ近いほど、２つのベクトルが対向する短い時間幅で急激に位相が変化する。逆に、２つのベクトルの強度が異なるほど、強度が小さいベクトルの影響が小さくなるため、対向する時間幅で緩やかに位相が変化する。したがって、位相変化量の大小と、位相とび補正領域の時間幅（変化幅）とは、反比例の関係にある。これを利用し、位相変化量の大きい領域から順に位相とびを補正するのではなく、前記変化幅の小さい領域から順に、位相とび補正の影響がなくなるまで、位相とび補正を行うよう構成してもよい。この場合、位相折返し接続処理後、ノイズによる影響を軽減するためスムージング処理を行う。
　また、位相とび補正を行う領域を決定するに当たり、位相変化量の大きさと、上述の変化幅とを併用するよう構成してもよい。例えば、予めそれぞれ閾値を定めておき、位相変化量が予め定めた第一の閾値以上であり、かつ、変化幅が、予め定めた第二の閾値以下である場合、位相とび発生領域と判別し、位相とび補正を行うよう構成してもよい。
１：被検体、２：静磁場コイル、３：傾斜磁場コイル、４：シムコイル、５：送信コイル、６：受信コイル、７：送信機、８：受信機、９：計算機、１０：ディスプレイ、１１：外部記憶装置、１２：傾斜磁場用電源部、１３：シム用電源部、１４：シーケンス制御装置、１５：入力装置、１００：ＭＲＩ装置、１１０：ＭＲＩ装置、１２０：ＭＲＩ装置、３００：ＭＲＳＩパルスシーケンス、４０１：断面、４０２：断面、４０３：断面、４０４：関心領域、８１０：変曲点、８１１：境界点、８２１：小領域、８２２：小領域、８２３：小領域、８２４：小領域、８２５：小領域、８２６：小領域、８２７：小領域、８３１：グラフ、８４１：位相とび発生領域、８４２：位相とび補正領域、９１０：計測制御部、９２０：表示情報生成部、９３０：渦電流補正部、９４０：位相データ算出部、９５０：位相データ補正部、９５１：位相とび発生領域抽出部、９５２：位相とび補正部、９５３：領域分割部

Claims (15)

  1. 　被検体に静磁場を印加する静磁場印加手段と、前記被検体に傾斜磁場を印加する傾斜磁場印加手段と、前記被検体に高周波磁場パルスを照射する高周波磁場パルス照射手段と、前記被検体から核磁気共鳴信号を受信する受信手段と、前記核磁気共鳴信号に渦電流補正を実施する渦電流補正手段と、を備える磁気共鳴撮影装置であって、
    　前記渦電流補正手段は、
    　計測対象の代謝物質より信号強度の大きい物質の核磁気共鳴信号の位相データを算出する位相データ算出手段と、
    　前記位相データを、時間方向に複数の小領域に分割する位相データ分割手段と、
    　前記小領域の中から位相とびが存在する位相とび発生領域を抽出する位相とび発生領域抽出手段と、
    　前記位相とび発生領域内の位相とびを補正する位相とび補正手段と、を備えること
    　を特徴とする磁気共鳴撮影装置。
  2. 　請求項１記載の磁気共鳴撮影装置であって、
    　前記位相とび発生領域抽出手段は、前記位相データの位相変化量が、予め定められた条件を満たす領域を、前記位相とび発生領域として抽出すること
    　を特徴とする磁気共鳴撮影装置。
  3. 　請求項１記載の磁気共鳴撮影装置であって、
    　前記位相とび補正手段は、前記位相とび発生領域抽出手段により位相とび発生領域が抽出される毎に前記位相とびの補正を行い、
    　前記位相とび発生領域抽出手段は、
    　前記位相とび補正手段によって位相とびが補正される毎に、当該補正後の前記位相データの評価値を算出する評価値算出手段と、
    　前記評価値が算出される毎に、前記評価値が所定の閾値以上であるか否かを判別し、閾値以上である場合、前記位相とび発生領域の抽出を継続するよう決定し、閾値未満である場合、前記抽出を終了するよう決定する判別手段と、を備え、
    　前記判別手段が抽出を継続すると判別する間、前記位相データの位相変化量が最も大きい小領域から順に１の前記小領域を前記位相とび発生領域として抽出し、
    　を特徴とする磁気共鳴撮影装置。
  4. 　請求項１記載の磁気共鳴撮影装置であって、
    　前記位相データ分割手段は、
    　前記位相データの変曲点を算出する変曲点算出手段を備え、
    　前記各小領域に１の前記変曲点算出手段が算出した変曲点が含まれるよう分割すること
    　を特徴とする磁気共鳴撮影装置。
  5. 　請求項３記載の磁気共鳴撮影装置であって、
    　前記位相とび発生領域抽出手段は、
    　前記位相データの１次微分値を算出する１次微分算出手段と、
    　前記小領域ごとに、前記１次微分値の最大値と最小値との差分の絶対値を算出する位相変化量算出手段と、を備え、
    　前記絶対値が最大の小領域を、前記位相変化量が最も大きい小領域とすること
    　を特徴とする磁気共鳴撮影装置。
  6. 　請求項３記載の磁気共鳴撮影装置であって、
    　前記位相とび発生領域抽出手段は、
    　前記位相データの１次微分値を算出する変化量算出手段と、
    　前記小領域毎に、前記一次微分値をモデル関数でフィッティングするフィッティング手段と、
    　前記小領域毎に、前記フィッティング後のモデル関数の高さを算出する高さ算出手段と、を備え、
    　前記高さが最も高い小領域を、前記位相変化量が最も大きい小領域とすること
    　を特徴とする磁気共鳴撮影装置。
  7. 　請求項５記載の磁気共鳴撮影装置であって、
    　前記位相とび補正手段は、
    　前記位相とび発生領域として抽出された小領域における、前記１次微分値のピーク値の半値幅を、前記位相とびを補正する位相とび補正領域と設定する位相とび補正領域設定手段を備えること
    　を特徴とする磁気共鳴撮影装置。
  8. 　請求項６記載の磁気共鳴撮影装置であって、
    　前記位相とび補正手段は、
    　前記フィッティング後のモデル関数の半値幅の２倍の領域を、前記位相とびを補正する位相とび補正領域と設定する位相とび補正領域設定手段を備えること
    　を特徴とする磁気共鳴撮影装置。
  9. 　請求項７記載の磁気共鳴撮影装置であって、
    　前記位相とび補正手段は、前記位相とび補正領域の前記１次微分値を補間により接続し、前記位相とびを補正すること
    　を特徴とする磁気共鳴撮影装置。
  10. 　請求項９記載の磁気共鳴撮影装置であって、
    　前記補間は、前記１次微分値の前記位相とび補正領域の両端を直線で接続する直線補間であること
    　を特徴とする磁気共鳴撮影装置。
  11. 　請求項９記載の磁気共鳴撮影装置であって、
    　前記位相とび補正手段は、前記補間により接続後の前記１次微分値を積分し、前記位相とび補正後の位相データとすること
    　を特徴とする磁気共鳴撮影装置。
  12. 　請求項９記載の磁気共鳴撮影装置であって、
    　前記位相とび補正手段は、前記補間により接続後の前記１次微分値の累乗和を、前記位相とび補正後の位相データとすること
    　を特徴とする磁気共鳴撮影装置。
  13. 　請求項３記載の磁気共鳴撮影装置であって、
    　前記評価値算出手段は、前記補正後の位相データから算出した補正スペクトルのピークの高さと、補正前の位相データから算出した補正スペクトルのピークの高さとの差分を前記評価値とすること
    　を特徴とする磁気共鳴撮影装置。
  14. 　請求項３記載の磁気共鳴撮影装置であって、
    　前記評価値算出手段は、前記補正後の位相データから算出した補正スペクトルの標準偏差と、補正前の位相データから算出した補正スペクトルの標準偏差との差分を前記評価値とすること、
    　を特徴とする磁気共鳴撮影装置。
  15. 　請求項３記載の磁気共鳴撮影装置であって、
    　前記評価値算出手段は、前記補正後の位相データから算出した信号値であって、予め定めた値以上の値を有する信号値の数と、補正前の位相データから算出した信号値であって、予め定めた値以上の値を有する信号値の数との差を、前記評価値とすること
    　を特徴とする磁気共鳴撮影装置。
PCT/JP2010/058564 2009-05-27 2010-05-20 磁気共鳴撮影装置 WO2010137516A1 (ja)

Priority Applications (1)

Application Number Priority Date Filing Date Title
JP2011515997A JP5165791B2 (ja) 2009-05-27 2010-05-20 磁気共鳴撮影装置

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
JP2009127221 2009-05-27
JP2009-127221 2009-05-27

Publications (1)

Publication Number Publication Date
WO2010137516A1 true WO2010137516A1 (ja) 2010-12-02

Family

ID=43222626

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2010/058564 WO2010137516A1 (ja) 2009-05-27 2010-05-20 磁気共鳴撮影装置

Country Status (2)

Country Link
JP (1) JP5165791B2 (ja)
WO (1) WO2010137516A1 (ja)

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2013027539A1 (ja) * 2011-08-24 2013-02-28 株式会社日立メディコ 磁気共鳴撮影装置、位相値補正方法およびプログラム
CN106651981A (zh) * 2016-12-07 2017-05-10 深圳先进技术研究院 一种环形伪影修正的方法及装置
CN108279393A (zh) * 2017-01-05 2018-07-13 上海康达卡勒幅医疗科技有限公司 一种磁共振成像全自动预加重校正涡流的方法
JP2019048152A (ja) * 2011-07-07 2019-03-28 ザ ボード オブ トラスティーズ オブ ザ レランド スタンフォード ジュニア ユニバーシティー ボリューム測定位相コントラストｍｒｉによる総合的心血管解析

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPH05269094A (ja) * 1992-03-23 1993-10-19 Toshiba Corp 磁気共鳴診断装置
JP2002224083A (ja) * 2001-02-06 2002-08-13 Hitachi Medical Corp 磁気共鳴画像診断装置
WO2006040866A1 (ja) * 2004-10-13 2006-04-20 Hitachi Medical Corporation 磁気共鳴撮影装置及び磁気共鳴撮影方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPH05269094A (ja) * 1992-03-23 1993-10-19 Toshiba Corp 磁気共鳴診断装置
JP2002224083A (ja) * 2001-02-06 2002-08-13 Hitachi Medical Corp 磁気共鳴画像診断装置
WO2006040866A1 (ja) * 2004-10-13 2006-04-20 Hitachi Medical Corporation 磁気共鳴撮影装置及び磁気共鳴撮影方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
L.N.RYNER: "Effect of Oscillatory Eddy Currents on MR Spectroscopy", PROCEEDINGS OF THE INTERNATIONAL SOCIETY FOR MAGNETIC RESONANCE IN MEDICINE, 10 March 1998 (1998-03-10) *

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2019048152A (ja) * 2011-07-07 2019-03-28 ザ ボード オブ トラスティーズ オブ ザ レランド スタンフォード ジュニア ユニバーシティー ボリューム測定位相コントラストｍｒｉによる総合的心血管解析
WO2013027539A1 (ja) * 2011-08-24 2013-02-28 株式会社日立メディコ 磁気共鳴撮影装置、位相値補正方法およびプログラム
US9678187B2 (en) 2011-08-24 2017-06-13 Hitachi, Ltd. Magnetic resonance imaging device, phase value correction method and program
CN106651981A (zh) * 2016-12-07 2017-05-10 深圳先进技术研究院 一种环形伪影修正的方法及装置
CN108279393A (zh) * 2017-01-05 2018-07-13 上海康达卡勒幅医疗科技有限公司 一种磁共振成像全自动预加重校正涡流的方法

Also Published As

Publication number Publication date
JPWO2010137516A1 (ja) 2012-11-15
JP5165791B2 (ja) 2013-03-21

Similar Documents

Publication Publication Date Title
JP5398149B2 (ja) 磁気共鳴イメージング装置
US9317917B2 (en) Method, reconstruction device, and magnetic resonance apparatus for reconstructing magnetic resonance raw data
JP6014770B2 (ja) 磁気共鳴撮影装置および温度情報計測方法
US10481232B2 (en) Magnetic resonance imaging apparatus and magnetic resonance imaging method
US9977105B2 (en) Magnetic resonance imaging device and water-fat separation method
JP4889791B2 (ja) 磁気共鳴撮影装置
JP2005270285A (ja) 磁気共鳴イメージング装置および磁気共鳴イメージングデータ処理方法
US9696400B2 (en) Method and apparatus for optimization of a pulse sequence for a magnetic resonance system
JP4564015B2 (ja) 磁気共鳴撮影装置及び磁気共鳴撮影方法
JP5829687B2 (ja) 磁気共鳴撮影装置
US9746537B2 (en) Magnetic resonance imaging apparatus and magnetic resonance imaging method
JP5165791B2 (ja) 磁気共鳴撮影装置
Meineke et al. Data consistency‐driven determination of‐fluctuations in gradient‐echo MRI
US9851424B2 (en) Magnetic resonance imaging apparatus
US20140253121A1 (en) Magnetic resonance imaging apparatus and irradiation magnetic field distribution measurement method
JP5004967B2 (ja) 核磁気共鳴撮影装置，高周波磁場パルスの強度調整方法及び核磁気共鳴撮影方法
US10684340B1 (en) Systems and methods for predicting errors and optimizing protocols in quantitative magnetic resonance imaging
US10302720B2 (en) Method and apparatus for determining dephasing factors in magnetic resonance imaging and spectroscopy
JP5718148B2 (ja) 磁気共鳴イメージング装置及びＤｕａｌＳｌｉｃｅ計測方法
JP5127291B2 (ja) 核磁気共鳴撮影装置
US10203388B2 (en) Magnetic resonance imaging apparatus and magnetic resonance imaging method
US20180011160A1 (en) Reduction of eddy currents during flow encoded magnetic resonance imaging
JP6776217B2 (ja) 画像処理装置、画像処理方法、画像処理プログラム及び磁気共鳴イメージング装置
JP3945918B2 (ja) 磁気共鳴装置
JP4827545B2 (ja) 磁気共鳴撮影装置および方法

Legal Events

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

Ref document number: 10780471

Country of ref document: EP

Kind code of ref document: A1

WWE Wipo information: entry into national phase

Ref document number: 2011515997

Country of ref document: JP

NENP Non-entry into the national phase

Ref country code: DE

122 Ep: pct application non-entry in european phase

Ref document number: 10780471

Country of ref document: EP

Kind code of ref document: A1