WO2014185323A1 - 磁気共鳴イメージング装置および磁気共鳴イメージング方法 - Google Patents

磁気共鳴イメージング装置および磁気共鳴イメージング方法 Download PDF

Info

Publication number
WO2014185323A1
WO2014185323A1 PCT/JP2014/062306 JP2014062306W WO2014185323A1 WO 2014185323 A1 WO2014185323 A1 WO 2014185323A1 JP 2014062306 W JP2014062306 W JP 2014062306W WO 2014185323 A1 WO2014185323 A1 WO 2014185323A1
Authority
WO
WIPO (PCT)
Prior art keywords
phase
blade
phase correction
data
correction amount
Prior art date
Application number
PCT/JP2014/062306
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 JP2015517044A priority Critical patent/JP6410715B2/ja
Priority to CN201480022131.8A priority patent/CN105120745B/zh
Priority to US14/785,454 priority patent/US10048343B2/en
Publication of WO2014185323A1 publication Critical patent/WO2014185323A1/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/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
    • 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/4818MR characterised by data acquisition along a specific k-space trajectory or by the temporal order of k-space coverage, e.g. centric or segmented coverage of k-space
    • G01R33/4824MR characterised by data acquisition along a specific k-space trajectory or by the temporal order of k-space coverage, e.g. centric or segmented coverage of k-space using a non-Cartesian trajectory
    • 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
    • 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/56509Correction of image distortions, e.g. due to magnetic field inhomogeneities due to motion, displacement or flow, e.g. gradient moment nulling
    • 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/56563Correction of image distortions, e.g. due to magnetic field inhomogeneities caused by a distortion of the main magnetic field B0, e.g. temporal variation of the magnitude or spatial inhomogeneity of B0
    • 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/56572Correction of image distortions, e.g. due to magnetic field inhomogeneities caused by a distortion of a gradient magnetic field, e.g. non-linearity of a gradient magnetic field

Definitions

  • the present invention relates to a magnetic resonance imaging (hereinafter referred to as “MRI”) technique, and more particularly to an MRI technique using a non-orthogonal sampling method.
  • MRI magnetic resonance imaging
  • the MRI device measures NMR signals generated by the spins of the subject, especially the tissues of the human body, and visualizes the form and function of the head, abdomen, limbs, etc. in two or three dimensions Device.
  • different phase encoding is given by a gradient magnetic field, and a frequency-encoded NMR signal is measured as time-series data.
  • the measured NMR signal is reconstructed into an image by two-dimensional or three-dimensional Fourier transform.
  • Non-orthogonal sampling methods such as the radial sampling method and the hybrid radial method perform sampling by scanning the measurement space radially at various rotation angles with the rotation point as the center of the measurement space (generally the origin). The data necessary for image reconstruction for one sheet is obtained. Since sampling is performed in a radial manner, it is known to be strong against artifacts due to body movement. However, since scanning trajectories (blades) overlap in the measurement space, the image quality of the reconstructed image deteriorates if the positional relationship between the blades is inappropriate or if there is a phase difference at the intersection between the blades. .
  • one linear locus in the radial sampling method, a plurality of parallel linear locus in the hybrid radial method, and these are collectively referred to as a blade.
  • Patent Document 1 does not correct the phase difference at the intersection between the blades. If there is a phase difference at the location where the blades overlap in the measurement space, the signals cancel each other out, resulting in uneven pixel values in the image and a decrease in image formation.
  • the present invention has been made in view of the above circumstances, and an object of the present invention is to provide a technique for reducing deterioration in image quality due to a phase difference of intersections between scanning trajectories (blades) in measurement by a non-orthogonal sampling method. .
  • correction is performed to reduce (reduce) the phase difference of intersections between a plurality of scanning trajectories (blades) measured by a non-orthogonal sampling method.
  • the reduction of the phase difference is made by matching the phases of the intersections between the blades, and by matching the phase of each blade at the position including the shift amount in the frequency direction between all the blades. This is done by canceling the amount.
  • the present invention it is possible to reduce image quality deterioration due to a phase difference between scanning trajectories (blades) in measurement using a non-orthogonal sampling method.
  • Functional block diagram of the control processing system of the first embodiment (a) is an explanatory diagram for explaining the peak shift of the echo signal of the first embodiment, (b) is an enlarged view of the vicinity of the k-space origin of (a). (a) is an explanatory diagram for explaining the phase correction amount calculation method of the first embodiment, (b) is an enlarged view of the vicinity of the k-space origin of (a).
  • FIG. 1 is a block diagram showing the overall configuration of one embodiment of the MRI apparatus of this embodiment.
  • the MRI apparatus 100 of the present embodiment obtains a tomographic image of a subject using an NMR phenomenon.
  • the static magnetic field generator 120 generates a uniform static magnetic field in the direction perpendicular to the body axis in the space around the subject 101 if the vertical magnetic field method is used, and in the body axis direction if the horizontal magnetic field method is used.
  • the apparatus includes a permanent magnet type, normal conducting type or superconducting type static magnetic field generating source disposed around the subject 101.
  • the gradient magnetic field generation unit 130 includes a gradient magnetic field coil 131 wound in the three-axis directions of X, Y, and Z, which is a coordinate system (device coordinate system) of the MRI apparatus 100, and a gradient magnetic field power source that drives each gradient magnetic field coil 132, and by applying the gradient magnetic field power supply 132 of each gradient coil 131 in accordance with a command from the sequencer 140 described later, gradient magnetic fields Gx, Gy, and Gz are applied in the X, Y, and Z axis directions. .
  • the transmitter 150 irradiates the subject 101 with a high-frequency magnetic field pulse (hereinafter referred to as “RF pulse”) in order to cause nuclear magnetic resonance to occur in the nuclear spins of the atoms constituting the biological tissue of the subject 101.
  • RF pulse high-frequency magnetic field pulse
  • the high frequency oscillator 152 generates an RF pulse.
  • the modulator 153 amplitude-modulates the output RF pulse in accordance with a command from the sequencer 140.
  • the high-frequency amplifier 154 amplifies the amplitude-modulated RF pulse and supplies the amplified RF pulse to the transmission coil 151 disposed close to the subject 101.
  • the transmission coil 151 irradiates the subject 101 with the supplied RF pulse.
  • the receiving unit 160 detects a nuclear magnetic resonance signal (echo signal, NMR signal) emitted by nuclear magnetic resonance of the nuclear spin constituting the biological tissue of the subject 101, and receives a high-frequency coil (receiving coil) on the receiving side. 161, a signal amplifier 162, a quadrature detector 163, and an A / D converter 164.
  • the reception coil 161 is disposed in the vicinity of the subject 101 and detects an NMR signal in response to the subject 101 induced by the electromagnetic wave irradiated from the transmission coil 151.
  • the detected NMR signal is amplified by the signal amplifier 162 and then divided into two orthogonal signals by the quadrature phase detector 163 at the timing according to the command from the sequencer 140, and each is digitally converted by the A / D converter 164. The amount is converted and sent to the control processing unit 170.
  • the sequencer 140 applies an RF pulse and a gradient magnetic field pulse in accordance with an instruction from the control processing unit 170. Specifically, in accordance with instructions from the control processing unit 170, various commands necessary for collecting tomographic image data of the subject 101 are transmitted to the transmission unit 150, the gradient magnetic field generation unit 130, and the reception unit 160.
  • the control processing unit 170 controls the entire MRI apparatus 100, performs various data processing operations, displays and stores processing results, and includes a CPU 171, a storage device 172, a display device 173, and an input device 174.
  • the storage device 172 includes an internal storage device such as a hard disk and an external storage device such as an external hard disk, an optical disk, and a magnetic disk.
  • the display device 173 is a display device such as a CRT or a liquid crystal.
  • the input device 174 is an interface for inputting various control information of the MRI apparatus 100 and control information of processing performed by the control processing unit 170, and includes, for example, a trackball or a mouse and a keyboard.
  • the input device 174 is disposed in the vicinity of the display device 173. The operator interactively inputs instructions and data necessary for various processes of the MRI apparatus 100 through the input device 174 while looking at the display device 173.
  • the CPU 171 realizes each process of the control processing unit 170 such as control of the operation of the MRI apparatus 100 and various data processes by executing a program stored in the storage device 172 in advance according to an instruction input by the operator.
  • the command to the sequencer 140 is performed according to a pulse sequence stored in the storage device 172 in advance.
  • the CPU 171 executes signal processing, image reconstruction processing, and the like, and displays the tomographic image of the subject 101 as a result thereof on the display device 173. And stored in the storage device 172.
  • the transmission coil 151 and the gradient magnetic field coil 131 are opposed to the subject 101 in the static magnetic field space of the static magnetic field generation unit 120 into which the subject 101 is inserted if the vertical magnetic field method is used, and if the horizontal magnetic field method is used. It is installed so as to surround the subject 101. Further, the receiving coil 161 is installed so as to face or surround the subject 101.
  • the nuclide to be imaged by the MRI apparatus which is widely used clinically, is a hydrogen nucleus (proton) which is a main constituent material of the subject 101.
  • the MRI apparatus 100 by imaging information on the spatial distribution of proton density and the spatial distribution of relaxation time in the excited state, the form or function of the human head, abdomen, limbs, etc. can be expressed two-dimensionally or three-dimensionally. Take an image.
  • the control processing unit 170 of the present embodiment measures echo signals along a plurality of scanning trajectories in a predetermined measurement space (k space) according to a pulse sequence of a non-orthogonal sampling method, A measurement unit 210 that is arranged on the scanning locus as a data string, and an image reconstruction unit 220 that reconstructs an image from a plurality of echo signals and obtains a reconstructed image.
  • the measurement unit 210 controls the operations of the static magnetic field generation unit 120, the gradient magnetic field generation unit 130, the high frequency magnetic field generation unit 150, and the high frequency magnetic field detection unit 160 according to the pulse sequence of the non-orthogonal sampling method, and samples the echo signal.
  • the obtained data string is arranged on a scanning locus in k space.
  • each scanning locus is hereinafter referred to as a blade.
  • the image reconstruction unit 220 reduces (reduces) the phase difference in the low spatial frequency region between the data strings on each blade in the k space, and obtains a reconstructed image.
  • the image reconstruction unit 220 includes a phase correction amount calculation unit 221 that calculates a phase correction amount (phase correction amount) for each data sequence on a plurality of blades, and a calculated phase correction.
  • a phase correction unit 222 that corrects the phase of each of the data sequences on the plurality of blades using the amount, and a reconstruction unit 223 that calculates a reconstructed image from the corrected data sequence on the plurality of blades.
  • the phase correction amount calculated by the phase correction amount calculation unit 221 is calculated so as to reduce the phase difference between the data strings on the plurality of blades.
  • the phase correction amount calculation unit 221 calculates the phase correction amount so that the phases of the data at predetermined positions on the blade are aligned in each of the plurality of data strings.
  • the predetermined position is determined by reflecting a peak shift amount that is a shift amount from the k-space origin of the center (peak position) of the echo signal.
  • the predetermined position is an intersection of a plurality of blades.
  • the phase correction amount calculation unit 221 uses the data on the blades other than the reference blade as the intersection of a predetermined reference blade (reference blade) among a plurality of blades and a blade other than the reference blade. The phase correction amount is calculated so that the phase at the intersection of the rows matches the phase at the intersection of the data rows on the reference blade.
  • the phase correction amount calculation unit 221 of the present embodiment calculates the position of the intersection using the peak shift amount for each data string on the blades other than the reference blade, and uses the calculated information on the position of the intersection from the intersection. Calculate the distance to the midpoint of the data string and the distance to the midpoint of the data string on the reference blade (reference data string), and use the calculated distances to determine the phase at the intersection of the data string and the reference data. Calculate the phase at the intersection of the columns.
  • the peak shift amount of the echo signal is used when calculating the phase correction amount.
  • the peak shift which is the shift of the peak position of the echo signal, is caused by the nonuniformity of the static magnetic field and the output error of the gradient magnetic field. Further, the peak position of the echo signal changes according to the area of the dephase pulse applied before the readout gradient magnetic field pulse is applied. This is because the timing at which the phases of the echo signals are aligned by the application of the readout gradient magnetic field pulse predetermined by the pulse sequence after dephasing is the peak position.
  • the deviation amount (peak shift amount) of the echo signal is calculated separately by performing pre-scanning in advance. Thereby, it can prevent that imaging time is prolonged. For example, as described in Japanese Patent Application Laid-Open No. 2005-152175, the calculation is performed by acquiring an echo signal by a pulse sequence for acquiring only a specific echo signal and using the echo signal.
  • FIG. 3 (a) is a diagram for explaining the shift of the blade 301, which is the scanning locus, due to the peak shift of the echo signal in the k space.
  • the white circle position 305 is the position of the data constituting each blade 301 when there is no shift (displacement) in the blade 301
  • the black circle position 304 is data when there is a shift (actual). Position.
  • the blade center 301c of the blade 301 arranged on the kx axis is the position 303 of the origin of the k space. However, since there is a shift, the blade center 301c is actually the position 301c.
  • the shift amount of the blade 301 the value (coordinate value) of the coordinate Pdn on the k-space coordinates of the shifted blade center 301c is calculated.
  • the shift amount in the X-axis direction is ⁇ D'x
  • the shift amount in the Y-axis direction is ⁇ D'y due to the non-uniformity of the static magnetic field, the peak shift caused by the output error of the gradient magnetic field, etc.
  • Fig. 3 (b) As shown, the shift amount ⁇ D'x_d in the kx direction and the shift amount ⁇ D'y_d in the ky direction of the d-th blade (d is an integer equal to or greater than 1) are expressed by the following equations (1) and (2). expressed.
  • ⁇ d is an angle formed by the d-th blade 301 with the X axis.
  • FIG. 4A and FIG. 4B are diagrams for explaining a method of calculating the phase correction amount according to the present embodiment.
  • FIG. 4 (b) is an enlarged view of the vicinity of the center of the k space in FIG. 4 (a).
  • the first blade 311 is arranged on the kx axis in the k space, and this is used as the reference blade 311.
  • the reference blade 311 is not limited to the blade on the kx axis, and may be any blade.
  • intersection point of the b-th blade (where b is an integer of 2 or more) 312 and the reference blade 311 is defined as an intersection point 313.
  • a position 315 indicated by a white circle is an ideal data position when there is no positional shift (shift) in the blade
  • a position 314 indicated by a black circle is actual data reflecting the shift of the blade. Is the position.
  • the b-th blade 312 intersects the reference blade 311 at the origin 316 of the k space.
  • the two intersect at a position (intersection) 313 indicated by a white cross (cross).
  • phase correction is performed so that the phase at the intersection point 313 of the b-th blade 312 is aligned with the phase at the intersection point 313 of the reference blade.
  • the phase correction amount calculation unit 221 of the present embodiment obtains the coordinates of the intersection point 313, obtains the distances ⁇ Db1 and ⁇ Dbb between the intersection point 313 and the midpoints (blade centers) 311c and 312c of the blades 311 and 312, and each blade
  • the phase values (Phase_1, Phase_b) of 311 and 312 at the intersection 313 are obtained.
  • the difference (phase difference) between the phase value of the b-th blade 312 at the intersection point 313 and the phase value at the intersection point 313 of the reference blade 311 is calculated as the phase correction amount PhC_b.
  • the distances ⁇ Db1 and ⁇ Dbb between the intersection point 313 and the respective blade centers 311c and 312c are obtained from the peak shift amounts ⁇ D′x and ⁇ D′y in the X-axis direction and the Y-axis direction and the blade angle ⁇ b of the b-th blade 312.
  • the blade angle ⁇ b is an angle formed by the b-th blade 312 and the X axis.
  • the coordinates PIb ⁇ x, y ⁇ of the intersection point 313 can be expressed by the following equation (3) when expressed using the peak shift amounts ⁇ D′x and ⁇ D′y and the blade angle ⁇ b: 4) and formula (5).
  • ⁇ Db1 ⁇ D'x ⁇ 1-cos ( ⁇ b) ⁇ + ⁇ D'ycos ( ⁇ b) (7)
  • ⁇ Dbb ⁇ D'y (8) ⁇
  • Blade_1 () and Blade_b () are data strings of the reference blade 311 and the b-th blade 312
  • Real () is real part data
  • Imaginary () is imaginary part data
  • N is each blade.
  • CENTER represents the positions of the middle points 311c and 312c, respectively.
  • phase correction amount (phase difference) PhC_b is obtained by the following equation (12).
  • PhC_b Phase_b-Phase_1 (12) By applying this phase correction amount PhC_b to the data string of the b-th blade 312 and correcting the phase value, the phase values of the first blade 311 and the b-th blade 312 at the intersection 313 are matched.
  • the number of b-th blade 312 is 1 or more. If there is a peak shift (shift at the center of the blade), the intersection of the reference blade 311 and each b-th blade 312 does not concentrate on one point. Therefore, the phase correction amount calculation unit 221 calculates the position of the intersection point 313 with the reference blade indicated by the white cross for each b-th blade 312, the distance ⁇ Db1 from the midpoint 311c of the reference blade 311, and the blade The distance ⁇ Dbb from the midpoint 312c of 312 is calculated by the above procedure, the phase values Phase_1 and Phase_b of the midpoints 311c and 312c are obtained, and the phase correction amount PhC_b is calculated.
  • phase correction unit 222 corrects the phases of all data constituting the b-th blade 312 using the calculated phase correction amount PhC_b.
  • phase correction is performed on all data constituting the blade 312 using the same phase correction amount. That is, the corrected complex data sequence Blade_b (x) of the b-th blade 312 is expressed by the following equation (13).
  • Blade_b (x)
  • represents the amplitude of the complex data string of the b-th blade. I is an imaginary unit.
  • the reconstruction unit 223 reconstructs an image using the corrected data sequence of each blade.
  • FIG. 6 is a processing flow of the phase correction processing of this embodiment.
  • the first acquired blade (first blade) is used as a reference blade.
  • b is used as a counter for counting blades.
  • the total number of blades is M (M is an integer of 1 or more).
  • the phase correction amount calculation unit 221 determines whether or not the phase correction processing has been completed for all the blades (step S1104). Here, it is determined whether or not the counter b has exceeded the blade number M. If it is determined that all the blades 312 have been processed, the processing ends.
  • the phase correction amount calculation unit 221 first calculates the coordinates PIb ⁇ x of the intersection point 313 with the reference blade 311 for the b-th blade 312 acquired by the measurement unit 210. , Y ⁇ is calculated (step S1105).
  • phase correction amount calculation unit 221 calculates the phase Phase_1 of the intersection point 313 of the reference blade 311 and the phase Phase_b of the intersection point 313 of the b-th blade 312 (step S1106). Then, the phase correction amount PhC_b of the b-th blade is calculated as both phase differences (step S1107).
  • phase correction unit 222 corrects the phase of the b-th blade according to the above equation (13) using the calculated phase correction amount PhC_b (step S1108).
  • the phase correction amount calculation unit 221 returns to step S1103 and repeats the process.
  • the phase correction amount is determined so that the phases of the intersections of the blades are aligned for the data sequence of the plurality of blades acquired according to the non-orthogonal pulse sequence, and the phase correction is performed.
  • a reference blade is defined among a plurality of blades, and the phase correction amount is determined so that the phase of each of the other blades coincides with the reference blade, and each blade is configured with the phase correction amount. Correct the phase of the data. As a result, the phase difference is reduced at the point where the blades overlap with each other, and signal cancellation due to this is also reduced. At this time, the phase is corrected with the amount of phase change taking into account the deviation of the center of the reception echo of the blade.
  • each blade intersects the reference blade in a low spatial frequency region. Therefore, the phase difference in the low spatial frequency region that most affects the image quality is reduced by matching the phase of the intersection with the reference blade for each blade to the phase of the reference blade. Thereby, image quality degradation can be reduced.
  • the phase at the intersection point 313 is matched, the interpolation value of only one intersection point is used, but the present invention is not limited to this.
  • the phase may be matched using the average value of the intersection 313 and the surrounding points.
  • the phase correction is performed in the k space, but the space in which the phase correction is performed may be an image space.
  • the following equation (14) is followed.
  • FT [Blade_b (x)]
  • FT [] represents a Fourier transform.
  • the data to be held may be all data of the reference blade, or may be data for several points around the origin of the reference blade in order to reduce the memory used.
  • Second Embodiment a second embodiment of the present invention will be described.
  • the reference blade is determined, and the phase of other blades is corrected so as to match the phase of the reference blade at the intersection with the reference blade.
  • correction is performed so that the phases at the origin offset position POb of each blade described later match. Therefore, in this embodiment, setting of the reference blade is unnecessary.
  • the MRI apparatus of the present embodiment has basically the same configuration as the MRI apparatus 100 of the first embodiment. However, the calculation method of the phase correction amount is different as described above. Therefore, the processing of the phase correction amount calculation unit 221 is different. Hereinafter, the present embodiment will be described focusing on the configuration different from the first embodiment. Also in the present embodiment, the measurement unit 210 measures echo signals along a plurality of scanning trajectories in a predetermined measurement space (k space) in accordance with a pulse sequence of a non-orthogonal sampling method, as a data string They are arranged on the scanning locus.
  • k space predetermined measurement space
  • the phase correction amount is calculated so that the phases of data at predetermined positions on the blade of each of the plurality of data strings are aligned.
  • the predetermined position is an origin offset position that is an intersection of perpendicular lines drawn from the origin of the k space to each of the plurality of blades, and the phase correction amount calculation unit 221 has all phases of the origin offset positions of the plurality of data strings.
  • the phase correction amount is calculated as follows.
  • the phase correction amount calculation unit 221 of the present embodiment calculates an offset position using a peak shift amount for each of a plurality of data strings, and uses the calculated offset position information to calculate the data string from the offset position. The distance to the midpoint is obtained, and the phase at the offset position of the data string is calculated using the distance.
  • the shift amount (peak shift amount) of the echo signal is separately calculated by performing pre-scanning in advance.
  • phase correction amount calculation processing by the phase correction amount calculation unit 221 of the present embodiment will be described with reference to the drawings.
  • FIG. 7 (a) shows the positional relationship between the b-th blade 322 of this embodiment and the kx axis and ky axis.
  • FIG. 7 (b) is an enlarged view of the vicinity of the center of the k space in FIG. 7 (a).
  • b is an integer of 1 or more.
  • a position 325 indicated by a white circle is an ideal data position when there is no positional shift (shift) in the blade 322, and a position 324 indicated by a black circle is actual data reflecting the shift of the blade 322. Is the position.
  • each blade 322 is located at the origin 326 of the k space indicated by a white circle.
  • it shifts due to the output response of the gradient magnetic field, and becomes a position 322c indicated by a black circle.
  • the intersection point 323 when a perpendicular is drawn from the origin of the k space to the b-th blade 322 is referred to as an origin offset position. In the figure, it is indicated by white cross.
  • the phase of each blade 322 is corrected so that the phases of the origin offset positions of the blades 322 are aligned.
  • the phase of the origin offset position 323 of each blade 322 is calculated by obtaining a distance (error) ⁇ Db in the k space between the blade center 322c of each blade 322 and the origin offset position 323.
  • the distance ⁇ Db in the k space from the origin offset position 323 is referred to as a shift amount in the frequency direction.
  • ⁇ D′x and ⁇ D′y are blade shift amounts in the X-axis direction and the Y-axis direction, as in the first embodiment.
  • ⁇ b is an angle formed by the b-th blade 322 and the X axis.
  • the coordinate value ⁇ x, y ⁇ of the coordinate POb of the origin offset 323 has the relationship of the following expressions (16) and (17).
  • the shift amount ⁇ Db in the frequency direction of the b-th blade 322 is calculated using the following equation (19) from the equation (20) using the coordinate position POb of the origin offset position 323 and the coordinate Pbn of the blade center 322c of the b-th blade 322. ).
  • ⁇ Db ⁇ D'xcos 2 ( ⁇ b) + ⁇ D'ysin 2 ( ⁇ b) (19) ⁇
  • Phase_b is calculated by the following equation (21).
  • Blade_b is the data string of the b-th blade 322
  • Real () is the real part data
  • Imaginary () is the imaginary part data
  • N is the number of data points of each blade
  • CENTER is the middle point.
  • Each position is represented.
  • phase correction amount PhC_b is expressed by the following equation (22).
  • PhC_b Phase_b- ⁇ (22) For example, 0 is used for ⁇ .
  • the phase correction amount calculation unit 221 obtains the phase value Phase_b of the origin offset position 323 for each b-th blade and calculates the position correction amount.
  • phase correction unit 222 corrects the phases of all the data constituting the b-th blade 322 using the calculated phase correction amount PhC_b.
  • the corrected complex data sequence Blade_b (x) of the b-th blade 322 is expressed by the following equation (23), as in the first embodiment.
  • Blade_b (x)
  • represents the amplitude of the complex data string of the b-th blade. I is an imaginary unit.
  • FIG. 9 is a processing flow of the phase correction processing of this embodiment.
  • b is used as a counter for counting blades.
  • the total number of blades is M (M is an integer of 1 or more).
  • the phase correction amount calculation unit 221 determines whether or not the phase correction processing has been completed for all the blades (step S1202). Here, it is determined whether or not the counter b has exceeded the blade number M. If it is determined that all the blades 312 have been processed, the processing ends.
  • the phase correction amount calculation unit 221 calculates the shift amount ⁇ Db in the frequency direction for the b-th blade 322 acquired by the measurement unit 210 (step S1203).
  • the phase correction amount calculation unit 221 calculates the phase of the origin offset position 323 using the calculated shift amount in the frequency direction (step S1204), and calculates the phase correction amount PhC_b of the b-th blade 322 (step S1205).
  • the phase correction unit 222 corrects the b-th blade according to the above equation (23) using the calculated phase correction amount PhC_b (step S1206).
  • the phase correction amount calculation unit 221 increments the counter b by 1 (step S1207), returns to step S1202, and repeats the process.
  • the data sequence of each blade is corrected so that the phase at the origin offset position of each blade is aligned for the data sequence of a plurality of blades acquired according to the non-orthogonal pulse sequence. To do. Thereby, in the low spatial frequency region, the phase shift between the blades is eliminated, and cancellation of signals generated by the phase shift can be reduced. Therefore, image quality deterioration can be suppressed.
  • the phase may be adjusted using an average value of positions around the position 324 as in the first embodiment.
  • the phase correction may be performed in the image space as in the first embodiment.
  • the phase correction amount is determined by calculating the phase value at the center of each blade from the change in the reception frequency during off-center imaging.
  • the MRI apparatus of the present embodiment has basically the same configuration as the MRI apparatus 100 of the first embodiment. However, the calculation method of the phase correction amount is different as described above. Therefore, the processing of the phase correction amount calculation unit 221 is different. Hereinafter, the present embodiment will be described focusing on the configuration different from the first embodiment. Also in the present embodiment, the measurement unit 210 measures echoes respectively along a plurality of scanning trajectories in a predetermined k space according to the pulse sequence of the non-orthogonal sampling method, and as a data string on the scanning trajectory (blade). To place. At this time, off-center imaging is performed by changing the reception frequency.
  • the phase correction amount calculation processing by the phase correction amount calculation unit 221 of the present embodiment will be described.
  • the phase correction amount calculation unit 221 of the present embodiment calculates the phase change amount at the midpoint of each of the plurality of data strings using the off-center distance in the k space, and uses the calculated phase change amount as the phase correction amount. That is, the phase value of the blade center is calculated from the amount of phase change during echo signal acquisition (between A / D) during off-center imaging.
  • phase correction amount calculation unit 221 of the present embodiment uses the calculated phase rotation amount as the phase correction amount.
  • the specific calculation procedure is as follows.
  • the reception frequency Rf is obtained by the following equation (24).
  • Rf BW ⁇ OffcD / FOV ... (24)
  • BW is a reception frequency band [Hz]
  • OffcD is an off-center distance [m]
  • FOV is a visual field size [m].
  • phase rotation amount (phase change amount) in the blade is expressed by the following (25).
  • n is an element number of a data string constituting the blade.
  • phase change amount ⁇ x in the kx axis direction of the k space and the phase change amount ⁇ y in the ky axis direction at the center of the blade are expressed by the following equations (26) and (27), respectively.
  • OffcD_X and OffcD_Y are the off-center distances [m] in the X-axis direction and Y-axis direction, respectively
  • FOV_X and FOV_Y are the field sizes [m] in the X-axis direction and Y-axis direction, respectively
  • CENTER is the midpoint Position.
  • phase change amount ⁇ b at the b-th blade center is expressed by the following equation (28).
  • the phase correction amount calculation unit 221 of the present embodiment obtains the blade center phase change amount ⁇ b of each blade according to the reception frequency shift amount by the above procedure, and determines the phase correction amount PhC_b so that all of them are aligned. To do. For example, if the blade center phases are all ⁇ [rad], the phase correction amount PhC_b is expressed by the following equation (29).
  • PhC_b ⁇ b- ⁇ ... (29) For example, 0 is used for ⁇ .
  • phase of all data constituting the b-th blade is corrected using the calculated phase correction amount PhC_b.
  • the corrected complex data string Blade_b (x) of the b-th blade 322 is expressed by the following equation (30), as in the first embodiment.
  • Blade_b (x)
  • the phase rotation amount at the blade center of each blade is obtained by calculation, and the phase of each blade is corrected so as to reduce the phase difference between the blades. Thereby, the phase shift between the blades is reduced. Accordingly, it is possible to suppress image quality degradation caused by the phase shift between the blades.
  • the phase correction may be performed in the image space as in the first embodiment.
  • the peak shift amount of the echo signal is not used, but the peak shift due to the static magnetic field inhomogeneity and the area error of the gradient magnetic field waveform can be considered.
  • the phase change amounts ⁇ ′x and ⁇ ′y of each blade in this case are expressed by the following equations (31) and (32).
  • phase change amount ⁇ ′b at the blade center of the b-th blade is expressed by the following equation (33).
  • phase correction amount PhC_b is determined using the phase change amount ⁇ ′b.
  • the imaging cross section is the XY plane of the apparatus coordinate system
  • the imaging cross section may be an arbitrary cross section.
  • the peak shift amount in each axis of the apparatus coordinate system is developed in the measurement coordinate system according to the angle of each blade, and the distance (deviation amount) ⁇ Db1 from the intersection with the reference blade to the midpoint of each blade ⁇ Dbb is calculated.
  • the specific development and calculation methods are as follows.
  • the coordinate k R in the measurement coordinate system reflecting the peak shift amount of each of the XYZ axes is expressed by the following equation (34).
  • k RA is a coordinate in the apparatus coordinate system
  • G RR (bl) is a readout gradient magnetic field [T] of each axis of XYZ
  • d A is a peak shift amount of each axis of XYZ
  • b is a blade number
  • ⁇ b is a blade number.
  • n is the data point number in the blade
  • Nc is the data number of the middle point position
  • e is the unit matrix
  • is the magnetic rotation ratio [rad / T]
  • ⁇ t is the sampling interval [s] of the data string
  • G is the gradient magnetic field strength [T].
  • ROM is a rotation matrix for converting the measurement coordinate system to the apparatus coordinate system, and is defined by the following equation (35).
  • sx, sy, and sz are projection components of the gradient magnetic field in the slice axis direction in the measurement coordinate system onto the XYZ axes of the apparatus coordinate system
  • px, py, and pz are in the phase encode axis (ky) direction in the measurement coordinate system.
  • Projection components of the gradient magnetic field on the XYZ axis of the apparatus coordinate system fx, fy, and fz are projection components on the XYZ axis of the apparatus coordinate system of the gradient magnetic field in the readout gradient magnetic field axis (kx) direction in the measurement coordinate system.
  • K R (b, n) is the coordinate in the measurement coordinate system of the nth data point of the bth blade. The same shall apply hereinafter.
  • equations (10) to (13) are calculated to perform phase correction.
  • the embodiment of the present invention is not limited to the above-described embodiment.
  • the radial sampling method that scans the k-space radially among the non-orthogonal sampling methods is illustrated and described as an example, but the sampling method used is not limited thereto. Any sampling method that draws a locus in which each blade overlaps in the k space may be used.
  • a hybrid radial method in which phase encoding is combined with a radial sampling method may be used.
  • the phase correction amount PhC_b is calculated by the method of each of the above embodiments using one central locus among a plurality of parallel linear loci, and the phase correction amount PhC_b is used to calculate 1
  • the phase of data on all the parallel linear trajectories constituting the blade is corrected.
  • 100 MRI apparatus 101 subject, 120 static magnetic field generation unit, 130 gradient magnetic field generation unit, 131 gradient magnetic field coil, 132 gradient magnetic field power supply, 140 sequencer, 150 high frequency magnetic field generation unit, 150 transmission unit, 151 transmission coil, 152 high frequency oscillator , 153 modulator, 154 high frequency amplifier, 160 high frequency magnetic field detection unit, 160 reception unit, 161 reception coil, 162 signal amplifier, 163 quadrature phase detector, 164 A / D converter, 170 control processing unit, 171 CPU, 172 storage Device, 173 display device, 174 input device, 210 measurement unit, 220 image reconstruction unit, 221 phase correction amount calculation unit, 222 phase correction unit, 223 reconstruction unit, 301 blade, 301c blade center, 303 position, 304 position, 305 position, 311 reference blade, 311c midpoint, 312 blade, 312c midpoint, 313 intersection, 314 position, 315 position, 316 origin, 322 blade, 322c blade center, 323 origin offset Position, 324 position, 325 position, 3

Landscapes

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

Abstract

 非直交系サンプリング法による計測において、走査軌跡(ブレード)間の位相差による画質の劣化を低減する技術を提供することを目的とする。そのために、本発明は、画像再構成時に、非直交系サンプリング法により計測した複数の走査軌跡(ブレード)間の位相差を縮小する補正を行う。例えば、位相差の縮小は、ブレード間の交点の位相を合致させる、ブレードの、周波数方向におけるシフト量を加味した位置での位相を全ブレード間で合致させる、計算で求めた各ブレードの位相変化量をキャンセルするといった手法で行う。

Description

磁気共鳴イメージング装置および磁気共鳴イメージング方法
 本発明は、磁気共鳴イメージング(以下、「MRI」という)技術に関し、特に、非直交系サンプリング法を用いたMRI技術に関する。
 MRI装置は、被検体、特に人体の組織を構成する原子核スピンが発生するNMR信号を計測し、その頭部、腹部、四肢等の形態や機能を2次元的に或いは3次元的に画像化する装置である。撮像においては、予め定められたシーケンスに従って、傾斜磁場により異なる位相エンコードが付与され、周波数エンコードされたNMR信号が時系列データとして計測される。そして、計測されたNMR信号は、2次元又は3次元フーリエ変換されることにより画像に再構成される。
 ラディアルサンプリング法、ハイブリッドラディアル法などの、非直交系サンプリング法は、計測空間の略一点(一般的には原点)を回転中心として、様々な回転角で放射状に計測空間を走査してサンプリングを行い、一枚の画像再構成に必要なデータを得る。放射状にサンプリングを行うことから、体動によるアーチファクトには強いことが知られている。しかし、計測空間で走査軌跡(ブレード)同士が重なりあうことから、ブレード間の位置関係が不適切であったり、ブレード間の交点に位相差が生じていたりすると、再構成画像の画質が劣化する。
 なお、以下、本明細書では、ラディアルサンプリング法における1本の直線軌跡、ハイブリッドラディアル法における平行な複数の直線軌跡、これらをまとめてブレードと呼ぶ。
 ところが、実際の撮像では、静磁場の不均一や傾斜磁場の出力誤差により、各ブレードの計測空間内の配置位置が計算上の配置位置(座標)と異なったり、ブレード間の交点に位相差が発生したりする。ブレードの配置位置の誤りを補正する技術として、傾斜磁場の出力誤差によって生じるブレード位置のシフト量を算出するためのデータを取得し、各ブレードのk空間上におけるシフト量を算出し、画像再構成処理の中で補正する手法がある(例えば、特許文献1参照)。
国際公開第2008/152937号
 しかしながら、特許文献1に開示の技術では、ブレード間の交点の位相差は補正されない。計測空間のブレードとブレードとが重なり合う箇所において、位相差があると、信号が打ち消しあい、画像に画素値のムラや結像性の低下が発生する。
 本発明は、上記事情に鑑みてなされたもので、非直交系サンプリング法による計測において、走査軌跡(ブレード)間の交点の位相差による画質の劣化を低減する技術を提供することを目的とする。
 本発明は、画像再構成時に、非直交系サンプリング法により計測した複数の走査軌跡(ブレード)間の交点の位相差を縮小(低減)する補正を行う。例えば、位相差の縮小は、ブレード間の交点の位相を合致させる、周波数方向におけるシフト量を加味した位置での各ブレードの位相を全ブレード間で合致させる、計算で求めた各ブレードの位相変化量をキャンセルするといった手法で行う。
 本発明によれば、非直交系サンプリング法による計測において、走査軌跡(ブレード)間の位相差による画質の劣化を低減できる。
第一の実施形態の磁気共鳴イメージング装置の全体構成を示すブロック図 第一の実施形態の制御処理系の機能ブロック図 (a)は、第一の実施形態のエコー信号のピークシフトを説明するための説明図であり、(b)は、(a)のk空間原点付近の拡大図 (a)は、第一の実施形態の位相補正量算出手法を説明するための説明図であり、(b)は、(a)のk空間原点付近の拡大図 第一の実施形態の位相補正量算出手法を説明するための説明図 第一の実施形態の位相補正処理のフローチャート (a)は、第二の実施形態の位相補正量算出手法を説明するための説明図であり、(b)は、(a)のk空間原点付近の拡大図 第二の実施形態の位相補正量算出手法を説明するための説明図 第二の実施形態の位相補正処理のフローチャート 第三の実施形態のブレード内の位相変化量を説明するための説明図
 <<第一の実施形態>>
 以下、添付図面に従って本発明に係る好ましい実施形態について詳説する。なお、発明の実施形態を説明するための全図において、特に明示しない限り、同一機能を有するものは同一符号を付け、その繰り返しの説明は省略する。
 まず、本実施形態のMRI装置の一例の全体概要を図1に基づいて説明する。図1は、本実施形態のMRI装置の一実施形態の全体構成を示すブロック図である。
 本実施形態のMRI装置100は、NMR現象を利用して被検体の断層画像を得るもので、図1に示すように、静磁場発生部120と、傾斜磁場発生部130と、高周波磁場発生部(以下、送信部)150と、高周波磁場検出部(以下、受信部)160と、制御処理部170と、シーケンサ140と、を備える。
 静磁場発生部120は、垂直磁場方式であれば、被検体101の周りの空間にその体軸と直交する方向に、水平磁場方式であれば、体軸方向に、均一な静磁場を発生させるもので、被検体101の周りに配置される永久磁石方式、常電導方式あるいは超電導方式の静磁場発生源を備える。
 傾斜磁場発生部130は、MRI装置100の座標系(装置座標系)であるX、Y、Zの3軸方向に巻かれた傾斜磁場コイル131と、それぞれの傾斜磁場コイルを駆動する傾斜磁場電源132とを備え、後述のシーケンサ140からの命令に従ってそれぞれの傾斜磁場コイル131の傾斜磁場電源132を駆動することにより、X、Y、Zの3軸方向に傾斜磁場Gx、Gy、Gzを印加する。
 送信部150は、被検体101の生体組織を構成する原子の原子核スピンに核磁気共鳴を起こさせるために、被検体101に高周波磁場パルス(以下、「RFパルス」と呼ぶ。)を照射するもので、高周波発振器(シンセサイザ)152と変調器153と高周波増幅器154と送信側の高周波コイル(送信コイル)151とを備える。高周波発振器152はRFパルスを生成する。変調器153は、シーケンサ140からの指令に従って、出力されたRFパルスを振幅変調する。高周波増幅器154は、この振幅変調されたRFパルスを増幅し、被検体101に近接して配置された送信コイル151に供給する。送信コイル151は供給されたRFパルスを被検体101に照射する。
 受信部160は、被検体101の生体組織を構成する原子核スピンの核磁気共鳴により放出される核磁気共鳴信号(エコー信号、NMR信号)を検出するもので、受信側の高周波コイル(受信コイル)161と信号増幅器162と直交位相検波器163と、A/D変換器164とを備える。受信コイル161は、被検体101に近接して配置され、送信コイル151から照射された電磁波によって誘起された被検体101の応答のNMR信号を検出する。検出されたNMR信号は、信号増幅器162で増幅された後、シーケンサ140からの指令によるタイミングで直交位相検波器163により直交する二系統の信号に分割され、それぞれがA/D変換器164でディジタル量に変換されて、制御処理部170に送られる。
 シーケンサ140は、制御処理部170からの指示に従って、RFパルスと傾斜磁場パルスとを印加する。具体的には、制御処理部170からの指示に従って、被検体101の断層画像のデータ収集に必要な種々の命令を送信部150、傾斜磁場発生部130、および受信部160に送信する。
 制御処理部170は、MRI装置100全体の制御、各種データ処理等の演算、処理結果の表示及び保存等を行うもので、CPU171と記憶装置172と表示装置173と入力装置174とを備える。記憶装置172は、ハードディスクなどの内部記憶装置と、外付けハードディスク、光ディスク、磁気ディスクなどの外部記憶装置とにより構成される。表示装置173は、CRT、液晶などのディスプレイ装置である。入力装置174は、MRI装置100の各種制御情報や制御処理部170で行う処理の制御情報の入力のインタフェースであり、例えば、トラックボールまたはマウスとキーボードとを備える。入力装置174は、表示装置173に近接して配置される。操作者は、表示装置173を見ながら入力装置174を通してインタラクティブにMRI装置100の各種処理に必要な指示、データを入力する。
 CPU171は、操作者が入力した指示に従って、記憶装置172に予め保持されるプログラムを実行することにより、MRI装置100の動作の制御、各種データ処理等の制御処理部170の各処理を実現する。上記シーケンサ140への指令は、予め記憶装置172に保持されたパルスシーケンスに従って行われる。また、例えば、受信部160からのデータが制御処理部170に入力されると、CPU171は、信号処理、画像再構成処理等を実行し、その結果である被検体101の断層像を表示装置173に表示するとともに、記憶装置172に記憶する。
 送信コイル151と傾斜磁場コイル131とは、被検体101が挿入される静磁場発生部120の静磁場空間内に、垂直磁場方式であれば被検体101に対向して、水平磁場方式であれば被検体101を取り囲むようにして設置される。また、受信コイル161は、被検体101に対向して、或いは取り囲むように設置される。
 現在、MRI装置の撮像対象核種で、臨床で普及しているものは、被検体101の主たる構成物質である水素原子核(プロトン)である。MRI装置100では、プロトン密度の空間分布や、励起状態の緩和時間の空間分布に関する情報を画像化することで、人体頭部、腹部、四肢等の形態または機能を、二次元もしくは三次元的に撮像する。
 本実施形態の制御処理部170は、図2に示すように、非直交系サンプリング法のパルスシーケンスに従って、予め定めた計測空間(k空間)の複数の走査軌跡に沿ってエコー信号を計測し、データ列として走査軌跡上にそれぞれ配置する計測部210と、複数のエコー信号から画像を再構成し、再構成画像を得る画像再構成部220と、を備える。
 計測部210は、非直交系サンプリング法のパルスシーケンスに従って、静磁場発生部120、傾斜磁場発生部130、高周波磁場発生部150および高周波磁場検出部160の動作を制御し、エコー信号をサンプリングして得たデータ列をk空間の走査軌跡上に配置する。なお、本実施形態では、以後、各走査軌跡をブレードと呼ぶ。
 また、画像再構成部220は、k空間の各ブレード上のデータ列間の、低空間周波数領域の位相差を低減(縮小)し、再構成画像を得る。これを実現するため、本実施形態の画像再構成部220は、複数のブレード上のデータ列各々の位相の補正量(位相補正量)を算出する位相補正量算出部221と、算出した位相補正量を用い、複数のブレード上のデータ列各々の位相を補正する位相補正部222と、補正後の複数のブレード上のデータ列から再構成画像を算出する再構成部223と、を備える。位相補正量算出部221が算出する位相補正量は、複数のブレード上のデータ列間の位相差を低減するよう算出される。
 位相補正量算出部221は、複数のデータ列各々の、ブレード上の所定位置のデータの位相が揃うよう位相補正量を算出する。このとき、所定位置は、エコー信号の中心(ピーク位置)のk空間原点からのシフト量であるピークシフト量を反映して決定される。
 本実施形態では、所定の位置は、複数のブレードの交点とする。具体的には、複数のブレードの中の予め定めた基準とするブレード(基準ブレード)と、基準ブレード以外のブレードとの交点とし、位相補正量算出部221は、基準ブレード以外のブレード上のデータ列の交点における位相が、基準ブレード上のデータ列の交点における位相に合致するよう位相補正量を算出する。
 本実施形態の位相補正量算出部221は、基準ブレード以外のブレード上のデータ列各々について、ピークシフト量を用いて交点の位置を算出し、算出した交点の位置の情報を用いて交点から当該データ列の中点までの距離および基準ブレード上のデータ列(基準データ列)の中点までの距離をそれぞれ算出し、算出した各々の距離を用いて当該データ列の交点における位相と前記基準データ列の交点における位相とを算出する。
 上述のように、本実施形態では、位相補正量を算出する際、エコー信号のピークシフト量を用いる。このエコー信号のピーク位置のシフトであるピークシフトは、静磁場の不均一、傾斜磁場の出力誤差により発生する。また、エコー信号のピーク位置は、読み出し傾斜磁場パルス印加前に印加されるディフェーズパルスの面積に応じて変化する。ディフェーズ後、パルスシーケンスによって予め定められた読み出し傾斜磁場パルスの印加によってエコー信号の位相が揃うタイミングが、ピーク位置となるためである。
 エコー信号のずれ量(ピークシフト量)は、予めプリスキャンなどを行い、別途算出しておく。これにより、撮像時間が長引くのを防ぐことができる。算出は、例えば、特開2005-152175号公報に記載されているように、特定のエコー信号のみを取得するパルスシーケンスによりエコー信号を取得し、当該エコー信号を用いて行う。
 本実施形態の位相補正量算出の詳細の説明に先立ち、エコー信号のピークシフトによる各ブレードのシフトについて簡単に説明する。なお、以下、本明細書では、説明を簡単にするため、撮像断面がMRI装置100の装置座標系のXY面である場合を例にあげて説明する。
 図3(a)は、k空間上で、エコー信号のピークシフトによる、走査軌跡であるブレード301のシフトを説明するための図である。ここでは、簡単のため装置座標系のX軸,Y軸が、k空間上の座標系(計測座標系)のkx軸,ky軸と一致している場合について図示する.以下、本明細書の走査軌跡を図示する全図において同様とする。本図において、白丸の位置305が、ブレード301にシフト(位置ずれ)がない場合の、各ブレード301を構成するデータの位置であり、黒丸の位置304が、シフトがある場合(実際)のデータ位置である。例えば、kx軸上に配置されるブレード301のブレード中心301cは、シフトがなければ、k空間の原点の位置303となるが、シフトがあるため、実際は、301cの位置となる。
 本実施形態では、ブレード301のシフト量として、シフト後のブレード中心301cの、k空間座標上の座標Pdnの値(座標値)を算出する。上記静磁場の不均一、傾斜磁場の出力誤差等により発生するピークシフトによる、X軸方向のシフト量をΔD’x、Y軸方向のシフト量をΔD’yとすると、図3(b)に示すように、d番目(dは1以上の整数)のブレード301の、kx方向のシフト量ΔD’x_dおよびky方向のシフト量ΔD’y_dは、以下の式(1)および式(2)で表される。
   ΔD’x_d=ΔD’xcos(θd)・・・(1)
   ΔD’y_d=ΔD’ysin(θd)・・・(2)
 よって、d番目のブレード中心301cの、k空間上での座標Pdn{x,y}は、以下の式(3)で表される。
   Pdn{x,y}={ΔD’xcos(θd),ΔD’ysin(θd)}・・・(3)
 ここで、θdは、d番目のブレード301が、X軸と成す角度である。
 次に、本実施形態の位相補正量算出部221による位相補正量算出の詳細を、図面を用いて説明する。図4(a)および図4(b)は、本実施形態の位相補正量の算出手法を説明するための図である。図4(b)は、図4(a)のk空間中心付近の拡大図である。ここでは、1番目のブレード311をk空間のkx軸上に配置し、これを基準ブレード311とする。
 なお、基準ブレード311は、kx軸上のブレードに限られず、どのブレードであってもよい。
 b番目(bは2以上の整数)のブレード312と、基準ブレード311との交点を交点313とする。なお、本図において、白丸で示した位置315は、ブレードに位置ずれ(シフト)が無い場合の理想的なデータの位置、黒丸で示した位置314は、ブレードのシフトを反映した、実際のデータの位置である。
 シフトが無い場合、b番目のブレード312は、基準ブレード311とk空間の原点316で交差する。しかし、実際は、各ブレードにシフトがあるため、両者は、白バツ(クロス)で示した位置(交点)313で交差する。本実施形態では、b番目のブレード312の交点313における位相を、基準ブレードの交点313における位相に揃えるよう位相補正を行う。
 本実施形態の位相補正量算出部221は、交点313の座標を求め、交点313とそれぞれのブレード311、312の中点(ブレード中心)311c、312cとの距離ΔDb1およびΔDbbを求め、それぞれのブレード311、312の、交点313における位相値(Phase_1、Phase_b)を求める。そして、交点313における、b番目のブレード312の位相値の、基準ブレード311の交点313における位相値に対する差分(位相差分)を、位相補正量PhC_bとして算出する。
 交点313とそれぞれのブレード中心311c、312cとの距離ΔDb1及びΔDbbは、X軸方向およびY軸方向のピークシフト量ΔD’xおよびΔD’yとb番目のブレード312のブレード角度θbとから求める。なお、本実施形態では、ブレード角度θbは、b番目のブレード312が、X軸と成す角度とする。
 図4(b)および上記式(3)より、交点313の座標PIb{x、y}は、ピークシフト量ΔD’xとΔD’yとブレード角度θbとを用いて表すと、以下の式(4)および式(5)の関係を有する。
   y+ΔD’ysin(θb)=tan(θb){x+ΔD’xcos(θb)}・・・(4)
   y=0・・・(5)
 従って、交点313の座標PIb{x、y}は、b番目のブレード312のブレード角度θbとピークシフト量とを用い、以下の式(6)で表される。
   PIb{x,y}={cos(θb)・(ΔD’y-ΔD’x),0}・・・(6)
 交点313の座標値を用い、基準ブレード311の中点311cから交点313までの距離ΔDb1と、b番目のブレード312の中点312cから交点313までの距離ΔDbbとは、それぞれ、以下の式(7)、式(8)および式(9)で求められる。
   ΔDb1=ΔD’x{1-cos(θb)}+ΔD’ycos(θb) ・・・(7)
   ΔDbb=ΔD’y ・・・(8)
   ∵|ΔDbb|2=[cos(θb)・(ΔD’y-ΔD’x)+ΔD’xcos(θb)]2+[ΔD’ysin(θb)]2=(ΔD’y)2 ・・・(9)
 そして、算出した交点313の位置の実虚の信号値を補間により取得し、基準ブレード311の交点313における位相値Phase_1と、b番目のブレード312の交点313における位相値Phase_bとを、以下の式(10)および式(11)に従って得る。
 Phase_1=tan-1(Blade_1(Imaginary(CENTER+ΔDb1))/Blade_1(Real(N/2+ΔDb1)))・・・(10)
 Phase_b=tan-1(Blade_b(Imaginary(CENTER+ΔDbb))/Blade_b(Real(N/2+ΔDbb)))・・・(11)
 ここで、Blade_1()およびBlade_b()は、それぞれ、基準ブレード311およびb番目のブレード312のデータ列を、Real()は実部データを、Imaginary()は虚部データを、Nは各ブレードのデータ点数を、CENTERは中点311c,312cの位置を、それぞれ表す。
 従って、位相補正量(位相差分)PhC_bは、以下の式(12)で得られる。
   PhC_b=Phase_b-Phase_1・・・(12)
 この位相補正量PhC_bをb番目のブレード312のデータ列に適用し、その位相値を補正することにより、交点313における1番目のブレード311およびb番目のブレード312それぞれの位相値が合致する。
 なお、図5に示すように、b番目のブレード312の数は1以上である。また、ピークシフト(ブレード中心のシフト)があると、基準ブレード311と、それぞれのb番目のブレード312との交点は、1点に集中しない。このため、位相補正量算出部221は、各b番目のブレード312に関し、白バツで示す基準ブレードとの交点313の位置を算出し、基準ブレード311の中点311cとの距離ΔDb1と、当該ブレード312の中点312cとの距離ΔDbbとを、上記手順で算出し、各々の中点311c、312cの位相値Phase_1およびPhase_bを得、位相補正量PhC_bを算出する。
 位相補正部222は、上述したように、算出した位相補正量PhC_bを用いて、b番目のブレード312を構成する全データの位相を補正する。本実施形態では、ブレード312を構成する全てのデータに対し、同じ位相補正量を用いて位相補正を行う。すなわち、補正後のb番目のブレード312の複素データ列Blade_b(x)は、以下の式(13)で表される。
   Blade_b(x)=|Blade_b(x)|・exp(i(PhC_b))・・・(13)
 ここで、|Blade_b(x)|は、b番目のブレードの複素データ列の振幅を表す。またiは虚数単位である。
 再構成部223は、補正後の各ブレードのデータ列を用い、画像を再構成する。
 以下、本実施形態の位相補正量算出部221および位相補正部222による位相補正処理の流れを説明する。図6は、本実施形態の位相補正処理の処理フローである。ここでは、最初に取得したブレード(1番目のブレード)を、基準ブレードとする。bを、ブレードをカウントするカウンタに用いる。また、全ブレード数をM(Mは1以上の整数)とする。
 位相補正量算出部221は、カウンタbを初期化(b=1)する(ステップS1101)。
 次に、位相補正量算出部221は、計測部210が取得したブレードのデータ列を基準ブレードデータとして保存する(ステップS1102)。そして、カウンタbを1インクリメント(b=1+1)する(ステップS1103)。
 位相補正量算出部221は、全ブレードについて位相補正処理を終えたか判別する(ステップS1104)。ここでは、カウンタbがブレード数Mを越えたか否かを判別する。
全てのブレード312について処理を終えたと判別された場合、処理を終了する。
 一方、未処理のブレード312があると判別された場合、位相補正量算出部221は、計測部210が取得したb番目のブレード312について、まず、基準ブレード311との交点313の座標PIb{x,y}を算出する(ステップS1105)。
 次に、位相補正量算出部221は、基準ブレード311の交点313の位相Phase_1と、b番目のブレード312の交点313の位相Phase_bと、をそれぞれ算出する(ステップS1106)。そして、両位相差として、b番目のブレードの位相補正量PhC_bを算出する(ステップS1107)。
 そして、位相補正部222は、算出した位相補正量PhC_bを用い、上記式(13)に従って、b番目のブレードの位相を補正する(ステップS1108)。位相補正量算出部221は、ステップS1103へもどり、処理を繰り返す。
 以上説明したように、本実施形態によれば、非直交系パルスシーケンスに従って取得した複数のブレードのデータ列について、ブレードとブレードとの交点の位相が揃うよう位相補正量を決定し、位相補正を行う。具体的には、複数のブレードの中で基準ブレードを定め、他のブレード各々について基準ブレードとの交点において位相が合致するよう位相補正量を決定し、当該位相補正量で各ブレードを構成する全データの位相を補正する。これにより、ブレードとブレードとが重なりあう箇所において、位相差が低減し、それによる信号の打ち消しあいも低減する。このとき、ブレードの受信エコーの中心のずれを考慮した位相変化量で位相を補正する。
 また、一般に非直交サンプリング法では、各ブレードは、低空間周波数領域で基準ブレードと交差する。従って、ブレード毎に基準ブレードとの交点の位相を基準ブレードの位相に合致させることにより、画質に最も影響する低空間周波数領域での位相差が低減する。これにより、画質劣化を低減させることができる。
 なお、本実施形態では、交点313での位相を合わせる際、交点1点のみの補間値を用いて合わせているが、これに限られない。例えば、交点313とその周辺の点の平均値などを用いて位相を合わせても良い。
 また、本実施形態では、k空間で位相補正を行っているが、位相補正を行う空間は、画像空間であってもよい。画像空間で行う場合は、以下の式(14)に従う。
   FT[Blade_b(x)]=|FT[Blade_1(x)]|・exp(i(PhC_b))・・・(14)
 ここでFT[]はフーリエ変換を表す。
 なお、本実施形態では、上述のように、基準ブレードのデータを保持する必要がある。
このとき、保持するデータは、基準ブレードの全データであってもよいし、また、使用メモリを削減するため、基準ブレードの原点周りの数点分のデータであってもよい。
 <<第二の実施形態>>
 次に、本発明の第二の実施形態を説明する。第一の実施形態では、基準ブレードを決定し、その他のブレードについて、基準ブレードとの交点で基準ブレードの位相と合致するよう、位相を補正している。これに対し、本実施形態では、後述する各ブレードの原点オフセット位置PObにおける位相が合致するよう補正を行う。従って、本実施形態では、基準ブレードの設定は、不要である。
 本実施形態のMRI装置は、基本的に第一の実施形態のMRI装置100と同様の構成を有する。ただし、上述のように位相補正量の算出手法が異なる。従って、位相補正量算出部221の処理が異なる。以下、本実施形態について、第一の実施形態と異なる構成に主眼をおいて説明する。なお、本実施形態においても、計測部210は、非直交系のサンプリング法のパルスシーケンスに従って、予め定めた計測空間(k空間)の複数の走査軌跡に沿ってエコー信号を計測し、データ列として走査軌跡上にそれぞれ配置する。
 以下、本実施形態の位相補正量算出部221による位相補正量算出処理を説明する。本実施形態では、第一の実施形態同様、複数のデータ列各々の、ブレード上の所定位置のデータの位相が揃うよう位相補正量を算出する。所定位置は、k空間の原点から前記複数のブレード各々に下ろした垂線の交点である原点オフセット位置であり、位相補正量算出部221は、複数のデータ列の原点オフセット位置の位相が、全て揃うよう位相補正量を算出する。
 本実施形態の位相補正量算出部221は、複数のデータ列各々について、ピークシフト量を用いて、オフセット位置を算出し、算出したオフセット位置の情報を用いて、当該オフセット位置から当該データ列の中点までの距離を得、当該距離を用い、当該データ列のオフセット位置における位相を算出する。
 なお、本実施形態においても、第一の実施形態同様、エコー信号のずれ量(ピークシフト量)は、予めプリスキャンなどを行い、別途算出しておく。
 以下、本実施形態の位相補正量算出部221による位相補正量算出処理について、図を用いて説明する。
 図7(a)は、本実施形態のb番目のブレード322と、kx軸およびky軸との位置関係を示したものである。また、図7(b)は、図7(a)のk空間中心付近の拡大図である。
 なお、本実施形態では、bは1以上の整数である。本図において、白丸で示した位置325は、ブレード322に位置ずれ(シフト)が無い場合の理想的なデータの位置、黒丸で示した位置324は、ブレード322のシフトを反映した、実際のデータの位置である。
 各ブレード322の中点(ブレード中心)は、理論上は、白丸で示すk空間の原点326に位置する。しかしながら、実際には、傾斜磁場の出力応答のためシフトし、黒丸で示す位置322cとなる。
 上述のように、k空間の原点からb番目のブレード322に垂線を下ろした際の交点323を原点オフセット位置と呼ぶ。図中では、白バツで示す。本実施形態では、各ブレード322の原点オフセット位置の位相が揃うよう、各ブレード322の位相を補正する。
 各ブレード322の、原点オフセット位置323の位相は、各ブレード322のブレード中心322cと、原点オフセット位置323とのk空間上での距離(誤差)ΔDbを求め、算出する。なお、以下、本実施形態では、原点オフセット位置323とのk空間上での距離ΔDbを、周波数方向のシフト量と呼ぶ。
 各b番目のブレード322の中点(ブレード中心)322cの座標Pbnは、第一の実施形態の上記式(3)同様、以下の式(15)で表される。
   Pbn{x,y}={-ΔD’xcos(θb),-ΔD’ysin(θb)}・・・(15)
 ΔD’x、ΔD’yは、第一の実施形態同様、X軸方向およびY軸方向のブレードのシフト量である。また、θbは、b番目のブレード322がX軸と成す角度である。
 また、原点オフセット323の座標PObの座標値{x,y}は、以下の式(16)および式(17)の関係を有する。
   y+ΔD’ysin(θb)=tan(θb){x+ΔD’xcos(θb)}・・・(16)
   y=-x/tan(θb)・・・(17)
 従って、原点オフセット位置323の座標PObは、以下の式(18)で表される。
 POb{x,y}={cos(θb)・sin2(θb)・(ΔD’y-ΔD’x),cos2(θb)・sin(θb)・(ΔD’x-ΔD’y)}・・・(18)
 b番目のブレード322の、周波数方向のシフト量ΔDbは、原点オフセット位置323の座標位置PObとb番目のブレード322のブレード中心322cの座標Pbnとを用い式(20)から、以下の式(19)のように算出される。
 ΔDb=ΔD’xcos2(θb)+ΔD’ysin2(θb)・・・(19)
 ∵|ΔDb|2=[cos(θb)・sin2(θb)・(ΔD’y-ΔD’x)+ΔD’xcos(θb)]2+[cos2(θb)・sin(θb)・(ΔD’x-ΔD’y)+ΔD’ysin(θb)]2
     =(ΔD’x)2・cos4(θb)+2ΔD’x・ΔD’ycos2(θb)・sin2(θb)+(ΔD’y)2・sin4(θb)
     =[ΔD’xcos2(θb)+ΔD’ysin2(θb)]2・・・(20)
 そして、第一の実施形態同様、算出した原点オフセット位置323の実虚の信号値を補間により取得し、b番目のブレード322の原点オフセット位置323の位相値Phase_bを算出する。Phase_bは、以下の式(21)で算出される。
 Phase_b=tan-1(Blade_b(Imaginary(CENTER+ΔDb))/Blade_b(Real(CENTER+ΔDb)))・・・(21)
 ここで、Blade_b()は、b番目のブレード322のデータ列を、Real()は実部データを、Imaginary()は虚部データを、Nは各ブレードのデータ点数を、CENTERは中点の位置を、それぞれ表す。
 本実施形態では、全てのブレード322の原点オフセット位置323の位相値を揃える。従って、例えば、α[rad]に揃えるとすると、位相補正量PhC_bは、以下の式(22)となる。
   PhC_b=Phase_b-α・・・(22)
 αには、例えば、0を用いる。
 この位相補正量PhC_bをb番目のブレード322のデータ列に適用することにより、各ブレード322の原点オフセット位置323おける位相値が合致する。
 なお、図8に示すように、b番目のブレード322は複数ある。また、ブレード中心にシフトがあるため、各ブレード322の原点オフセット位置323は、1点に集中しない。このため、本実施形態においても、位相補正量算出部221は、各b番目のブレードに関し、原点オフセット位置323の位相値Phase_bを得、位置補正量を算出する。
 なお、本実施形態においても、位相補正部222は、算出した位相補正量PhC_bを用いて、b番目のブレード322を構成する全データの位相を補正する。補正後のb番目のブレード322の複素データ列Blade_b(x)は、第一の実施形態同様、以下の式(23)で表される。
   Blade_b(x)=|Blade_b(x)|・exp(i(PhC_b))・・・(23)
 ここで、|Blade_b(x)|は、b番目のブレードの複素データ列の振幅を表す。またiは虚数単位である。
 以下、本実施形態の位相補正量算出部221および位相補正部222による位相補正処理の流れを説明する。図9は、本実施形態の位相補正処理の処理フローである。bを、ブレードをカウントするカウンタに用いる。全ブレード数はM(Mは1以上の整数)とする。
 位相補正量算出部221は、カウンタbを初期化(b=1)する(ステップS1201)。
 本実施形態では、次に、位相補正量算出部221は、全ブレードについて位相補正処理を終えたか判別する(ステップS1202)。ここでは、カウンタbがブレード数Mを越えたか否かを判別する。全てのブレード312について処理を終えたと判別された場合、処理を終了する。
 一方、未処理のブレード322があると判別された場合、位相補正量算出部221は、計測部210が取得したb番目のブレード322について、周波数方向のシフト量ΔDbを算出する(ステップS1203)。位相補正量算出部221は、算出した周波数方向のシフト量を用い、原点オフセット位置323の位相を算出し(ステップS1204)、b番目のブレード322の位相補正量PhC_bを算出する(ステップS1205)。
 位相補正部222は、算出した位相補正量PhC_bを用い、上記式(23)に従って、b番目のブレードを補正する(ステップS1206)。位相補正量算出部221は、カウンタbを1インクリメントし(ステップS1207)、ステップS1202へもどり、処理を繰り返す。
 以上説明したように、本実施形態によれば、非直交系パルスシーケンスに従って取得した複数のブレードのデータ列について、各ブレードの原点オフセット位置における位相が揃うよう、各ブレードのデータ列の位相を補正する。これにより低空間周波数領域において、ブレード間の位相ずれが解消し、位相ずれにより発生する信号の打ち消し合いを低減することができる。従って、画質劣化を抑えることができる。
 さらに、本実施形態によれば、基準とするブレードを選択する必要がない。従って、基準とするブレードのデータを保持するためのメモリ領域の確保が不要となる。さらに、位相補正結果が、基準とするブレードの精度の影響を受けない。
 なお、本実施形態においても、第一の実施形態同様、位置324の周辺の位置の平均値を用いて位相を合わせてもよい。
 また、本実施形態においても、第一の実施形態同様、位相補正は、画像空間で行ってもよい。
 <<第三の実施形態>>
 次に、本発明の第三の実施形態を説明する。本実施形態では、オフセンタ撮像時の受信周波数の変化から、計算によって各ブレード中心の位相値を算出し、位相補正量を決定する。ただし、本実施形態では、静磁場の不均一がなく、傾斜磁場波形の面積誤差が時間と共に変化しないものとする。
 本実施形態のMRI装置は、基本的に第一の実施形態のMRI装置100と同様の構成を有する。ただし、上述のように位相補正量の算出手法が異なる。従って、位相補正量算出部221の処理が異なる。以下、本実施形態について、第一の実施形態と異なる構成に主眼をおいて説明する。本実施形態においても、計測部210は、非直交系サンプリング法のパルスシーケンスに従って、予め定めたk空間の複数の走査軌跡に沿ってそれぞれエコーを計測し、データ列として当該走査軌跡(ブレード)上に配置する。このとき、受信周波数を変更することにより、オフセンタ撮像を行う。
 本実施形態の位相補正量算出部221による位相補正量算出処理を説明する。本実施形態の位相補正量算出部221は、k空間上のオフセンタ距離を用いて複数のデータ列各々の中点の位相変化量を算出し、算出した位相変化量を、位相補正量とする。すなわち、オフセンタ撮像時の、エコー信号取得間(A/D間)の位相変化量から、ブレード中心の位相値を算出する。
 静磁場の不均一がなく、傾斜磁場波形の面積誤差が時間と共に変化しないとの前提で、受信周波数をシフトさせてオフセンタ撮像を行うと、受信周波数のシフトにより、各ブレード内で位相が回転する。この位相回転量は、オフセンタ距離、受信帯域幅および視野サイズから計算できる。本実施形態の位相補正量算出部221は、算出した位相回転量を位相補正量とする。具体的な算出手順は、以下のとおりである。
 まずオフセンタ撮像時の受信周波数Rfを求める。受信周波数Rfは、以下の式(24)で求められる。
   Rf=BW・OffcD/FOV・・・(24)
ここで、BWは受信周波数帯域[Hz]、OffcDはオフセンタ距離[m]、FOVは視野サイズ[m]である。
 検出中(A/D中)の位相は、オフセンタ撮像の場合、図10に示すように変化する。従って、ブレード内の位相回転量(位相変化量)は、以下の(25)で表される。
   ΔΦf(n)=Rf・π・n/BW・・・(25)
 ここで、nはブレードを構成するデータ列の要素番号である。
 従って、ブレード中心の、k空間のkx軸方向の位相変化量ΔΦxと、ky軸方向の位相変化量ΔΦyは、それぞれ、以下の式(26)と式(27)で表される。
   ΔΦx=OffcD_X/FOV_X・π・CENTER     ・・・(26)
   ΔΦy=OffcD_Y/FOV_Y・π・CENTER     ・・・(27)
 ここで、OffcD_XおよびOffcD_Yは、それぞれ、X軸方向およびY軸方向のオフセンタ距離[m]、FOV_XおよびFOV_Yは、それぞれ、X軸方向およびY軸方向の視野サイズ[m]、CENTERは中点の位置である。
 これを用い、b番目のブレード中心の位相変化量ΔΦbは、以下の式(28)で表される。
   ΔΦb=ΔΦx・cos(θb)-ΔΦy・sin(θb)・・・(28)
 本実施形態の位相補正量算出部221は、以上の手順により受信周波数のシフト量に応じて、各ブレードの、ブレード中心の位相変化量ΔΦbを求め、これが全て揃うよう、位相補正量PhC_bを決定する。例えば、ブレード中心の位相が、全てα〔rad〕になるとすると、位相補正量PhC_bは、以下の式(29)で表される。
   PhC_b=ΔΦb-α・・・(29)
 αには、例えば、0を用いる。
 本実施形態においても、算出した位相補正量PhC_bを用いて、b番目のブレードを構成する全データの位相を補正する。補正後のb番目のブレード322の複素データ列Blade_b(x)は、第一の実施形態同様、以下の式(30)で表される。
   Blade_b(x)=|Blade_b(x)|・exp(i(PhC_b))・・・(30)
 以上説明したように、本実施形態では、各ブレードのブレード中心の位相回転量を計算により求め、ブレード間の当該位相差を低減するよう各ブレードの位相を補正する。これにより、ブレード間の位相ずれが低減する。従って、ブレード間の位相ずれにより発生する画質劣化を抑えることができる。
 なお、本実施形態においても、第一の実施形態同様、位相補正は、画像空間で行ってもよい。
 また、本実施形態では、エコー信号のピークシフト量を用いていないが、静磁場不均一、傾斜磁場波形の面積誤差による、ピークシフトを考慮することもできる。この場合の各ブレードの、ブレード中心における位相変化量ΔΦ’xおよびΔΦ’yは、以下の式(31)および式(32)で表される。
   ΔΦ'x=OffcD_X/FOV_X・π・(CENTER+ΔD'x)・・・(31)
   ΔΦ'y=OffcD_Y/FOV_Y・π・(CENTER+ΔD'y)・・・(32)
 従って、b番目のブレードのブレード中心の位相変化量ΔΦ’bは、以下の式(33)のように表される。
   ΔΦ'b=ΔΦ'x・cos(θb)-ΔΦ'y・sin(θb)・・・(33)
 ピークシフトを考慮する場合は、この位相変化量ΔΦ’bを用いて、位相補正量PhC_bを決定する。
 なお、上記各実施形態では、撮像断面が装置座標系のXY面である場合を例にあげて説明したが、撮像断面は任意の断面であってよい。この場合は、装置座標系の各軸におけるピークシフト量を各ブレードの角度に応じて計測座標系に展開し、基準ブレードとの交点からそれぞれのブレードの中点までの距離(ずれ量)ΔDb1およびΔDbbを算出する。
 具体的な展開、算出手法は以下のとおりである。
 XYZ軸それぞれのピークシフト量を反映させた計測座標系における座標kRは、以下の式(34)で表される。
Figure JPOXMLDOC01-appb-I000001
 ここで、kRAは装置座標系における座標、GRR(bl)はXYZ各軸の読み出し傾斜磁場[T]、dAはXYZ各軸のピークシフト量、bはブレード番号、θbはブレード番号がbのブレード(ブレードb)の、kx軸に対するブレード角度、nはブレード内のデータ点番号、Ncは中点の位置のデータ番号、eは単位行列、γは磁気回転比[rad/T]、Δtはデータ列のサンプリング間隔[s]、Gは傾斜磁場強度[T]である。また、ROMは計測座標系を装置座標系に変換するための回転行列であり、以下の式(35)で定義される。
Figure JPOXMLDOC01-appb-I000002
 ここで,sx,sy,szは計測座標系におけるスライス軸方向の傾斜磁場の、装置座標系のXYZ軸への投影成分、px,py,pzは計測座標系における位相エンコード軸(ky)方向の傾斜磁場の装置座標系のXYZ軸への投影成分、fx,fy,fzは計測座標系における読み出し傾斜磁場軸(kx)方向の傾斜磁場の装置座標系のXYZ軸への投影成分である。
 なお、上記式(34)のKR(b、n)は、b番目のブレードのn番目のデータ点の、上記計測座標系の座標である。以下、同様とする。
 従って、ゼロ位相エンコードデータの計測座標系の座標kS,kP,kFは、以下の式(36)で求まる。なお、以降では、(γ/2π)×Δt×G=1と規格化して扱う。
Figure JPOXMLDOC01-appb-I000003
 ブレードを2次元平面(kP軸-kF軸平面)に配置する場合,スライス方向の座標(kS軸)は考慮する必要がない.そのため,以下kP軸及びkF軸のみ扱う。1番目のブレード(b=1)におけるゼロ位相エンコードの中点D′0n(n=Nc)の座標(kF,kP)は、以下の式(37)で表される。
Figure JPOXMLDOC01-appb-I000004
 また、b番目のブレードにおけるゼロ位相エンコードの中点D′bn(n=Nc)の計測座標系における座標(kF,kP)は、以下の式(38)で表される。
Figure JPOXMLDOC01-appb-I000005
 よって、b番目のブレードの計測座標系における式は、その傾きがtanθbであるため、以下の式(39)で表される。
Figure JPOXMLDOC01-appb-I000006
 1番目のブレードとの交点PIbの座標(kF,kP)は、1番目のブレード(b=1)におけるゼロ位相エンコードの中点D′0nが上記式(37)で表されるため、θb≠πとして、以下の式(40)で求まる。
Figure JPOXMLDOC01-appb-I000007
 従って、交点PIbの、1番目のブレードの中点D′0nからのずれ量ΔDb1は、以下の式(41)で表される。
Figure JPOXMLDOC01-appb-I000008
 また、交点PIbの、b番目のブレードの中点D′bnからのずれ量|ΔDbb|は、以下の式(42)で求まる。
Figure JPOXMLDOC01-appb-I000009
 符号を考慮すると、交点Ibの、b番目のブレードの中点D′bnからのずれ量ΔDbbは、以下の式(43)で表される。
Figure JPOXMLDOC01-appb-I000010
 以上のように算出したずれ量ΔDb1およびΔDbbを用いて式(10)~(13)を計算し、位相補正を行う。
 なお、本発明の実施形態は、上記形態に限られない。例えば、上記各実施形態では、非直交サンプリング法の中でもk空間を放射状に走査するラディアルサンプリング法を例にあげて図示し、説明しているが、用いるサンプリング法はこれに限られない。k空間上で各ブレードが重なり合う軌跡を描くサンプリング法であればよく、例えば、ラディアルサンプリング法に位相エンコードを組み合わせた、ハイブリッドラディアル法であってもよい。
 ハイブリッドラディアル法の場合、平行な複数の直線軌跡の中の、中心の一つの軌跡を用いて、上記各実施形態の手法で位相補正量PhC_bを算出し、この位相補正量PhC_bを用いて1のブレードを構成する全ての平行な複数の直線軌跡上のデータの位相を補正する。
 100 MRI装置、101 被検体、120 静磁場発生部、130 傾斜磁場発生部、131 傾斜磁場コイル、132 傾斜磁場電源、140 シーケンサ、150 高周波磁場発生部、150 送信部、151 送信コイル、152 高周波発振器、153 変調器、154 高周波増幅器、160 高周波磁場検出部、160 受信部、161 受信コイル、162 信号増幅器、163 直交位相検波器、164 A/D変換器、170 制御処理部、171 CPU、172 記憶装置、173 表示装置、174 入力装置、210 計測部、220 画像再構成部、221 位相補正量算出部、222 位相補正部、223 再構成部、301 ブレード、301c ブレード中心、303 位置、304 位置、305 位置、311 基準ブレード、311c 中点、312 ブレード、312c 中点、313 交点、314 位置、315 位置、316 原点、322 ブレード、322c ブレード中心、323 原点オフセット位置、324 位置、325 位置、326 原点

Claims (12)

  1.  非直交系サンプリング法のパルスシーケンスに従って、予め定めたk空間の複数の走査軌跡に沿ってそれぞれエコーを計測し、データ列として当該走査軌跡上に配置する計測部と、
     前記複数の走査軌跡上のデータ列から画像を再構成し、再構成画像を得る画像再構成部と、を備え、
     前記画像再構成部は、
     前記複数の走査軌跡上のデータ列各々の位相補正量を算出する位相補正量算出部と、
     前記算出した位相補正量を用い、前記複数の走査軌跡上のデータ列各々の位相を補正する位相補正部と、
     前記補正後の前記複数の走査軌跡上のデータ列から前記再構成画像を生成する再構成部と、を備え、
     前記位相補正量算出部は、前記複数の走査軌跡上のデータ列間の位相差を低減するよう前記位相補正量を算出すること
     を特徴とする磁気共鳴イメージング装置。
  2.  請求項1記載の磁気共鳴イメージング装置であって、
     前記位相補正量算出部は、前記複数の走査軌跡各々の所定位置の位相が揃うよう前記位相補正量を算出し、
     前記所定位置は、前記エコーのエコー中心のk空間原点からのシフト量であるピークシフト量を反映して決定されること
     を特徴とする磁気共鳴イメージング装置。
  3.  請求項2記載の磁気共鳴イメージング装置であって、
     前記所定位置は、前記複数の走査軌跡の中の予め定めた基準とする走査軌跡である基準走査軌跡と、その他の走査軌跡との交点であり、
     前記位相補正量算出部は、前記その他の走査軌跡上のデータ列の前記交点における位相が、前記基準走査軌跡上のデータ列の前記交点における位相に合致するよう前記位相補正量を算出すること
     を特徴とする磁気共鳴イメージング装置。
  4.  請求項2記載の磁気共鳴イメージング装置であって、
     前記所定位置は、k空間の原点から前記複数の走査軌跡各々に下ろした垂線の交点である原点オフセット位置であり、
     前記位相補正量算出部は、前記複数のデータ列の前記原点オフセット位置の位相が、全て揃うよう前記位相補正量を算出すること
     を特徴とする磁気共鳴イメージング装置。
  5.  請求項1記載の磁気共鳴イメージング装置であって、
     前記計測部は、オフセンタ撮像を行い、
     前記位相補正量算出部は、k空間上のオフセンタ距離を用いて前記複数のデータ列各々の中点の位相変化量を算出し、当該位相変化量を、前記位相補正量とすること
     を特徴とする磁気共鳴イメージング装置。
  6.  請求項3記載の磁気共鳴イメージング装置であって、
     前記位相補正量算出部は、前記その他の走査軌跡上のデータ列各々について、前記ピークシフト量を用いて前記交点の位置を算出し、算出した交点の位置の情報を用いて前記交点から当該データ列の中点までの距離および前記基準走査軌跡上のデータ列の中点までの距離をそれぞれ得、各々の前記距離を用い、当該データ列の前記交点における位相と前記基準走査軌跡上のデータ列の前記交点における位相とを算出すること
     を特徴とする磁気共鳴イメージング装置。
  7.  請求項4記載の磁気共鳴イメージング装置であって、
     前記位相補正量算出部は、前記複数のデータ列各々について、前記ピークシフト量を用いて前記オフセット位置を算出し、算出した前記オフセット位置の情報を用いて、当該オフセット位置から当該データ列の中点までの距離を得、当該距離を用い、当該データ列の前記オフセット位置における位相を算出すること
     を特徴とする磁気共鳴イメージング装置。
  8.  請求項5記載の磁気共鳴イメージング装置であって、
     前記位相補正量算出部は、前記複数のデータ列各々の中点の前記位相変化量に、前記ピークシフト量を反映させること
     を特徴とする磁気共鳴イメージング装置。
  9.  請求項1記載の磁気共鳴イメージング装置であって、
     前記非直交系サンプリング法は、ラディアルサンプリング法、ハイブリッドラディアル法のいずれかであること
     を特徴とする磁気共鳴イメージング装置。
  10.  非直交系サンプリング法のパルスシーケンスに従って、予め定めたk空間の複数の走査軌跡それぞれに沿って計測したエコーから得、当該k空間の前記走査軌跡上に配置されたデータ列各々について、各データ列間のk空間の低空間周波数帯域の位相差を低減するよう、前記データ列各々の位相を補正する位相補正処理を行い、
     補正後のエコー信号から画像を再構成すること
     を特徴とする磁気共鳴イメージング方法。
  11.  請求項10記載の磁気共鳴イメージング方法であって、
     前記複数のデータ列の中の1のデータ列を基準データ列として保存し、
     その他のデータ列各々について、当該データ列の走査軌跡と前記基準データ列の走査軌跡との交点を算出し、当該データ列の前記交点の位相と前記基準データ列の前記交点の位相とを算出し、前記両位相の位相差を位相補正量として算出し、算出した位相補正量で、当該データ列の位相を補正することにより、前記位相補正処理を行うこと
     を特徴とする磁気共鳴イメージング方法。
  12.  請求項10記載の磁気共鳴イメージング方法であって、
     前記複数のデータ列各々について、当該データ列の中点の原点オフセット位置の位相を算出し、算出した前記位相が所定の値となるよう位相補正量を算出し、算出した位相補正量で当該データ列の位相を補正することにより前記位相補正処理を行い、
     前記原点オフセット位置は、k空間の原点から前記走査軌跡に下ろした垂線の交点位置であること
     を特徴とする磁気共鳴イメージング方法。
PCT/JP2014/062306 2013-05-17 2014-05-08 磁気共鳴イメージング装置および磁気共鳴イメージング方法 WO2014185323A1 (ja)

Priority Applications (3)

Application Number Priority Date Filing Date Title
JP2015517044A JP6410715B2 (ja) 2013-05-17 2014-05-08 磁気共鳴イメージング装置および磁気共鳴イメージング方法
CN201480022131.8A CN105120745B (zh) 2013-05-17 2014-05-08 磁共振成像装置以及磁共振成像方法
US14/785,454 US10048343B2 (en) 2013-05-17 2014-05-08 Magnetic resonance imaging apparatus and magnetic resonance imaging method

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
JP2013-104955 2013-05-17
JP2013104955 2013-05-17

Publications (1)

Publication Number Publication Date
WO2014185323A1 true WO2014185323A1 (ja) 2014-11-20

Family

ID=51898301

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2014/062306 WO2014185323A1 (ja) 2013-05-17 2014-05-08 磁気共鳴イメージング装置および磁気共鳴イメージング方法

Country Status (4)

Country Link
US (1) US10048343B2 (ja)
JP (1) JP6410715B2 (ja)
CN (1) CN105120745B (ja)
WO (1) WO2014185323A1 (ja)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2021515648A (ja) * 2018-03-13 2021-06-24 コーニンクレッカ フィリップス エヌ ヴェKoninklijke Philips N.V. スパイラル獲得によるmr像形成

Families Citing this family (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP6410715B2 (ja) * 2013-05-17 2018-10-24 株式会社日立製作所 磁気共鳴イメージング装置および磁気共鳴イメージング方法
JP6464088B2 (ja) * 2013-09-03 2019-02-06 株式会社日立製作所 磁気共鳴イメージング装置および磁気共鳴イメージング方法
WO2016125572A1 (ja) * 2015-02-06 2016-08-11 株式会社日立製作所 磁気共鳴イメージング装置および磁気共鳴イメージング方法
KR102259846B1 (ko) * 2018-07-03 2021-06-03 가천대학교 산학협력단 자기공명 영상장치의 기계 학습 기반의 경사자계 오차 보정 시스템 및 방법
CN110244246B (zh) * 2019-07-03 2021-07-16 上海联影医疗科技股份有限公司 磁共振成像方法、装置、计算机设备和存储介质
US10884086B1 (en) * 2019-07-29 2021-01-05 GE Precision Healthcare LLC Systems and methods for accelerated multi-contrast propeller
DE102020209382A1 (de) * 2020-07-24 2022-01-27 Siemens Healthcare Gmbh Verfahren zur Aufnahme von Messdaten mittels einer Magnetresonanzanlage mit einer Korrektur der verwendeten k-Raumtrajektorien

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2005021691A (ja) * 2003-07-02 2005-01-27 Ge Medical Systems Global Technology Co Llc 位相エンコード配置のためのシステム及び方法
WO2005023108A1 (ja) * 2003-09-05 2005-03-17 Hitachi Medical Corporation 磁気共鳴イメージング装置
WO2007013423A1 (ja) * 2005-07-27 2007-02-01 Hitachi Medical Corporation 磁気共鳴イメージング装置
WO2008152937A1 (ja) * 2007-06-14 2008-12-18 Hitachi Medical Corporation 磁気共鳴イメージング装置及び傾斜磁場に起因する誤差補正方法
WO2009093517A1 (ja) * 2008-01-23 2009-07-30 Hitachi Medical Corporation 磁気共鳴イメージング装置及びマルチコントラスト画像取得方法

Family Cites Families (14)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6188219B1 (en) * 1999-01-22 2001-02-13 The Johns Hopkins University Magnetic resonance imaging method and apparatus and method of calibrating the same
US7023207B1 (en) * 2005-02-16 2006-04-04 General Electric Company Method and system of MR imaging with reduced radial ripple artifacts
DE102005046732B4 (de) * 2005-04-18 2010-04-15 Siemens Ag Verbessertes Rekonstruktionsverfahren bei der Propellerbildgebung in der Magnetresonanztomographie
DE102005019214B4 (de) * 2005-04-25 2007-01-11 Siemens Ag Kalibrier-Verfahren zur artefaktreduzierten MRT-Bildgebung bei Verschiebung des FOV sowie Computersoftwareprodukt
US7382127B2 (en) * 2006-09-15 2008-06-03 General Electric Company System and method of accelerated MR propeller imaging
US7482806B2 (en) * 2006-12-05 2009-01-27 Siemens Aktiengesellschaft Multi-coil magnetic resonance data acquisition and image reconstruction method and apparatus using blade-like k-space sampling
DE102007044463B4 (de) * 2007-09-18 2009-05-14 Bruker Biospin Mri Gmbh Verfahren zur Bestimmung der räumlichen Verteilung von Magnetresonanzsignalen durch mehrdimensionale HF-Anregungspulse
JP5575385B2 (ja) * 2007-11-02 2014-08-20 株式会社東芝 磁気共鳴イメージング装置
EP2526439B1 (en) * 2010-01-18 2014-07-30 Max-Planck-Gesellschaft zur Förderung der Wissenschaften e.V. Method and device for magnetic resonance spectroscopic imaging
WO2012043311A1 (ja) * 2010-09-27 2012-04-05 株式会社 日立メディコ 磁気共鳴イメージング装置および磁気共鳴イメージング方法
DE102011077197B4 (de) * 2011-06-08 2013-05-16 Siemens Aktiengesellschaft Verzeichnungskorrektur bei einer Magnetresonanz-Bildgebung
WO2013047275A1 (ja) * 2011-09-29 2013-04-04 株式会社 日立メディコ 磁気共鳴イメージング装置および磁気共鳴イメージング方法
DE102012204434B3 (de) * 2012-03-20 2013-07-11 Siemens Aktiengesellschaft Mehrschicht-MRI-Anregung mit simultaner Refokussierung aller angeregten Schichten
JP6410715B2 (ja) * 2013-05-17 2018-10-24 株式会社日立製作所 磁気共鳴イメージング装置および磁気共鳴イメージング方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2005021691A (ja) * 2003-07-02 2005-01-27 Ge Medical Systems Global Technology Co Llc 位相エンコード配置のためのシステム及び方法
WO2005023108A1 (ja) * 2003-09-05 2005-03-17 Hitachi Medical Corporation 磁気共鳴イメージング装置
WO2007013423A1 (ja) * 2005-07-27 2007-02-01 Hitachi Medical Corporation 磁気共鳴イメージング装置
WO2008152937A1 (ja) * 2007-06-14 2008-12-18 Hitachi Medical Corporation 磁気共鳴イメージング装置及び傾斜磁場に起因する誤差補正方法
WO2009093517A1 (ja) * 2008-01-23 2009-07-30 Hitachi Medical Corporation 磁気共鳴イメージング装置及びマルチコントラスト画像取得方法

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2021515648A (ja) * 2018-03-13 2021-06-24 コーニンクレッカ フィリップス エヌ ヴェKoninklijke Philips N.V. スパイラル獲得によるmr像形成
JP7209007B2 (ja) 2018-03-13 2023-01-19 コーニンクレッカ フィリップス エヌ ヴェ スパイラル獲得によるmr像形成

Also Published As

Publication number Publication date
CN105120745A (zh) 2015-12-02
US20160054415A1 (en) 2016-02-25
US10048343B2 (en) 2018-08-14
JP6410715B2 (ja) 2018-10-24
JPWO2014185323A1 (ja) 2017-02-23
CN105120745B (zh) 2017-12-19

Similar Documents

Publication Publication Date Title
JP6410715B2 (ja) 磁気共鳴イメージング装置および磁気共鳴イメージング方法
US7622926B2 (en) Magnetic resonance imaging apparatus
US7372269B2 (en) Magnetic resonance imaging method and apparatus
JP6464088B2 (ja) 磁気共鳴イメージング装置および磁気共鳴イメージング方法
US8473028B2 (en) K-space sample density compensation for magnetic resonance image reconstruction
CN106419917B (zh) 磁共振成像中用体积导航做前瞻性运动校正的方法和装置
JP2000060821A (ja) 系統誤差を補正する方法及びシステム
US9513356B2 (en) Magnetic resonance imaging apparatus and reconstructed image acquisition method
JP5683984B2 (ja) 磁気共鳴イメージング装置および非線形性歪み補正方法
EP4092439A1 (en) System and method for reconstruction using a high-resolution phase in magnetic resonance images
JP5978430B2 (ja) 磁気共鳴イメージング装置、及びエコー信号計測方法
WO2018020905A1 (ja) 磁気共鳴イメージング装置及びその制御方法
JP4707558B2 (ja) 磁気共鳴イメージング装置
US20040160221A1 (en) Method to excite planar slices in a magnetic resonance tomography device, accounting for nonlinear gradient fields
JP5637694B2 (ja) 磁気共鳴イメージング装置及び非直交座標系走査法
JP5064685B2 (ja) 磁気共鳴イメージング装置
JP6579908B2 (ja) 磁気共鳴イメージング装置及び拡散強調画像計算方法
JP5484001B2 (ja) 磁気共鳴イメージング装置及び画像補正方法
US11914017B2 (en) Magnetic resonance imaging apparatus and image processing method
JP5638324B2 (ja) 磁気共鳴イメージング装置及び画像補正方法
JP4318837B2 (ja) 磁気共鳴イメージング装置
JP2018186892A (ja) 磁気共鳴イメージング装置
JP2018186942A (ja) 磁気共鳴イメージング装置及びパルス設計方法

Legal Events

Date Code Title Description
WWE Wipo information: entry into national phase

Ref document number: 201480022131.8

Country of ref document: CN

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

Ref document number: 14797920

Country of ref document: EP

Kind code of ref document: A1

ENP Entry into the national phase

Ref document number: 2015517044

Country of ref document: JP

Kind code of ref document: A

WWE Wipo information: entry into national phase

Ref document number: 14785454

Country of ref document: US

NENP Non-entry into the national phase

Ref country code: DE

122 Ep: pct application non-entry in european phase

Ref document number: 14797920

Country of ref document: EP

Kind code of ref document: A1