WO2013047583A1 - 画像解析装置、画像解析方法及び画像解析プログラム - Google Patents

画像解析装置、画像解析方法及び画像解析プログラム Download PDF

Info

Publication number
WO2013047583A1
WO2013047583A1 PCT/JP2012/074701 JP2012074701W WO2013047583A1 WO 2013047583 A1 WO2013047583 A1 WO 2013047583A1 JP 2012074701 W JP2012074701 W JP 2012074701W WO 2013047583 A1 WO2013047583 A1 WO 2013047583A1
Authority
WO
WIPO (PCT)
Prior art keywords
phase difference
image
distribution
image analysis
difference distribution
Prior art date
Application number
PCT/JP2012/074701
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 KR1020147011275A priority Critical patent/KR101685377B1/ko
Priority to JP2013536329A priority patent/JP6041356B2/ja
Priority to CN201280046660.2A priority patent/CN103826535B/zh
Priority to AU2012317518A priority patent/AU2012317518B2/en
Priority to US14/347,703 priority patent/US9330457B2/en
Priority to EP12836937.8A priority patent/EP2762071B1/en
Publication of WO2013047583A1 publication Critical patent/WO2013047583A1/ja

Links

Images

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/0002Inspection of images, e.g. flaw detection
    • G06T7/0012Biomedical image inspection
    • 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
    • AHUMAN NECESSITIES
    • A61MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61BDIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00Measuring for diagnostic purposes; Identification of persons
    • A61B5/40Detecting, measuring or recording for evaluating the nervous system
    • A61B5/4076Diagnosing or monitoring particular conditions of the nervous system
    • A61B5/4088Diagnosing of monitoring cognitive diseases, e.g. Alzheimer, prion diseases or dementia
    • 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/54Signal processing systems, e.g. using pulse sequences ; Generation or control of pulse sequences; Operator console
    • G01R33/56Image enhancement or correction, e.g. subtraction or averaging techniques, e.g. improvement of signal-to-noise ratio and resolution
    • G01R33/5608Data processing and visualization specially adapted for MR, e.g. for feature analysis and pattern recognition on the basis of measured MR data, segmentation of measured MR data, edge contour detection on the basis of measured MR data, for enhancing measured MR data in terms of signal-to-noise ratio by means of noise filtering or apodization, for enhancing measured MR data in terms of resolution by means for deblurring, windowing, zero filling, or generation of gray-scaled images, colour-coded images or images displaying vectors instead of pixels
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01RMEASURING ELECTRIC VARIABLES; MEASURING MAGNETIC VARIABLES
    • G01R33/00Arrangements or instruments for measuring magnetic variables
    • G01R33/20Arrangements or instruments for measuring magnetic variables involving magnetic resonance
    • G01R33/24Arrangements or instruments for measuring magnetic variables involving magnetic resonance for measuring direction or magnitude of magnetic fields or magnetic flux
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/30Subject of image; Context of image processing
    • G06T2207/30004Biomedical image processing

Definitions

  • the present invention relates to an image analysis apparatus, an image analysis method, and an image analysis program, and more particularly, to an image analysis apparatus, an image analysis method, and an image analysis program for analyzing a nuclear magnetic resonance image obtained from a living body.
  • Alzheimer's disease one of the dementias, is currently a major problem, and research is progressing all over the world.
  • One of the main causes of Alzheimer's disease is known to be a protein called amyloid ⁇ that accumulates in the brain.
  • the diameter of amyloid ⁇ is about 0.1 mm or less, and it is considered that direct depiction is impossible even with any current insurance medical image detection device.
  • ultra-high magnetic field MRI for example, 7 Tesla MRI
  • diagnostic imaging methods that require expensive equipment that is out of the scope of insurance medical treatment are not solutions to the current medical problem, but are realistic within the scope of current insurance medical treatment.
  • a method for detecting tissues and lesions that cannot be rendered by medical imaging equipment including
  • the MR system described in Patent Document 1 captures a high-resolution MR image at a magnetic field strength of 3 Tesla or more using a three-dimensional (3D) gradient double echo pulse sequence, and 3D MR image data is acquired at two different echo times so that a phase image is acquired.
  • the local variation of the magnetic field in the region of interest is obtained by subtracting the smoothed spherical harmonics from the measured values collected during MRI to estimate the constant and linear components related to the magnetic field inhomogeneity. Since the ability to measure magnetic field fluctuations in the brain by MRI increases according to the square of B0, according to the technique according to Patent Document 1, the measurement of brain magnetism and iron becomes more sensitive.
  • Patent Document 1 is intended to improve the overall sensitivity in MRI, and is intended only for direct rendering of MR images.
  • MRI that does not employ the technique described in Patent Document 1 described above cannot receive the benefit of improved sensitivity.
  • the present invention has been made in view of the above-described problems, and is an image that can realistically detect tissues and lesions that could not be rendered by a medical imaging device including MRI within the scope of current insurance medical treatment.
  • An object is to provide an analysis apparatus, an image analysis method, and an image analysis program.
  • One aspect of the image analysis device is an image analysis device for analyzing a nuclear magnetic resonance image obtained from a living body, and includes a phase difference of a nuclear magnetic resonance signal obtained from a predetermined region of the living body or A phase difference distribution creating unit that creates a distribution of phases (hereinafter collectively referred to as a phase difference), a fitting unit that fits a plurality of functions to the phase difference distribution, and a phase difference distribution A verification unit that verifies the normality of the living body included in the predetermined region based on the magnetic susceptibility of the tissue included in the predetermined region specified by the parameters of the plurality of fitted function groups. .
  • the above-described image analysis device includes various modes such as being implemented in another device or implemented together with another method.
  • the present technology provides an image analysis system including the image analysis device, an image analysis method having a process corresponding to the configuration of the above-described device, an image analysis program for causing a computer to realize a function corresponding to the configuration of the above-described device, and the image It can also be realized as a computer-readable recording medium on which an analysis program is recorded.
  • minute tissue changes can be detected. That is, by capturing a phase change derived from a partial volume effect of a tissue (partial volume effect) as a phase difference distribution change of a predetermined region of interest (region of interest in an embodiment described later), the imaging time is short, Even if MRI with a low magnetic field strength is used, the presence of the target tissue can be confirmed stably.
  • This method which can be detected using short-time imaging and low magnetic field devices, is not used for direct and highly accurate pixel-by-pixel image creation that takes a long imaging time and is created by an unrealistic method. It is thought that it is easy to be accepted in the medical field. Furthermore, since this method is processing by software, it can be used without modifying conventional devices or introducing additional devices.
  • FIG. 1 is a diagram showing a schematic configuration of an MRI apparatus 1 (magnetic resonance imaging apparatus) according to the present embodiment.
  • the MRI apparatus 1 is an apparatus that images internal information in the subject 2 using an NMR phenomenon.
  • the MRI apparatus 1 visualizes the rotation angle of the magnetization vector in addition to the intensity image obtained by imaging the intensity component of the nuclear magnetic resonance signal (MR signal) as the nuclear magnetic resonance image (MR image).
  • MR signal nuclear magnetic resonance signal
  • MR image nuclear magnetic resonance image
  • This is a new type of MRI apparatus that draws a morphological image using a phase image.
  • the MRI apparatus 1 and any one of the control system 20 and the information processing apparatus 26 described later constitute an image analysis apparatus.
  • the MRI apparatus 1 shown in FIG. 1 includes a coil system 10 and a control system 20.
  • the coil system 10 includes, for example, a static magnetic field coil unit 11, a gradient magnetic field coil unit 12, and an RF (Radio Frequency) coil unit 13. These have, for example, a substantially cylindrical shape, and are arranged so that their central axes (not shown) are coaxial with each other.
  • a bed portion 30 that supports the subject 2 is provided in a plane including the central axis.
  • the couch 30 is installed in the bore 10A (internal space) of the coil system 10.
  • the subject 2 on the bed unit 30 is carried into or out of the bore 10A by the movement of the bed unit 30 by a conveying means (not shown).
  • a conveying means not shown.
  • the direction parallel to the central axis is taken as the Z axis
  • the two directions perpendicular to the Z axis and perpendicular to each other are taken as the X axis and the Y axis.
  • the static magnetic field coil unit 11 forms a static magnetic field in the bore 10A.
  • the static magnetic field coil unit 11 is constituted by, for example, a superconducting coil or a normal conducting coil.
  • the direction of the static magnetic field formed by the static magnetic field coil unit 11 is substantially parallel to the Z-axis direction. 1 illustrates the case where the body axis direction of the subject 2 is parallel to the direction of the static magnetic field, but may be a direction orthogonal to the direction of the static magnetic field.
  • the gradient magnetic field coil unit 12 forms, for example, gradient magnetic fields (gradient magnetic fields) in directions of three axes perpendicular to each other, that is, a slice axis, a phase axis, and a frequency axis.
  • the gradient magnetic field coil unit 12 includes, for example, three types of coils: a slice axis direction coil, a phase axis direction coil, and a frequency axis direction coil.
  • any of the X, Y, and Z axes can be a slice axis.
  • the Z axis is a slice axis
  • the X axis can be a phase axis
  • the Y axis can be a frequency axis.
  • MR signals can be collected in coordinate systems other than the Cartesian coordinate system assumed above (for example, a polar coordinate system).
  • Cartesian coordinate system assumed above
  • An axis suitable for the coordinate system for example, radial direction, angular direction
  • the RF coil unit 13 forms an RF magnetic field that excites spins in the subject 2 in the space of the static magnetic field, and receives MR signals generated with the spins excited by the RF magnetic field.
  • the coil that receives the MR signal may be used also as a coil that forms the RF magnetic field, or may be provided separately from the coil that forms the RF magnetic field.
  • the RF coil unit 13 may be composed of a single coil, or may be composed of a plurality of coils 13-1 to 13-8 (multi-channel) as shown in FIG. Good.
  • an MR signal is obtained for each channel.
  • the MR signal is obtained by, for example, a GE (gradient echo) pulse sequence, and becomes a sampling signal in the frequency domain, that is, Fourier space (k-space). Yes.
  • the GE system includes, for example, a steady state in addition to the GE.
  • the pulse sequence may be, for example, BalancedBSSFP (Steady State Free Precession), TrueFISP (True Fast Imaging with Steady-state Precession), and other than GE such as SE (Spin Echo). Also good.
  • control system 20 includes a static magnetic field power source 21, a gradient magnetic field power source 22, a transmission unit 23, a reception unit 24, and a sequence control unit 25.
  • the static magnetic field power source 21 supplies power to the static magnetic field coil unit 11 to drive the coil system 10. When this electric power is supplied to the static magnetic field coil unit 11, a static magnetic field is formed in the bore 10A.
  • the gradient magnetic field power supply 22 supplies power to the gradient magnetic field coil unit 12 based on a control signal input from the sequence control unit 25. By supplying power to the gradient coil unit 12, desired gradient magnetic fields (gradient magnetic fields) are formed in the respective directions of the slice axis, the phase axis, and the frequency axis.
  • the transmission unit 23 applies the RF signal to the RF coil unit 13 based on the control signal input from the sequence control unit 25, for example.
  • the receiving unit 24 receives MR signals generated by driving the coil system 10. For example, MR signals received by the RF coil unit 13 are detected, and necessary signal processing is performed, and A / D (Analog-to-Digital) conversion is performed, thereby digitizing complex data (raw data). Is generated. Of course, the receiving unit 24 may directly perform A / D conversion on the detected MR signal to generate raw data. The raw data generated by the receiving unit 24 is output to the sequence control unit 25, for example.
  • the sequence control unit 25 drives the gradient magnetic field power source 22, the transmission unit 23, and the reception unit 24 for driving the MRI apparatus 1.
  • the sequence control unit 25 drives the gradient magnetic field power source 22, the transmission unit 23, and the reception unit 24 by applying control signals to the gradient magnetic field power source 22, the transmission unit 23, and the reception unit 24.
  • the control signal is generated in accordance with a pulse sequence that defines the magnitude of the pulse current applied to the gradient magnetic field power source 22, the transmitter 23, and the receiver 24, the application time, the application timing, and the like.
  • Information about this pulse sequence is input to the sequence control unit 25 from an information processing device 26 described later.
  • the sequence control unit 25 outputs the raw data input from the reception unit 24 to the information processing device 26, for example.
  • the control system 20 further includes an information processing device 26 as shown in FIG. 1, for example.
  • the information processing device 26 includes, for example, a calculation unit 26A, an input unit 26B, a display unit 26C, and a storage unit 26D.
  • the input unit 26B is, for example, a device that captures information from the user as digital data into the information processing device 26, and includes, for example, a keyboard, a mouse, and a scanner.
  • the display unit 26C is for displaying a result (for example, morphological image) processed by the computing unit 26A, a dialog for data input such as an imaging condition, and the like, and is configured by a display device such as a liquid crystal display device. ing.
  • Various programs for controlling the MRI apparatus 1 are stored in the storage unit 26D.
  • a program 27 used for determination processing described later, a phase difference enhanced imaging program used for creating morphological images described later, Etc. are stored.
  • the computing unit 26A is, for example, for interpreting and executing instructions of various programs, and is configured by, for example, a CPU (Central Processing Unit).
  • a program 27 stored in the storage unit 26D is loaded on the calculation unit 26A in conjunction with the activation of the MRI apparatus 1, so that the calculation unit 26A, for example, in response to an instruction from the user, the program 27 The command is interpreted and executed.
  • the computing unit 26A may be configured by hardware corresponding to the functions of various programs (for example, the program 27).
  • various programs for example, the program 27.
  • a phase image is created by executing various instructions of the program 27 in the arithmetic unit 26A.
  • phase difference image creation of a phase difference image performed by the MRI apparatus 1 will be described.
  • an intensity image and a phase image are required.
  • another pulse sequence included in the GE system or a pulse sequence other than the GE system may be used.
  • the phase image obtained by the GE pulse sequence is the product of the amount of change ⁇ B of the local magnetic field (local magnetic field compared to the external magnetic field) created by the tissue included in each pixel and the echo time (TE) required for imaging. It is proportional to ( ⁇ B ⁇ TE). Therefore, in order to extract large phase (difference) information from the phase image, TE is increased or a function for enhancing ⁇ B (so-called enhancement function) is changed to a stronger one.
  • imaging is performed using a predetermined pulse sequence.
  • the number of times of imaging may be one or may be a plurality of times so that it can be statistically handled.
  • the RF coil unit 13 is composed of multiple channels and MR signals are obtained for each channel, the arithmetic average of the MR signals obtained for each channel is calculated, and the arithmetically averaged MR signal is used.
  • sensitivity correction may be performed on individual MR signals in advance.
  • FIG. 3 is a diagram showing a flow of data processing until a phase difference image is created.
  • the data processing shown in the figure is executed in the control system 20 under the control of the arithmetic unit 26A.
  • the calculation unit 26A starts calculation of data processing shown in FIG.
  • the computing unit 26A first outputs a control signal for requesting the sequence control unit 25 to acquire raw data using a predetermined pulse sequence. Then, a control signal according to a predetermined pulse sequence is output from the sequence control unit 25 to the gradient magnetic field power source 22, the transmission unit 23, and the reception unit 24.
  • the RF coil unit 13 detects the MR signal.
  • the MR signal detected here is converted into raw data R by predetermined signal processing in the receiving unit 24.
  • the receiving unit 24 inputs the raw data R to the sequence control unit 25, and the sequence control unit 25 transfers (inputs) the raw data R to the arithmetic unit 26A. In this way, the calculation unit 26A acquires data (raw data R) corresponding to the MR signal.
  • the calculation unit 26A arranges the raw data R input from the sequence control unit 25 in the k space set in an internal memory (not shown).
  • data arranged in the k space is referred to as k space data S (k).
  • the computing unit 26A performs inverse Fourier transform on the k-space data S (k) to reconstruct a complex image ⁇ (x).
  • the complex image ⁇ (x) is an image having a real image in the real part and an imaginary image in the imaginary part as shown by the following formula (1).
  • the computing unit 26A obtains an intensity image M (x) and a phase image P (x) from the complex image ⁇ (x).
  • phase wrapping occurs in the phase image P (x), and the phase exceeding 2 ⁇ takes a phase value obtained by subtracting 2 ⁇ n (n is an integer) from the actual phase. Therefore, the phase image P (x) becomes a striped pattern image and does not show the original phase value. Therefore, the arithmetic unit 26A performs processing for removing the phase wrapping and extracting the phase difference.
  • the computing unit 26A first performs a Fourier transform on the complex image ⁇ (x) to return the complex image ⁇ (x) once to the k-space data S (k). Alternatively, the k-space data S (k) arranged in the k-space is read out. Next, the computing unit 26A applies LPF (Low Pass Filter) to the k-space data S (k), and performs inverse Fourier on the data L (k) ⁇ S (k) obtained thereby. The transformation is performed to obtain a complex image ⁇ ′ (x). Note that L (k) is a function of LPF.
  • the computing unit 26A creates a phase difference image PD (x) using the complex images ⁇ (x) and ⁇ ′ (x). Specifically, the calculation unit 26A divides the complex image ⁇ (x) by the complex image ⁇ ′ (x), and performs a complex quotient calculation to create a phase difference image PD (x). This eliminates phase wrapping in the phase portion.
  • the phase difference included in the phase difference image PD (x) has a width of 2 ⁇ .
  • the phase difference included in the phase difference image PD (x) is ⁇ ⁇ PD ( x) Assuming ⁇ .
  • the sign of the phase difference included in the phase difference image PD (x) is determined by ⁇ ⁇ ⁇ B ⁇ TE.
  • phase difference image PD (x) is included even if the magnitude of the phase difference does not change.
  • the sign of the phase difference may change. Therefore, in order to cope with this, the above definition is made so that the sign is negative particularly for venous blood.
  • said (gamma) is a positive proportionality constant, for example, corresponds to the gyromagnetic ratio of hydrogen.
  • a region in which the user of the MRI apparatus 1 is interested (hereinafter referred to as a region of interest) is set for the phase difference image, and a predetermined calculation is performed on the phase difference image included in the region of interest. Process. By this calculation process, the normality (or non-normality) of the tissue included in the region of interest can be determined.
  • the data processing shown in the figure is executed in the control system 20 under the control of the arithmetic unit 26A.
  • the calculation unit 26A starts the data processing calculation shown in FIG.
  • Display area specification interface When the process is started, an interface for designating a section or a part of the space imaged by the MRI apparatus 1 is displayed on the screen of the display unit 26C (S1).
  • this interface for example, an intensity image or a phase image created based on the MR signal is displayed, and a predetermined area on the image can be designated by area designation means such as a closed curve (frame line, etc.) or coordinate values. It is like that.
  • the operator of the MRI apparatus 1 designates a region of a cross section or a part of the space imaged by the MRI apparatus 1 via the interface displayed on the display unit 26C (S2).
  • the region designated by the operator is the region of interest.
  • the operator designates, as a region of interest, a region that is considered to contain an abnormal tissue.
  • the region of interest may be a two-dimensional region or a three-dimensional region.
  • the calculation unit 26A acquires all phase data of the MR signal acquired from the tissue included in the region of interest, statistically analyzes these phase data, for example, uses the horizontal axis as the phase value, A phase difference distribution with the axis as the number of data is created (S3).
  • a function group used for fitting the phase difference distribution is selected (S4).
  • the operator selects an appropriate function group according to the target tissue, lesion, disease state, imaging method, and the like.
  • step S4 may be skipped and a predetermined function group corresponding to the purpose may be automatically applied.
  • a predetermined function group corresponding to the target tissue, lesion, disease state, imaging method, etc. is prepared, and an appropriate function group corresponding to the purpose is selected by selecting the purpose from the selection screen. It may be.
  • a single phase difference distribution created for one region of interest is simultaneously fitted by a plurality of functions (S5). That is, one phase difference distribution is approximated by a curve obtained by superposing a plurality of functions.
  • Various functions can be adopted as the plurality of functions, and examples include a Gaussian distribution, a Lorentz distribution, and a binomial distribution.
  • the plurality of functions need not be limited to distribution functions, and need not be orthogonal to each other.
  • the plurality of functions may be configured by combining different types of functions. Note that when performing fitting by computer computation, the plurality of functions must be finite.
  • At least one of the plurality of functions used for fitting uses a Gaussian distribution. This is because the signal distribution obtained by random organization is a random variable expressed as the sum of many independent factors, and is guaranteed to follow a Gaussian distribution by the central limit theorem.
  • each function used for fitting has one or more parameters for changing the function shape, when performing fitting using a plurality of functions, at least parameters equal to or more than the number of functions used for fitting are adjusted. There is a need.
  • a parameter set in which the phase difference distribution and the superposition of a plurality of functions are closest is searched while appropriately changing a plurality of parameters.
  • the fitting function is a Gaussian distribution
  • the number of parameters is as small as three, and the adjustment of the parameters is easy.
  • the degree of approximation between the plurality of functions and the phase difference distribution can be evaluated by, for example, the least square method.
  • a double Gaussian distribution model using two Gaussian distributions as a plurality of functions is adopted. Since the Gaussian distribution has three parameters: height, center position (average), and standard deviation ⁇ (or variance ⁇ ⁇ 2), in the double Gaussian distribution model, fitting is performed while adjusting the six parameters. become.
  • Equation (2) is a function used for fitting in the double Gaussian distribution model.
  • A1 corresponds to the height of the first Gaussian distribution
  • B1 corresponds to the reciprocal of the variance of the first Gaussian distribution
  • C1 corresponds to the center position of the first Gaussian distribution.
  • A2 corresponds to the height of the second Gaussian distribution
  • B2 corresponds to the reciprocal of the variance of the second Gaussian distribution
  • C2 corresponds to the center position of the second Gaussian distribution.
  • the parameter set obtained by the fitting described above is a combination of values characterizing the magnetic susceptibility of the tissue included in the region of interest. That is, the parameter set when the region of interest contains a non-normal tissue having a different magnetic susceptibility from the normal tissue is different from the parameter set when the region of interest contains only normal tissue.
  • the parameter set obtained by the fitting is displayed on the display unit 26C (S6), or based on the parameter set, the divergence degree of the tissue included in the region of interest from the normal tissue (the normality of the tissue (or non-normal) Normality)) is calculated, and this degree is displayed on the display unit 26C (S7).
  • the operator can obtain a reference for determining whether the tissue included in the region of interest is normal or abnormal, and can obtain a measure of the degree of deviation from the normal tissue.
  • imaging was performed using 3D-FLASH with a 7T-MRI apparatus, and an intensity image and a phase image were simultaneously acquired.
  • the main imaging conditions are TR / TE: 50 / 12.8 ms, Flip Angle: 20 °, Matrix Size: 194 ⁇ 128 ⁇ 82 (0.08 mm isovoxel), addition: 24 times.
  • intensity images and phase images obtained from the cortex, hippocampus and thalamus are used in the brain.
  • the cortex, hippocampus, and thalamus accumulate particularly large amounts of senile plaques, which is one of Alzheimer's pathological changes.
  • iron deposits in the senile plaque and it is considered that this change in magnetic susceptibility due to iron can be captured as a phase signal.
  • the region of interest is set so as to include the part considered to have senile plaques, and the phase difference distribution of each magnetic resonance signal collected from these brain parenchyma is converted into a double Gaussian distribution model. And fitting.
  • FIG. 5 is a graph plotting the phase acquired from the region of interest set in the control mouse
  • FIGS. 6 and 7 are graphs plotting the phase acquired from the region of interest set in the genetically modified mouse.
  • the horizontal axis indicates the phase (radian)
  • the vertical axis indicates the number of detections (for example, the number of pixels).
  • the phase value is adjusted to the range of ⁇ to ⁇ after removing the phase wrapping, but FIGS. 5 to 7 enlarge the vicinity of the center of the phase difference distribution (phase 0). It is shown.
  • the first Gaussian distribution is 7 to 8 times as high as the second Gaussian distribution.
  • the first Gaussian distribution and the second Gaussian distribution have substantially the same height.
  • the first Gaussian distribution is 8 to 9 times as high as the second Gaussian distribution. That is, it can be seen that the contribution rate (height) of each Gaussian distribution is different in each case even for the same lesion, and cannot be adopted as a criterion for the presence or absence of an abnormal tissue.
  • both the first Gaussian distribution and the second Gaussian distribution show sharp distributions, and in particular, the first Gaussian distribution compared to the first Gaussian distribution.
  • the Gaussian distribution of 2 has a narrow width (for example, half width), and any Gaussian distribution has a narrow width (for example, half width) compared to the phase difference distribution.
  • the width of the second Gaussian distribution (for example, the half-value width) is wider than the first Gaussian distribution, and the second Gaussian distribution is wider than the phase difference distribution. (For example, half width) is wide.
  • the standard deviation (width) of the Gaussian distribution shows a clear difference between the normal tissue of the control mouse and the tissue containing senile plaques of the genetically modified mouse, and can be a criterion for the presence or absence of an abnormal tissue. I understand.
  • the first Gaussian distribution constituting the fitting curves shown in FIGS. 5 to 7 is caused by noise components.
  • the MRI apparatus 1 it is naturally predicted that an error occurs in the phase difference due to electrical thermal noise, digitalized noise caused by cutoff, and the like. This is because these noises are known to follow a Gaussian distribution, and it is considered that a relatively high Gaussian distribution is generated. Therefore, when fitting the phase difference distribution with a plurality of functions, it is considered that a good fitting result can be obtained if a Gaussian distribution is adopted for at least one of the plurality of functions.
  • the second Gaussian distribution appearing in the fitting curves shown in FIGS. 6 and 7 is considered to be related to a tissue containing some magnetic substance.
  • the phase difference distribution shown in FIGS. 6 and 7 is acquired from a region of interest set to include senile plaques. This is because senile plaque is a tissue containing iron having magnetism as described above, and has a relatively large standard deviation when the phase difference distribution of the tissue containing a magnetic substance is fitted with a Gaussian distribution.
  • the phase difference distribution shown in FIG. 5 has a sharp rise relatively close to the Gaussian distribution, but the phase difference shown in FIGS.
  • the distribution is wider than the Gaussian distribution of the same height, suggesting that some kind of broad phase difference distribution is mixed. That is, when the bottom of the phase difference distribution is wider than the Gaussian distribution having the same height, it indicates that the region of interest includes a tissue containing some magnetic substance.
  • the center shift between the phase difference distribution and the second Gaussian distribution is larger in FIGS. 6 and 7 than in FIG.
  • the reason for this is that when normal tissue is considered to have a phase difference of 0, a Gaussian distribution having a center close to this phase difference of 0 is considered to be a Gaussian distribution indicated by normal tissue, and has a center away from this phase difference of 0.
  • the Gaussian distribution is considered to be a tissue that has a high possibility of containing many tissues that are not normal tissues. That is, the distance between the center of the Gaussian distribution corresponding to the normal tissue and the center of the Gaussian distribution corresponding to the abnormal tissue is considered to be a standard indicating the normality (or non-normality) of the tissue included in the region of interest. It is done.
  • phase difference enhancement imaging method for enhancing a specific phase difference, which is performed using a fitting function obtained using a double Gaussian distribution model, will be described.
  • the phase difference enhancement imaging method is selected by selecting a part of the obtained phase difference image PD (x), selecting a part or all of the phase difference image PD (x), and emphasizing by the enhancement function w ( ⁇ ).
  • This is a method of expressing a part of an image corresponding to the phase information on the intensity image M (x).
  • an image created by the enhancement function w ( ⁇ ) can be provided.
  • FIG. 8 shows the relationship between the fitting function obtained by the double Gaussian distribution model and the phase difference to be emphasized.
  • FIG. 8 shows the phase range to be emphasized by taking the experimental result of FIG. 6 as an example.
  • the phase where the second Gaussian distribution intersects the first Gaussian distribution is defined as ⁇ 1
  • the phase range in which the second Gaussian distribution is larger than the first Gaussian distribution ( ⁇ > ⁇ 1) is represented as the enhancement range. It is as.
  • the second Gaussian distribution has a larger standard deviation than the first Gaussian distribution, and has a wide base. Have. Therefore, the visibility of the non-normal tissue is improved by emphasizing the tissue in which the second Gaussian distribution indicating the distribution of the non-normal tissue has a phase exceeding the first Gaussian distribution indicating the distribution of the normal tissue. It is considered that morphological images can be provided.
  • FIG. 10 is a diagram showing a flow of data processing until a morphological image is created.
  • the computing unit 26A selects a phase where the phase ⁇ of the phase difference image PD (x) is larger than ⁇ 1, and emphasizes the selected phase ⁇ .
  • the phase difference distribution changes, for example, as in the models shown in FIGS. 9A to 9C by changing the LPF size. That is, as the LPF size increases, the width of the substantially symmetrical phase difference distribution centering on the phase difference zero decreases. For this reason, the phase difference image mainly has a phase difference of zero and values in the vicinity thereof, and it is difficult to add tissue contrast as the phase difference image. On the other hand, in a tissue having a fine structure, there is an aspect that it is easier to add contrast when using a large LPF than when using a small LPF.
  • an enhancement function that can be flexibly dealt with is selected as described later, and the same as shown in FIGS. 9A to 9C.
  • a device is devised so that the same contrast can be given to those showing different frequency distributions with respect to the phase difference.
  • the calculation unit 26A selects the width of the phase difference and the center value thereof so that the target tissue contrast has a desired magnitude. Note that the selection of the width of the phase difference and its center value may be left to the user of the MRI apparatus 1 instead of relying on the calculation unit 26A. In this case, the calculation unit 26A selects the target tissue phase ⁇ in consideration of the change in the phase difference distribution due to the filter processing.
  • the computing unit 26A obtains an enhanced image w (PD (x)) by enhancing the selected phase ⁇ with the enhancement function w ( ⁇ ).
  • an exponential function is used as the enhancement function w ( ⁇ ).
  • a ⁇ function is used as an example of an exponential function. This ⁇ function is expressed by the following two expressions.
  • the parameters a, b, and ⁇ all take real values.
  • the parameters a and b adjust the degree of phase difference enhancement, and are determined by the LPF filter size.
  • the parameters a and b are also determined so as to maximize the contrast C or contrast-to-noise ratio CNR between the target tissue and its background.
  • the parameter ⁇ reduces noise on the phase difference image PD (x), and is determined by, for example, the standard deviation of the tissue in which the phase average value on the phase difference image PD (x) takes a value near 0 (zero).
  • the parameter ⁇ can be obtained based on data obtained from many experiments.
  • the parameter ⁇ is determined by, for example, the contrast C or the contrast / noise ratio CNR.
  • the contrast C described above is obtained by emphasizing the signal w (PD (x1)) ⁇ M (x1) of the target tissue at the position x1 on the image and the image on the image, as shown in the following equation. It is determined by the absolute value of the difference from the background signal w (PD (x2)) ⁇ M (x2) of the emphasized tissue at position x2.
  • contrast / noise ratio CNR is expressed by C / ⁇ ′ as shown by the following equation.
  • ⁇ ′ is determined by the standard deviation on the emphasized image of the target tissue or the standard deviation on the emphasized image of the target tissue background.
  • any of the standard deviations described above may not be adopted.
  • the standard deviation of the external signal of the subject 2 or the standard deviation by the difference method is adopted.
  • Abs ( ⁇ ) represents the absolute value of ⁇ .
  • the ⁇ function is a function that does not emphasize PD (x) in the range of Abs ( ⁇ ) ⁇ ⁇ and emphasizes the phase difference image PD (x) in other ranges. Since this ⁇ function can approximate an arbitrary power function with arbitrary accuracy, more flexible emphasis can be performed as compared with the case where a polynomial is used as the emphasis function.
  • the filter size of the LPF when the filter size of the LPF is changed in accordance with the size of the imaging region, the distribution of the phase difference slightly changes accordingly.
  • the filter size of the LPF is changed in a state where the size of the imaging area is fixed, the distribution of the phase difference greatly changes accordingly.
  • the parameters a, b, and ⁇ of the enhancement function take real values, and these can be flexibly changed according to the imaging conditions and the like. Thereby, even when the imaging conditions are changed, it is possible to maintain the same or the same contrast.
  • the computing unit 26A masks the intensity image M (x) with the enhanced image w (PD (x)), for example, according to a predetermined mode (rule), thereby obtaining the morphological image I (x).
  • a predetermined mode for masking the intensity image M (x) with the emphasized image w (PD (x))
  • Specific conditions for masking the intensity image M (x) with the emphasized image w (PD (x)) can be set according to the object to be emphasized, and are basically exemplified below.
  • the four types of enhancement modes are set.
  • the total enhancement does not depend on the parameter ⁇ and the sign of the phase difference image PD (x), but in the case of structure enhancement, for example, the phase difference ⁇ (hereinafter simply referred to as the phase difference ⁇ ) generated by the intracortical structure.
  • the phase difference ⁇ (hereinafter simply referred to as the phase difference ⁇ ) generated by the intracortical structure.
  • a value is obtained in advance by an experiment, and a conditional expression is set according to the magnitude relationship between the phase difference ⁇ and the phase difference image PD (x).
  • Enhancement mode B blood vessel enhancement
  • I (x) w (PD (x)) ⁇ M (x) (PD (x) ⁇ 0)
  • I (x) M (x) (PD (x)> 0)
  • the calculation unit 26A extracts a portion where the phase difference image PD (x) is 0 (zero) or more, and creates a morphological image I (x) by performing enhancement only on that portion. At this time, no emphasis is performed on a portion where the phase difference image PD (x) is negative.
  • the morphological image I (x) is a tissue-emphasized image.
  • the computing unit 26A extracts a portion where the phase difference image PD (x) is 0 (zero) or less, and performs enhancement only on that portion, thereby creating a morphological image I (x). At this time, no emphasis is performed on a portion where the phase difference image PD (x) is positive.
  • the morphological image I (x) is a blood vessel enhanced image.
  • the calculation unit 26A creates the morphological image I (x) by enhancing the entire phase difference image PD (x).
  • the morphological image I (x) is an image in which the whole including tissues and blood vessels is enhanced.
  • the calculation unit 26A enhances a portion that satisfies ⁇
  • a morphological image I (x) is created. At this time, the portion satisfying PD (x) ⁇
  • the morphological image I (x) is a structure-enhanced image. Since the phase difference ⁇ is generated by the intracortical structure, the morphological image I (x) is actually an image in which the cortical structure is emphasized.
  • each morphological image I it is possible to specify a more accurate anatomical position of the brain function against the background of the tissue contrast indicated by x).
  • an MR signal phase difference distribution obtained from a region of interest is created, and this phase difference distribution is simultaneously fitted with a plurality of function groups.
  • the normality of the living body included in the region of interest is verified based on the magnetic susceptibility of the tissue included in the region of interest specified by the parameters of the plurality of function groups fitted to the phase difference distribution. Therefore, it becomes possible to realistically detect tissues and lesions that could not be drawn with medical imaging equipment including MRI within the scope of current insurance medical treatment.
  • This technology can be used for various diagnoses using MRI. In particular, it can greatly contribute to the detection of minute lesion tissues such as amyloid ⁇ .
  • This method can be applied in the same way to detect diffusely present in tissues, for example, quantitatively measuring the proportion of fat in the liver that is always examined in health examinations, etc. It is considered possible.

Landscapes

  • Health & Medical Sciences (AREA)
  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • General Health & Medical Sciences (AREA)
  • Radiology & Medical Imaging (AREA)
  • High Energy & Nuclear Physics (AREA)
  • Medical Informatics (AREA)
  • Neurology (AREA)
  • General Physics & Mathematics (AREA)
  • Heart & Thoracic Surgery (AREA)
  • Veterinary Medicine (AREA)
  • Biophysics (AREA)
  • Pathology (AREA)
  • Biomedical Technology (AREA)
  • Molecular Biology (AREA)
  • Surgery (AREA)
  • Animal Behavior & Ethology (AREA)
  • Public Health (AREA)
  • Signal Processing (AREA)
  • Condensed Matter Physics & Semiconductors (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Child & Adolescent Psychology (AREA)
  • Developmental Disabilities (AREA)
  • Hospice & Palliative Care (AREA)
  • Artificial Intelligence (AREA)
  • Psychiatry (AREA)
  • Psychology (AREA)
  • Neurosurgery (AREA)
  • Physiology (AREA)
  • Quality & Reliability (AREA)
  • Theoretical Computer Science (AREA)
  • Magnetic Resonance Imaging Apparatus (AREA)

Abstract

MRIを含む医療用画像機器で描出できなかった組織や病変を、現在の保険診療の範囲内で現実的に検知する。生体から得られたMR画像を解析するにあたり、関心領域から得たMR信号の位相差分布を作成し、この位相差分布を複数の関数群で同時にフィッティングし、位相差分布に対してフィッティングされた前記複数の関数群のパラメータに基づいて、前記関心領域に含まれる組織の磁化率に対する依存状態を決定する。

Description

画像解析装置、画像解析方法及び画像解析プログラム
 本発明は画像解析装置、画像解析方法及び画像解析プログラムに関し、特に、生体から得られた核磁気共鳴画像を解析するための画像解析装置、画像解析方法及び画像解析プログラムに関する。
 痴呆症の一つであるアルツハイマー病は、現在大きな問題であり、世界中でこぞって研究が進んでいる。このアルツハイマー病を引き起こす主原因の一つは、脳内に蓄積するアミロイドβと呼ばれるたんぱくであることが知られている。しかし、アミロイドβの直径は、約0.1mm以下であり、現行のいかなる保険診療用画像検出機器をもってしても、直接の描出は不可能であると考えられる。
 ところが、MRIを含む医療用画像機器のほとんどは、対象とする組織や病変を直接に描出することを目的としている。しかし、医療現場には撮像時間の制限があり、制限された時間の中で撮像できる画像には解像度の限界がある。このような制限下では、対象とする組織や病変を直接に描出できない。むろん、高解像度の画像機器であれば描出が可能であるが、このような画像機器は保険診療の範囲外であり、機器も高額である。
 むろん、超高磁場MRI(例えば、7テスラのMRI等)を用いれば、直接描出の可能性もあるが、このような超高磁場MRIは保険診療の範囲外である。保険診療の範囲外となる高額な器具を必要とする画像診断法は、現在の火急ともいえる医療問題への解決策とは言えず、現在の保険診療の範囲内で現実的に、現在のMRIを含む医療用画像機器では描出できない組織や病変を、検出する方法が求められている。
 このような背景から、例えば、特許文献1に記載のMRシステムは、3次元(3D)傾斜2重エコー・パルスシーケンスを用いて3テスラ以上の磁場強度において高分解能のMR画像を撮像し、3D位相画像が収集されるように2つの異なるエコー時間でMR画像の画像データを収集している。
 そして、磁場の不均一性に関する一定成分と線形成分を推定するためにMRIの間に収集される計測値から平滑化球面調和関数を差し引いて、関心領域内の磁場の局所変動を得ている。MRIによって脳内の磁場変動を計測できる能力はB0の二乗に従って増加するため、特許文献1に係る技術によれば、脳の磁性及び鉄の計測がより高感度になる。
特開2005-28151号公報
 しかしながら、上述した特許文献1に記載の技術は、MRIにおける全体的な感度向上を目的としたものであり、あくまで、MR画像の直接の描出を目的としている。また、現在の保険診療の範囲内で、現在のMRIを含む医療用画像機器で描出できない組織や病変を検出できるか不明である。さらに、上述した特許文献1に記載の技術を採用していないMRIは、感度向上の恩恵を受けられない。
 本発明は、前記課題に鑑みてなされたものであり、MRIを含む医療用画像機器で描出できなかった組織や病変を、現在の保険診療の範囲内で現実的に検知することが可能な画像解析装置、画像解析方法及び画像解析プログラムの提供を目的とする。
 本技術に係る画像解析装置の態様の1つは、生体から得られた核磁気共鳴画像を解析するための画像解析装置であって、生体の所定領域から得た核磁気共鳴信号の位相差もしくは位相(以下、統一して位相差と呼ぶ。)の分布を作成する位相差分布作成部と、前記位相差分布に対し、複数の関数群でフィッティングするフィッティング部と、前記位相差分布に対してフィッティングされた前記複数の関数群のパラメータによって特定される前記所定領域に含まれる組織の磁化率に基づいて、前記所定領域に含まれる生体の正常性を検証する検証部と、を備える構成としてある。
 なお、上述した画像解析装置は、他の機器に組み込まれた状態で実施されたり他の方法とともに実施されたりする等の各種の態様を含む。また、本技術は前記画像解析装置を備える画像解析システム、上述した装置の構成に対応した工程を有する画像解析方法、上述した装置の構成に対応した機能をコンピュータに実現させる画像解析プログラム、該画像解析プログラムを記録したコンピュータ読み取り可能な記録媒体、等としても実現可能である。
 組織の微小な磁性変化に鋭敏なMRI位相画像を用いて、微小な組織変化を検出することができる。即ち、組織の部分容積効果(partial volume effect)に由来する位相変化を、対象とする所定領域(後述の実施形態においては、関心領域)の位相差分布変化としてとらえることにより、撮像時間が短く、磁場強度が低いMRIを用いても、安定して対象組織の存在を確認することができる。
 長い撮像時間をかけ非現実的な手法により作成する直接かつ確度の高いpixel by pixelの画像作成ではなく、短時間の撮像や低磁場の機器を用いても検出可能な本手法は、産業界や医療現場に受け入れられやすいと考えられる。さらに、本方法はソフトウェアによる処理であるため、従来の機器を改変したり、追加の機器を導入したりせずに、利用可能である。
MRI装置(磁気共鳴画像化装置)の概略構成を表した図である。 RFコイル部の変形例を示す断面構成図である。 位相差画像を作成するまでのデータ処理の流れを示した図である。 関心領域内の組織が正常か否か判断する処理の流れを示すフローチャートである。 正常組織における位相差分布をフィッティングしたグラフである。 非正常組織における位相差分布をフィッティングしたグラフである。 非正常組織における位相差分布をフィッティングしたグラフである。 二重ガウス分布モデルにより得られるフィッティング関数と、強調すべき位相差との関係を示す図である。 位相差分布とLPFサイズとの関係について説明するための模式図である。 形態画像を作成するまでのデータ処理の流れを示した図である。
 以下、下記の順序に従って本技術を説明する。
(1)本実施形態の構成:
(2)位相差画像の作成:
(3)判断処理:
(4)形態画像の作成:
(5)まとめ:
(1)本実施形態の構成: 
 図1は、本実施形態に係るMRI装置1(磁気共鳴画像化装置)の概略構成を表した図である。MRI装置1は、NMR現象を利用して被検体2内の内部情報を画像化する装置である。MRI装置1は、後述するように、核磁気共鳴画像(MR画像)として、核磁気共鳴信号(MR信号)の強度成分を画像化した強度画像の他に、磁化ベクトルの回転角を画像化した位相画像を利用して形態画像を描画する、新しいタイプのMRI装置である。なお、本実施形態においては、MRI装置1と後述の制御システム20や情報処理装置26のいずれかが、画像解析装置を構成する。
 図1に示すMRI装置1は、コイルシステム10、制御システム20を備えている。
[コイルシステム10] 
 コイルシステム10は、例えば、静磁場コイル部11、傾斜磁場コイル部12、RF(Radio Frequency)コイル部13を含んで構成されている。これらは、例えば、概ね円筒状の形状となっており、それぞれの中心軸(図示せず)が互いに同軸となるように配置されている。その中心軸を含む面内に、被検体2を支持する寝台部30が設けられている。
 寝台部30は、コイルシステム10のボア10A(内部空間)に設置されている。寝台部30上の被検体2は、図示しない搬送手段による寝台部30の移動により、ボア10Aに搬入されたり、搬出されたりするようになっている。なお、本実施形態では、図1に示すように、中心軸と平行な方向をZ軸とし、Z軸と直交し、かつ互いに直交する2方向をX軸、Y軸とする。
 静磁場コイル部11は、ボア10Aに静磁場を形成する。静磁場コイル部11は、例えば、超伝導コイルや、常伝導コイル等によって構成されている。静磁場コイル部11によって形成される静磁場の方向は、概ねZ軸方向と平行となっている。なお、図1では、被検体2の体軸方向が、静磁場の方向と平行となっている場合が例示されているが、静磁場の方向と直交する方向となっていてもよい。
 傾斜磁場コイル部12は、例えば、互いに垂直な3つの軸、即ち、スライス軸、位相軸および周波数軸のそれぞれの方向に傾斜磁場(勾配磁場)を形成するものである。この傾斜磁場コイル部12は、例えば、スライス軸方向用コイル、位相軸方向用コイル、周波数軸方向用コイルの3種類のコイルからなる。このとき、X軸、Y軸およびZ軸のいずれの軸についてもスライス軸とすることが可能である。例えば、Z軸をスライス軸とした場合には、X軸を位相軸とし、Y軸を周波数軸とすることが可能である。
 なお、MR信号の収集は上で想定されているカーテシアン座標系以外の座標系(例えば極座標系)でも行うことは可能であり、そのような座標系でMR信号の収集を行う場合には、その座標系に適した軸(例えば動径方向、角度方向)が設定される。
 RFコイル部13は、静磁場の空間に、被検体2内のスピンを励起するRF磁場を形成すると共に、RF磁場によって励起されたスピンに伴い発生するMR信号を受信する。なお、RFコイル部13において、MR信号を受信するコイルは、RF磁場を形成するコイルと兼用されてもよいし、RF磁場を形成するコイルと別体に設けられてもよい。また、RFコイル部13は、単一のコイルで構成されていてもよいし、例えば、図2に示したように、複数のコイル13-1~13-8(多チャンネル)で構成されてもよい。なお、RFコイル部13が多チャンネルで構成されている場合には、チャンネルごとにMR信号が得られる。
 前記のMR信号は、例えば、GE(gradient echo;グラディエントエコー)系のパルスシーケンスによって得られるものであり、周波数ドメイン、即ち、フーリエ空間(k空間(k-space))についてのサンプリング信号となっている。GE系には、GEの他に、例えば、ステディステート(steady state)が含まれる。なお、パルスシーケンスは、例えば、Balanced SSFP(Steady State Free Precession)、TrueFISP(True Fast Imaging with Steady-state Precession)であってもよく、SE(Spin Echo;スピンエコー)などのGE以外ものであってもよい。
[制御システム20] 
 制御システム20は、例えば、図1に示したように、静磁場電源21、傾斜磁場電源22、送信部23と、受信部24と、シーケンス制御部25とを備えている。
 静磁場電源21は、静磁場コイル部11に電力を供給し、コイルシステム10を駆動する。この電力が静磁場コイル部11へ供給されると、ボア10Aに静磁場が形成される。傾斜磁場電源22は、シーケンス制御部25から入力された制御信号に基づいて傾斜磁場コイル部12に電力を供給するものである。傾斜磁場コイル部12への電力の供給により、所望の傾斜磁場(勾配磁場)がスライス軸、位相軸および周波数軸のそれぞれの方向に形成される。
 送信部23は、例えば、シーケンス制御部25から入力された制御信号に基づいてRF信号をRFコイル部13に印加する。
 受信部24は、コイルシステム10の駆動により発生したMR信号を受信する。例えば、RFコイル部13で受けたMR信号を検波して、所要の信号処理を実行すると共に、A/D(Analog-to-Digital)変換することにより、デジタル化された複素データ(生データ)を生成する。むろん、受信部24は、検波したMR信号を直接にA/D変換して生データを生成してもよい。受信部24が生成した生データは、例えば、シーケンス制御部25に出力される。
 シーケンス制御部25は、MRI装置1を駆動するための、傾斜磁場電源22と送信部23と受信部24とを駆動する。例えば、シーケンス制御部25は、制御信号を傾斜磁場電源22、送信部23および受信部24に印加することにより、傾斜磁場電源22と送信部23と受信部24とを駆動する。
 この制御信号は、例えば、傾斜磁場電源22、送信部23および受信部24に印加するパルス電流の大きさ、印加時間および印加タイミングなどを規定したパルスシーケンスに従って生成される。このパルスシーケンスについての情報は、後述の情報処理装置26からシーケンス制御部25に入力される。また、シーケンス制御部25は、例えば、受信部24から入力された生データを情報処理装置26に出力する。
 制御システム20は、さらに、例えば、図1に示したように、情報処理装置26を備えている。この情報処理装置26は、例えば、演算部26A、入力部26B,表示部26C、記憶部26Dを含んで構成されている。
 入力部26Bは、例えば、ユーザからの情報をデジタルデータとして情報処理装置26内部に取り込む装置であり、例えばキーボードやマウス、スキャナなどによって構成されている。
 表示部26Cは、演算部26Aにより処理された結果(例えば形態画像)や、撮像条件などのデータ入力用のダイヤログ等を表示するためのもので、例えば液晶表示装置等のディスプレイデバイスで構成されている。
 記憶部26Dには、MRI装置1を制御する種々のプログラムが記憶されており、例えば、後述する判断処理に用いられるプログラム27や、後述する形態画像の作成に用いられる位相差強調画像化プログラム、等が記憶されている。
 演算部26Aは、例えば、各種プログラムの命令を解釈し、実行するためのものであり、例えばCPU(Central Processing Unit)により構成されている。この演算部26Aには、例えば、記憶部26Dに記憶されたプログラム27がMRI装置1の起動に併せてロードされ、それにより、演算部26Aは、例えば、ユーザからの指示に応じて、プログラム27の命令を解釈し、実行するようになっている。
 なお、演算部26Aは、各種プログラム(例えばプログラム27)の機能に対応したハードウェアによって構成されていてもよい。以下、演算部26Aにおいてプログラム27の各種命令が実行されることによって位相画像が作成されるものとする。
(2)位相差画像の作成: 
 次に、MRI装置1にて行われる位相差画像の作成について説明する。
 位相差画像を得るためには、強度画像と位相画像が必要である。これら画像を得る際は、GEのパルスシーケンスを用いることが好ましいが、例えば、GE系に含まれる他のパルスシーケンスや、GE系以外のパルスシーケンスを用いてもよい。
 GE系のパルスシーケンスによって得られた位相画像は、各ピクセル内に含まれる組織が作る局所磁場(外部磁場に比べた局所磁場)の変化量ΔBと、撮像に要したエコータイム(TE)の積(ΔB×TE)に比例している。そのため、位相画像から、大きな位相(差)情報を取り出すためには、TEを大きくしたり、ΔBを強調する関数(いわゆる強調関数)をより強いものに変えたりする。
 具体的には、所定のパルスシーケンスを用いて撮像を行う。撮像回数は、1回でもよいし、統計的に扱えるように複数回としてもよい。また、RFコイル部13が多チャンネルで構成され、チャンネルごとにMR信号が得られる場合には、チャンネルごとに得られたMR信号の相加平均を算出し、相加平均されたMR信号を用いて強度画像および位相画像を作成する。なお、相加平均を算出するにあたって、事前に、個々のMR信号に対して感度補正を行っておいてもよい。
(生データRの取得)
 図3は、位相差画像を作成するまでのデータ処理の流れを示した図である。同図に示すデータ処理は、演算部26Aの制御の下、制御システム20において実行されるものである。演算部26Aは、ユーザからの指示を受けて、同図に示すデータ処理の演算を開始する。
 演算部26Aは、まず、シーケンス制御部25に対して、所定のパルスシーケンスを用いて、生データを取得することを要求する制御信号を出力する。すると、シーケンス制御部25から、傾斜磁場電源22、送信部23および受信部24に対して、所定のパルスシーケンスに従った制御信号が出力される。
 この制御信号が入力されると、傾斜磁場電源22および送信部23は、コイルシステム10に対して所定の電流パルスを出力するため、RFコイル部13は、MR信号を検波する。ここで検波されるMR信号は、受信部24における所定の信号処理によって生データRに変換される。
 受信部24は、この生データRをシーケンス制御部25に入力し、シーケンス制御部25は、この生データRを演算部26Aに転送(入力)する。このようにして、演算部26AはMR信号に対応したデータ(生データR)を取得する。
(強度画像M(x),位相画像P(x)の作成)
 次に、演算部26Aは、シーケンス制御部25から入力された生データRを、不図示の内部メモリに設定したk空間に配置する。以下では、k空間に配置したデータをk空間データS(k)と呼ぶことにする。
 演算部26Aは、k空間データS(k)に対して逆フーリエ変換を施して複素画像ρ(x)を再構成する。複素画像ρ(x)は、下記式(1)で示されるように、リアル画像を実部に持ち、イマジナリ画像を虚部に持つ画像である。
Figure JPOXMLDOC01-appb-M000001
 演算部26Aは、複素画像ρ(x)から、強度画像M(x)と位相画像P(x)を得る。
(位相差画像PD(x)の作成)
 本実施の形態では、MR信号の取得に際して、長いTEを用いている。そのため、位相画像P(x)にフェーズラッピング(phase wrapping)が発生し、位相が2πを超えるものが、実際の位相から2πn(nは整数)を差し引いた位相値を取る。そのため、位相画像P(x)が縞模様の画像となり、本来の位相値を示さなくなる。そこで、演算部26Aは、このフェーズラッピングを取り除くと共に、位相差を取り出す処理を行う。
 具体的には、演算部26Aは、まず、複素画像ρ(x)に対してフーリエ変換を施して、複素画像ρ(x)を一度、k空間データS(k)に戻す。もしくは、k空間に配置しておいたk空間データS(k)を読み出す。次に、演算部26Aは、k空間データS(k)に対してLPF(Low Pass Filter;ローパスフィルタ)をかけ、それにより得られたデータL(k)×S(k)に対して逆フーリエ変換を施して、複素画像ρ’(x)を得る。なお、前記のL(k)はLPFの関数である。
 次に、演算部26Aは、複素画像ρ(x)およびρ’(x)を利用して、位相差画像PD(x)を作成する。具体的には、演算部26Aは、複素画像ρ(x)を複素画像ρ’(x)で除算し、複素商の演算を行うことにより、位相差画像PD(x)を作成する。これにより、位相部分のフェーズラッピングが取り除かれる。
 このとき、位相差画像PD(x)に含まれる位相差は、2πの幅を持っており、本実施の形態では、位相差画像PD(x)に含まれる位相差を、-π≦PD(x)<πと仮定する。位相差画像PD(x)に含まれる位相差の符号は、-γ×ΔB×TEで決定されるものである。
 ただし、LPFの定義と、位相差画像PD(x)を取り出すときの複素商の定義と、を変えることにより、位相差の大きさは変わらなくても、位相差画像PD(x)に含まれる位相差の符号が変わることがある。そこで、それに対応するために、特に静脈血に対しては符号が負となるように、前記の定義がなされている。なお、前記のγは、正の比例定数であり、例えば、水素の磁気回転比に相当するものである。
(3)判断処理: 
 次に、この位相差画像に対し、MRI装置1の利用者が関心を持つ領域(以下、関心領域と記載する。)を設定し、当該関心領域に含まれる位相差画像に対して所定の演算処理を行う。この演算処理により、関心領域に含まれる組織の正常性(又は、非正常性)を判断することができる。
 以下、図4を参照しつつ、関心領域に含まれる組織の正常性(又は、非正常性)を判断する処理について説明する。同図に示すデータ処理は、演算部26Aの制御の下、制御システム20において実行されるものである。演算部26Aは、位相差画像の作成が完了したとき、又は、ユーザからの指示を受けて、同図に示すデータ処理の演算を開始する。
[領域指定インターフェースの表示] 
 処理が開始されると、表示部26Cの画面に、MRI装置1にて撮像した断面もしくは空間の一部を指定するためのインターフェースを表示する(S1)。このインターフェースには、例えば、MR信号に基づいて作成された強度画像や位相画像等を表示し、画像上の所定領域を、閉曲線(枠線等)や座標値等の領域指定手段によって、指定できるようになっている。
[関心領域の設定] 
 次に、MRI装置1の操作者は、表示部26Cに表示されたインターフェースを介して、MRI装置1にて撮像した断面もしくは空間の一部を、領域指定する(S2)。ここで操作者によって指定される領域が関心領域である。操作者は、例えば、表示部26Cに表示された画像を観察したときに、非正常組織を含む可能性があると考えた領域を、関心領域に指定する。なお、関心領域は、二次元領域であっても良いし、三次元領域であってもよい。
[位相差分布を作成] 
 関心領域が設定されると、演算部26Aは、関心領域に含まれる組織から取得されたMR信号の全位相データを取得し、これら位相データを統計し、例えば、横軸を位相値とし、縦軸をデータ数とした位相差分布を作成する(S3)。
 次に、位相差分布のフィッティングに用いる関数群を選択する(S4)。ここでは、操作者が、目的とする組織、病変、病態、撮像法等に応じて、適切な関数群を選択する。なお、目的が予め決まっている場合は、ステップS4をスキップし、目的に応じた既定の関数群が自動的に適用されるようにしてもよい。また、目的とする組織、病変、病態、撮像法等に応じた既定の関数群を用意しておき、選択画面から目的を選択することにより、目的に応じた適切な関数群が選択されるようにしてもよい。
[フィッティング] 
 次に、1つの関心領域に対して作成された1つの位相差分布に対し、複数の関数によって同時にフィッティングを行う(S5)。即ち、複数の関数を重ね合わせた曲線にて、1つの位相差分布を近似する。複数の関数には様々な関数を採用可能であり、ガウス分布、ローレンツ分布、二項分布等が例示される。
 むろん、複数の関数は、分布関数に限る必要はなく、互いに直交する必要もない。また、複数の関数は、異なる種類の関数を組み合わせて構成されてもよい。なお、コンピュータの演算によってフィッティングを行う際は、複数の関数は、有限個である必要がある。
 ここで、フィッティングに用いる複数の関数の少なくとも1つは、ガウス分布を用いることが好ましい。ランダムな組織により得られる信号分布は、独立な多数の因子の和として表される確率変数であり、中心極限定理によってガウス分布に従うことが保証されるからである。
 フィッティングに用いる各関数は、関数形状を変更するための1乃至複数のパラメータを有するため、複数の関数を用いてフィッティングを行う際は、少なくとも、フィッティングに使用する関数の数以上のパラメータを調整する必要がある。
 実際には、複数のパラメータを適宜に変更しつつ、位相差分布と複数の関数の重ね合わせとが最も近似するパラメータセットを探索することになる。このとき、フィッティング関数がガウス分布であれば、パラメータが3つと少なく、パラメータの調整が容易である。なお、複数の関数と位相差分布との近似度合いは、例えば、最小二乗法等によって評価することができる。
 後述する実施例では、複数の関数として2つのガウス分布を用いる二重ガウス分布モデルを採用してある。ガウス分布は、高さ、中心位置(平均)、標準偏差σ(又は分散σ^2)、の3つのパラメータを有するため、二重ガウス分布モデルでは、6つのパラメータを調整しつつフィッティングを行うことになる。
 下記(2)式は、二重ガウス分布モデルにおいて、フィッティングに用いる関数である。
Figure JPOXMLDOC01-appb-M000002
 前記(2)式において、A1は第1のガウス分布の高さに対応し、B1は第1のガウス分布の分散の逆数に対応し、C1は第1のガウス分布の中心位置に対応する。また、A2は第2のガウス分布の高さに対応し、B2は第2のガウス分布の分散の逆数に対応し、C2は第2のガウス分布の中心位置に対応する。
[表示] 
 以上説明したフィッティングによって求められるパラメータセットは、関心領域に含まれる組織の磁化率を特徴付ける値の組み合わせとなる。即ち、関心領域が正常組織とは磁化率の異なる非正常組織を含有している場合のパラメータセットと、関心領域が正常組織のみを含有している場合のパラメータセットとは異なるものとなる。
 そこで、フィッティングによって求められたパラメータセットを、表示部26Cに表示したり(S6)、パラメータセットに基づいて、関心領域に含まれる組織の正常組織からの乖離度合い(組織の正常性(又は、非正常性))を算出し、この度合いを表示部26Cに表示したりする(S7)。このような表示を行えば、操作者は、関心領域に含まれる組織が正常か非正常かを判断する基準を得ることができるし、正常組織からの乖離度合いの目安を得ることができる。
[実験結果の説明] 
 ここで、マウスの脳実質から得られた位相差分布に対し、二重ガウス分布モデルでフィッティングを行った実験の結果について説明する。本実験では、ヒトでアルツハイマーを起こすアミノ酸の置換を持つ変異型APP(Amyloid precursor protein)の遺伝子改変マウスと、いろいろな性質が一定であるように遺伝的にコントロールされたコントロールマウスと、で、脳内位相差分布の検討を行った。
 具体的には、遺伝子改変マウスおよびコントロールマウスの脳を取り出した後、7T-MRI装置で3D-FLASHを用いて撮像を行い、強度画像と位相画像を同時に取得した。主な撮像条件は、TR/TE:50/12.8ms、Flip Angle:20°、Matrix Size:194×128×82(0.08mm isovoxel)、加算:24回、である。
 本実験では、脳の中でも、皮質、海馬および視床から得た強度画像と位相画像を用いている。皮質、海馬および視床は、アルツハイマーの病理変化の1つである老人斑が特に多く集積するためである。なお、老人斑内には、鉄が沈着することが知られており、この鉄による磁化率の変化を位相信号として捉えることができるものと考えられる。具体的には、老人斑があると考えられる部位を含むように関心領域を設定し、これらの脳実質から収集した各磁気共鳴信号の位相値の位相差分布に対し、二重ガウス分布モデルにてフィッティングを行った。
 図5は、コントロールマウスに設定した関心領域から取得した位相をプロットしたグラフであり、図6、7は、遺伝子改変マウスに設定した関心領域から取得した位相をプロットしたグラフである。これらのグラフでは、横軸が位相(ラジアン)を示し、縦軸が検出数(例えば、画素数)を示している。なお、位相値は、上述したように、フェーズラッピングが取り除かれ、-πからπの範囲に調整されているが、図5~7には、位相差分布の中心(位相0)付近を拡大して示してある。
 ここで、フィッティングカーブにおける各ガウス分布の寄与率(高さ)を比較すると、図5では、第1のガウス分布が第2のガウス分布の7~8倍の高さであるのに対し、図7では、第1のガウス分布と第2のガウス分布とがほぼ同じ高さである。しかし、図6でも、図5と同様に、第1のガウス分布が第2のガウス分布の8~9倍の高さを有している。即ち、各ガウス分布の寄与率(高さ)は、同じ病変であっても各症例で異なり、非正常組織の有無の判断基準に採用できないことが分かる。
 次に、ガウス分布の標準偏差(幅)を比較すると、図5では、第1のガウス分布と第2のガウス分布の双方がシャープな分布を示し、特に、第1のガウス分布に比べて第2のガウス分布は幅(例えば、半値幅)が狭くなっており、また、いずれのガウス分布も、位相差分布に比べて幅(例えば、半値幅)が狭い。一方、図6,7では、いずれも第1のガウス分布に比べて第2のガウス分布の幅(例えば、半値幅)が広く、また、第2のガウス分布は、位相差分布に比べて幅(例えば、半値幅)が広い。即ち、ガウス分布の標準偏差(幅)は、コントロールマウスの正常組織と遺伝子改変マウスの老人斑を含む組織とで、明確な差異が認められ、非正常組織の有無の判断基準になり得ることが分かる。
 以上の実験結果に基づくと以下の考察が為される。
 まず、図5~7に示すフィッティングカーブを構成する第1のガウス分布は、ノイズ成分に起因すると考えられる。MRI装置1では、電気的な熱雑音やカットオフにより生じるデジタライズノイズ等により位相差に誤差が発生することが当然に予測される。これらノイズはガウス分布に従うことが知られており、比較的高いガウス分布が発生すると考えられるためである。従って、複数の関数で位相差分布をフィッティングする場合、複数の関数の少なくとも1つに、ガウス分布を採用すると良好なフィッティング結果が得られると考えられる。
 また、図6,7に示すフィッティングカーブに現れる第2のガウス分布は、何らかの磁性物質を含む組織に関係すると考えられる。図6,7に示す位相差分布は、老人斑を含むように設定された関心領域から取得されている。老人斑は、上述したように磁性を有する鉄を含む組織であり、磁性物質を含む組織の位相差分布をガウス分布でフィッティングすると、比較的大きな標準偏差を有するためである。
 さらに、図5と図6,7とで位相差分布の裾を比べると、図5に示す位相差分布は、ガウス分布に比較的近いシャープな立ち上がりを有するが、図6,7に示す位相差分布は、同じ高さのガウス分布と比べると裾が広く、何らかのブロードな位相差分布が混在していることを示唆している。即ち、位相差分布の裾が、同程度の高さのガウス分布に比べて広い場合は、関心領域に何らかの磁性物質を含む組織が含まれていることを示唆する。
 また、中心が+側にずれたガウス分布がフィッティングされる場合は、関心領域に、位相が+側にシフトした磁性物質を含む組織があると考えられ、中心が-側にずれたガウス分布がフィッティングされる場合は、関心領域に、位相が-側にシフトした磁性物質を含む組織があると考えられる。このことから、正常組織の磁化率の分布と、非正常組織の磁化率の分布の重ね合わせが、位相画像の分布になっていると考えることができる。
 さらに、位相差分布と第2のガウス分布との中心のずれは、図5に比べると図6,7の方が大きい。その理由は、正常組織が位相差0であると考えると、この位相差0に近い中心を有するガウス分布は正常組織の示すガウス分布であると考えられ、この位相差0から離れた中心を有するガウス分布ほど、正常組織でない組織を多く含む可能性の高い組織と考えられる。即ち、正常組織に対応するガウス分布の中心と、正常でない組織に対応するガウス分布の中心との距離は、関心領域に含まれる組織の正常性(又は、非正常性)を示す目安になると考えられる。
(4)形態画像の作成: 
 次に、二重ガウス分布モデルを用いて得られるフィッティング関数を利用して行う、特定の位相差を強調する位相差強調画像化法について説明する。なお、位相差強調画像化法とは、得られた位相差画像PD(x)の一部を選択し、その一部もしくは全部を取捨選択し、強調関数w(θ)によって強調することによって選択された位相情報に対応する画像の部分を強度画像M(x)上に表現する方法である。位相差強調画像化法によれば、強調関数w(θ)によって作成された画像を提供することができる。
 図8は、二重ガウス分布モデルにより得られるフィッティング関数と、強調すべき位相差との関係を示してある。なお、図8には、図6の実験結果を例に取り、強調すべき位相範囲を示している。図8では、第2のガウス分布と第1のガウス分布が交差する位相をθ1とし、第2のガウス分布が第1のガウス分布に比べて大きくなる位相範囲(θ>θ1)を、強調範囲としてある。
 二重ガウス分布モデルにおいては、上述したように、関心領域に磁性物質を含む非正常組織が含まれる場合、第2のガウス分布が第1のガウス分布に比べて標準偏差が大きく、広い裾野を有する。従って、非正常組織の分布を示唆する第2のガウス分布が正常組織の分布を示唆する第1のガウス分布を超える位相を有する組織を強調することにより、非正常組織の視認容易性を向上した形態画像を提供できると考えられる。
(強調画像w(PD(x))の作成)
 図10は、形態画像を作成するまでのデータ処理の流れを示した図である。まず、演算部26Aは、位相差画像PD(x)の位相θがθ1より大きくなる位相を選択し、その選択した位相θを強調する。
 位相差分布は、例えば、LPFサイズの変更により、例えば、図9(A)~(C)に示したモデルのように変化する。即ち、LPFサイズの増加とともに、位相差ゼロを中心としたほぼ対称な位相差分布の幅が減少する。このため、位相差画像は位相差ゼロおよびその近傍の値を主に持ち、位相差画像として組織のコントラストが付けづらくなる。一方、細かい構造を持つ組織においては、サイズの大きなLPFを用いた方が、サイズの小さなLPFを用いたときよりもコントラストが付け易くなるという側面もある。
 このような現象の単純な理解は、位相差画像を得るに際して用いられる複素商が指数上の減算に相当することを考慮すれば容易である。即ち、exp(φ-φ´)(φ:フェーズラップを含む全位相、φ´:フェーズラップを主に含むローパス成分(低周波成分)の位相)において、LPFサイズが全画像サイズとなっている場合には、全ての周波数を含むことを意味するので、φ=φ´となり、位相差はゼロに近づくはずである。このため、位相差分布の幅が狭くなっていく。
 しかし、これはあくまでも極限傾向であり、一般的な組織において、どのフィルタを用いた場合にどの分布になるかを予測するのは事実上不可能である。このような位相差分布の変化に対応するために、本実施の形態では、後述する強調関数として柔軟に対応可能なものを選び出し、図9(A)~(C)に示したような、同じ位相差に対して異なる度数分布を示すものに対して、同じコントラストを付けられるように工夫を行っている。
 例えば、演算部26Aは、目的とする組織のコントラストが所望の大きさとなるように、位相差の幅および、その中心値を選択する。なお、位相差の幅および、その中心値の選択については、演算部26Aに任せず、MRI装置1のユーザに委ねてもかまわない。その場合には、演算部26Aは、目的とする組織の位相θを、フィルタ処理による位相差分布の変化が加味された上で選択することになる。
 続いて、演算部26Aは、その選択した位相θを強調関数w(θ)で強調することにより、強調画像w(PD(x))を得る。
 ここで、強調関数w(θ)としては指数関数が用いられる。本実施の形態では、指数関数の一例として、π関数を用いる。このπ関数は、以下の2つの式で表される。
w(θ)=1…(-σ≦θ≦σ)
w(θ)=exp(-a×(Abs(θ)-σ))…(θが前記の範囲以外のとき)
(a、b、σの決定)
 パラメータa、b、σはいずれも、実数の値をとる。パラメータa、bは、位相差強調の度合いを調整するものであり、LPFのフィルタサイズによって決定される。パラメータa、bは、また、目的とする組織と、そのバックグラウンドとのコントラストCもしくはコントラスト・ノイズ比CNRを最大にするように決定される。パラメータσは、位相差画像PD(x)上のノイズを低減するものであり、例えば、位相差画像PD(x)上の位相平均値が0(ゼロ)付近をとる組織の標準偏差によって決定される。パラメータσは、多くの実験から得られるデータに基づいて求めることが可能である。ただし、位相平均値が0(ゼロ)付近をとる組織が一度に撮像された全ての位相差画像PD(x)上に存在しない場合もある。その場合には、パラメータσは、例えば、コントラストCまたはコントラスト・ノイズ比CNRによって決定される。
 なお、上述したコントラストCは、以下の式で示したように、画像上の位置x1にある目的とする組織の強調後の信号w(PD(x1))×M(x1)と、画像上の位置x2にある強調された組織のバックグラウンドの信号w(PD(x2))×M(x2)との差分を絶対値で表したものによって決定される。
 また、コントラスト・ノイズ比CNRについては、以下の式で示したように、C/σ’で表される。なお、σ’は、目的とする組織の強調画像上の標準偏差か、または目的とする組織のバックグラウンドの強調画像上の標準偏差によって決定される。
 ただし、σ’の決定に際して、上述したいずれの標準偏差も採用できないことがある。その場合には、被検体2の外部信号の標準偏差か、または差分法による標準偏差を採用するものとする。このようにして決定されたσによって規定されたπ関数を用いることにより、位相差が持つノイズ部分のみを除去することができ、高いS/N比の位相差強調画像を作成することが可能となる。
C=Abs(w(PD(x1))×M(x1)-w(PD(x2))×M(x2))
CNR=C/σ’
 また、Abs(θ)は、θの絶対値を表している。π関数は、Abs(θ)≦σの範囲でPD(x)を強調せず、その他の範囲で位相差画像PD(x)の強調を行う関数である。このπ関数は、任意の冪関数を任意の精度で近似することの可能なものであることから、強調関数として多項式を用いた場合と比べて、より柔軟な強調を行うことが可能である。
 例えば、撮像領域のサイズに合わせてLPFのフィルタサイズを変更すると、これに伴って位相差の分布も若干、変化する。一方、撮像領域のサイズを一定に固定した状態で、LPFのフィルタサイズを変更すると、これに伴って位相差の分布が大幅に変化する。このように、撮像条件などが互いに異なる場合には、位相差の分布も互いに異なる。従って、撮像条件などとは無関係に常に同一の強調関数を用いた場合には、それぞれのコントラストが変わってしまい、目的とした組織の強調を確実に行うことができなくなってしまう。一方、本実施の形態では、上述したように、強調関数のパラメータa、b、σは実数の値をとり、これらを撮像条件などに応じて柔軟に変更することが可能である。これにより、撮像条件を変えた場合であっても、コントラストを同一もしくは同等に維持することが可能である。
(形態画像I(x)の作成)
 次に、演算部26Aは、例えば、所定のモード(ルール)に従って、強調画像w(PD(x))で強度画像M(x)をマスクし、それにより形態画像I(x)を得る。強調画像w(PD(x))で強度画像M(x)をマスクする際の具体的な条件は、強調したい対象に応じて設定することが可能なものであり、基本的には以下に例示した4種類(組織強調、血管強調、全強調、構造強調)の強調モードに対応して設定される。
 例えば、パラメータσと位相差画像PD(x)とが共に正であるか負であるかで、組織強調および血管強調のいずれか一方を選択することが可能である。全強調については、パラメータσおよび位相差画像PD(x)の符号に依らないが、構造強調の場合には、例えば、皮質内構造が作り出す位相差α(以下、単に位相差αと称する)の値をあらかじめ実験で求めておき、位相差αと位相差画像PD(x)との大小関係に応じて条件式を設定する。
強調モードA(組織強調)
I(x)=w(PD(x))×M(x)…(PD(x)≧0)
I(x)=M(x)…(PD(x)<0)
強調モードB(血管強調)
I(x)=w(PD(x))×M(x)…(PD(x)≦0)
I(x)=M(x)…(PD(x)>0)
強調モードC(全強調)
I(x)=w(PD(x))×M(x)
強調モードD(構造強調)
|α|≦σ
PD(x)≦0のとき
I(x)=w(PD(x))×M(x)…(-|α|≦PD(x)≦-σ)
I(x)=M(x)…(PD(x)<-|α|)
PD(x)>0のとき
I(x)=w(PD(x))×M(x)
 強調モードAでは、演算部26Aは、位相差画像PD(x)が0(ゼロ)以上となる部分を取り出し、その部分についてだけ強調を行うことにより、形態画像I(x)を作成する。このとき、位相差画像PD(x)が負となる部分については、強調を行わない。強調モードAでは、形態画像I(x)は、組織強調された画像となる。
 強調モードBでは、演算部26Aは、位相差画像PD(x)が0(ゼロ)以下となる部分を取り出し、その部分についてだけ強調を行うことにより、形態画像I(x)を作成する。このとき、位相差画像PD(x)が正となる部分については、強調を行わない。強調モードBでは、形態画像I(x)は、血管強調された画像となる。
 強調モードCでは、演算部26Aは、位相差画像PD(x)全体を強調するにより、形態画像I(x)を作成する。強調モードCでは、形態画像I(x)は、組織や血管などを含む全体が強調された画像となる。
 強調モードDでは、演算部26Aは、位相差画像PD(x)が0(ゼロ)以下となっているときに、-|α|≦PD(x)≦-σを満たす部分を強調することにより、形態画像I(x)を作成する。このとき、PD(x)<-|α|を満たす部分については、強調を行わない。また、演算部26Aは、位相差画像PD(x)が0(ゼロ)よりも大きいときにも、前記の条件式を用いて、形態画像I(x)を作成する。強調モードDでは、形態画像I(x)は、構造強調された画像となる。なお、位相差αは、皮質内構造によって作り出されたものであることから、形態画像I(x)は、実際には、皮質構造が強調された画像となる。
 このように、本実施の形態では、強調画像w(PDr(x))、w(PDa(x))の、強度画像M(x)へのマスクの仕方を変えることにより、各形態画像I(x)が示す組織コントラストを背景とした、より正確な脳機能の解剖学的位置を特定することができる。
(5)まとめ: 
 以上説明した実施形態によれば、生体から得られたMR画像を解析するにあたり、関心領域から得たMR信号の位相差分布を作成し、この位相差分布を複数の関数群で同時にフィッティングし、位相差分布に対してフィッティングされた前記複数の関数群のパラメータによって特定される、関心領域に含まれる組織の磁化率に基づいて、関心領域に含まれる生体の正常性を検証する。よって、MRIを含む医療用画像機器で描出できなかった組織や病変を、現在の保険診療の範囲内で現実的に検知することが可能となる。
 なお、本技術は上述した実施形態や変形例に限られず、上述した実施形態や変形例の中で開示した各構成を相互に置換したり組み合わせを変更したりした構成、公知技術並びに上述した実施形態や変形例の中で開示した各構成を相互に置換したり組み合わせを変更したりした構成、等も含まれる。また,本発明の技術的範囲は上述した実施形態に限定されず、特許請求の範囲に記載された事項とその均等物まで及ぶものである。
 本技術は、MRIを用いた各種診断に用いることができると考えられる。特に、アミロイドβのような微小な病変組織を検出する際に、大きく貢献できる。組織中に散漫に存在するものを検出するときにも、本手法は同様に適用可能であり、例えば、健康診断等で必ず検査対象となる肝臓中に含まれる脂肪の割合などを定量的に計測できると考えられる。
1…MRI装置、2…被検体、10…コイルシステム、10A…ボア、11…静磁場コイル部、12…傾斜磁場コイル部、13…RF(Radio Frequency)コイル部、20…制御システム、21…静磁場電源、22…傾斜磁場電源、23…送信部、24…受信部、25…シーケンス制御部、26…情報処理装置、26A…演算部、26B…入力部、26C…表示部、26D…記憶部、27…プログラム、30…寝台部、13-1~13-8…コイル

Claims (5)

  1.  生体から得られた核磁気共鳴画像を解析するための画像解析装置であって、
     生体の所定領域から得た核磁気共鳴信号の位相差分布を作成する位相差分布作成部と、
     位相差分布作成部により作成された前記位相差分布に対し、複数の関数群でフィッティングするフィッティング部と、
     前記フィッティング部により前記位相差分布に対してフィッティングされた前記複数の関数群のパラメータによって特定される前記所定領域に含まれる組織の磁化率に基づいて、前記所定領域に含まれる生体の正常性を検証する検証部と、
    を備えることを特徴とする、画像解析装置。
  2.  前記複数の関数群は、少なくとも1つのガウス分布を含む請求項1に記載の画像解析装置。
  3.  前記複数の関数群は、2つのガウス分布で構成され、
     前記2つのガウス分布のうち、標準偏差の大きいガウス分布が他方のガウス分布に比べて大きい位相を有する画素を強調した核磁気共鳴画像を作成する画像生成部を、更に備える請求項1又は請求項2に記載の画像解析装置。
  4.  生体から得られた核磁気共鳴画像を解析するための画像解析方法であって、
     生体の所定領域から得た核磁気共鳴信号の位相差分布を作成する位相差分布作成工程と、
     前記位相差分布に対し、複数の関数群でフィッティングするフィッティング工程と、
     前記位相差分布に対してフィッティングされた前記複数の関数群のパラメータによって特定される前記所定領域に含まれる組織の磁化率に基づいて、前記所定領域に含まれる生体の正常性を検証する検証工程と、
    を備えることを特徴とする、画像解析方法。
  5.  生体から得られた核磁気共鳴画像を解析する機能をコンピュータに実現させるための画像解析プログラムであって、
     生体の所定領域から得た核磁気共鳴信号の位相差分布を作成する位相差分布作成機能と、
     前記位相差分布に対し、複数の関数群でフィッティングするフィッティング機能と、
     前記位相差分布に対してフィッティングされた前記複数の関数群のパラメータによって特定される前記所定領域に含まれる組織の磁化率に基づいて、前記所定領域に含まれる生体の正常性を検証する検証機能と、
    を備えることを特徴とする、画像解析プログラム。
PCT/JP2012/074701 2011-09-28 2012-09-26 画像解析装置、画像解析方法及び画像解析プログラム WO2013047583A1 (ja)

Priority Applications (6)

Application Number Priority Date Filing Date Title
KR1020147011275A KR101685377B1 (ko) 2011-09-28 2012-09-26 화상해석장치, 화상해석방법 및 화상해석프로그램을 기록한 기록매체
JP2013536329A JP6041356B2 (ja) 2011-09-28 2012-09-26 画像解析装置、画像解析装置の作動方法及び画像解析プログラム
CN201280046660.2A CN103826535B (zh) 2011-09-28 2012-09-26 图像分析装置、图像分析方法
AU2012317518A AU2012317518B2 (en) 2011-09-28 2012-09-26 Image analysis device, image analysis method, and image analysis programme
US14/347,703 US9330457B2 (en) 2011-09-28 2012-09-26 Device, method, and program for analyzing a magnetic resonance image using phase difference distribution
EP12836937.8A EP2762071B1 (en) 2011-09-28 2012-09-26 Image analysis device, image analysis method, and image analysis programme

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
JP2011213535 2011-09-28
JP2011-213535 2011-09-28

Publications (1)

Publication Number Publication Date
WO2013047583A1 true WO2013047583A1 (ja) 2013-04-04

Family

ID=47995613

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2012/074701 WO2013047583A1 (ja) 2011-09-28 2012-09-26 画像解析装置、画像解析方法及び画像解析プログラム

Country Status (7)

Country Link
US (1) US9330457B2 (ja)
EP (1) EP2762071B1 (ja)
JP (1) JP6041356B2 (ja)
KR (1) KR101685377B1 (ja)
CN (1) CN103826535B (ja)
AU (1) AU2012317518B2 (ja)
WO (1) WO2013047583A1 (ja)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2021177353A1 (ja) * 2020-03-04 2021-09-10 国立大学法人熊本大学 画像処理方法、画像処理装置、プログラム及び記録媒体
WO2022034691A1 (ja) * 2020-08-14 2022-02-17 株式会社Splink コンピュータプログラム、情報処理装置、情報処理方法、学習済みモデル生成方法及び相関画像出力装置

Families Citing this family (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20150134261A1 (en) * 2013-11-14 2015-05-14 J. Michael O'Connor Synchronization of patient motion detection equipment with medical imaging systems
JP6595393B2 (ja) * 2016-04-04 2019-10-23 株式会社日立製作所 磁気共鳴イメージング装置、及び、画像処理方法
EP3550318A1 (de) * 2018-04-06 2019-10-09 Siemens Healthcare GmbH Verfahren zur erzeugung einer b0-karte mittels eines magnetresonanz-fingerabdruckverfahrens
JP7177621B2 (ja) * 2018-08-02 2022-11-24 キヤノンメディカルシステムズ株式会社 磁気共鳴イメージング装置
JP7455508B2 (ja) * 2018-12-26 2024-03-26 キヤノンメディカルシステムズ株式会社 磁気共鳴イメージング装置および医用複素数画像処理装置
ES2773333B2 (es) * 2019-01-10 2021-07-27 Univ Valencia Politecnica Metodo y sistema de generacion de senales de resonancia magnetica por rotacion rapida con angulo magico de campos con codificacion espacial
CN114224298B (zh) * 2022-01-17 2023-12-01 中国科学院电工研究所 一种核磁共振下的磁声电成像系统及方法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2005028151A (ja) 2003-07-10 2005-02-03 General Electric Co <Ge> 磁気共鳴イメージングを用いて脳内鉄を検出するためのシステム及び方法
JP2009125582A (ja) * 2007-11-21 2009-06-11 Toshiba Corp 磁気共鳴イメージング装置
WO2010073923A1 (ja) * 2008-12-26 2010-07-01 国立大学法人 熊本大学 位相差強調画像化法(Phase Difference Enhanced Imaging;PADRE)、機能画像作成法、位相差強調画像化プログラム、位相差強調画像化装置、機能画像作成装置および磁気共鳴画像化(Magnetic Resonance Imaging;MRI)装置

Family Cites Families (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5909119A (en) * 1995-08-18 1999-06-01 Toshiba America Mri, Inc. Method and apparatus for providing separate fat and water MRI images in a single acquisition scan
JP4262518B2 (ja) 2003-05-21 2009-05-13 ジーイー・メディカル・システムズ・グローバル・テクノロジー・カンパニー・エルエルシー 磁気共鳴撮影装置
US10098563B2 (en) 2006-11-22 2018-10-16 Toshiba Medical Systems Corporation Magnetic resonance imaging apparatus
JP4936864B2 (ja) * 2006-11-22 2012-05-23 株式会社東芝 磁気共鳴イメージング装置
JP4945225B2 (ja) 2006-11-30 2012-06-06 株式会社東芝 画像処理装置及びプログラム
JP5624273B2 (ja) 2007-11-22 2014-11-12 株式会社東芝 磁気共鳴イメージング装置
JP5336731B2 (ja) 2007-12-17 2013-11-06 株式会社日立メディコ 磁気共鳴イメージング装置

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2005028151A (ja) 2003-07-10 2005-02-03 General Electric Co <Ge> 磁気共鳴イメージングを用いて脳内鉄を検出するためのシステム及び方法
JP2009125582A (ja) * 2007-11-21 2009-06-11 Toshiba Corp 磁気共鳴イメージング装置
WO2010073923A1 (ja) * 2008-12-26 2010-07-01 国立大学法人 熊本大学 位相差強調画像化法(Phase Difference Enhanced Imaging;PADRE)、機能画像作成法、位相差強調画像化プログラム、位相差強調画像化装置、機能画像作成装置および磁気共鳴画像化(Magnetic Resonance Imaging;MRI)装置

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
SHINGO KAKEDA ET AL.: "A novel tract imaging technique of the brainstem using phase difference enhanced imaging: normal anatomy and initial experience in multiple system atrophy", EUROPEAN RADIOLOGY, vol. 21, no. 10, 1 October 2011 (2011-10-01), pages 2202 - 2210, XP019946953 *
T.YONEDA ET AL.: "Triple-layer Appearance of Human Cerebral Cortices on Phase-Difference Enhanced Imaging using 3D Principle of Echo Shifting with a Train of Observations (PRESTO) Sequence", PROCEEDINGS OF INTERNATIONAL SOCIETY FOR MAGNETIC RESONANCE IN MEDICINE, 18 April 2009 (2009-04-18), XP055064621 *

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2021177353A1 (ja) * 2020-03-04 2021-09-10 国立大学法人熊本大学 画像処理方法、画像処理装置、プログラム及び記録媒体
JP7504340B2 (ja) 2020-03-04 2024-06-24 国立大学法人 熊本大学 画像処理方法、画像処理装置、プログラム及び記録媒体
WO2022034691A1 (ja) * 2020-08-14 2022-02-17 株式会社Splink コンピュータプログラム、情報処理装置、情報処理方法、学習済みモデル生成方法及び相関画像出力装置
US11972564B2 (en) 2020-08-14 2024-04-30 Splink, Inc. Recording medium, information processing device, information processing method, trained model generation method, and correlation image output device

Also Published As

Publication number Publication date
EP2762071A1 (en) 2014-08-06
KR20140084081A (ko) 2014-07-04
CN103826535A (zh) 2014-05-28
AU2012317518B2 (en) 2015-10-08
CN103826535B (zh) 2016-07-06
EP2762071A4 (en) 2015-04-29
EP2762071B1 (en) 2021-03-24
KR101685377B1 (ko) 2016-12-12
JPWO2013047583A1 (ja) 2015-03-26
AU2012317518A1 (en) 2014-05-15
US20140233825A1 (en) 2014-08-21
US9330457B2 (en) 2016-05-03
JP6041356B2 (ja) 2016-12-07

Similar Documents

Publication Publication Date Title
JP6041356B2 (ja) 画像解析装置、画像解析装置の作動方法及び画像解析プログラム
JP4982881B2 (ja) 位相差強調画像化法(PhaseDifferenceEnhancedImaging;PADRE)、機能画像作成法、位相差強調画像化プログラム、位相差強調画像化装置、機能画像作成装置および磁気共鳴画像化(MagneticResonanceImaging;MRI)装置
JP6434030B2 (ja) Dixonタイプ水/脂肪分離する磁気共鳴イメージング
US9928596B2 (en) Motion corrected imaging system
JP5982494B2 (ja) 磁気共鳴イメージング装置
WO2013054718A1 (ja) 磁気共鳴イメージング装置および磁化率強調画像生成方法
WO2012174177A2 (en) Systems and methods for imaging and quantifying tissue magnetism with magnetic resonance imaging
Koolstra et al. Cartesian MR fingerprinting in the eye at 7T using compressed sensing and matrix completion‐based reconstructions
US11796617B2 (en) System and method for reconstruction of magnetic resonance images acquired with partial Fourier acquisition
WO2013098060A1 (en) Mri with dixon-type water/fat separation and prior knowledge about inhomogeneity of the main magnetic field
WO2011139745A2 (en) Systems and methods for susceptibility tensor imaging
Hutton et al. Phase informed model for motion and susceptibility
WO2014038441A1 (ja) 磁気共鳴イメージング装置および磁気共鳴イメージング方法
JP2019122623A (ja) 磁気共鳴イメージング装置および医用画像処理装置
Tan et al. Denoising and Multiple Tissue Compartment Visualization of Multi‐b‐Valued Breast Diffusion MRI
US20160146918A1 (en) Corrected magnetic resonance imaging using coil sensitivities
JP5650724B2 (ja) 磁気共鳴イメージング装置
US8493067B2 (en) Magnetic resonance system and method to create a magnetic resonance image data set by radial scanning of a magnetic resonance system
WO2013147281A1 (ja) 発熱分布情報を生成する装置と方法、磁気共鳴画像装置及びプログラム
US10955510B2 (en) Magnetic resonance imaging apparatus
JP6855239B2 (ja) 磁気共鳴イメージング装置
Olsson et al. Bias field correction of MPRAGE by an external reference--The poor man's MP2RAGE
US20240331225A1 (en) Methods, apparatuses, systems, and computer-readable mediums for biomarker quantification using free-breathing stack-of-radial imaging
JP2013034661A (ja) 磁気共鳴イメージング装置および磁気共鳴スペクトロスコピー撮像方法
Lira Chemical Species Separation in MRI with Simultaneous Estimation of Field Map and T2* Using a k-Space Formulation

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: 12836937

Country of ref document: EP

Kind code of ref document: A1

ENP Entry into the national phase

Ref document number: 2013536329

Country of ref document: JP

Kind code of ref document: A

WWE Wipo information: entry into national phase

Ref document number: 14347703

Country of ref document: US

NENP Non-entry into the national phase

Ref country code: DE

WWE Wipo information: entry into national phase

Ref document number: 2012836937

Country of ref document: EP

ENP Entry into the national phase

Ref document number: 20147011275

Country of ref document: KR

Kind code of ref document: A

ENP Entry into the national phase

Ref document number: 2012317518

Country of ref document: AU

Date of ref document: 20120926

Kind code of ref document: A