Magnetic resonance imaging device and image processing device

Abstract

A magnetic resonance imaging device (100) according to the present invention comprises a sequence control unit (120), a first calculation unit (150d), a range selection unit (150h), a second calculation unit (150e), and a display control unit (150c). The sequence control unit (120) executes a pulse sequence and collects data from an imaging target. On the basis of the data, the first calculation unit (150d) calculates a first map indicating an amount corresponding to volume magnetic susceptibility. On the basis of an echo time (TE) value and a phase change value of the first map, the range selection unit (150h) determines a selection range for a phase slope value. On the basis of the first map and the selection range, the second calculation unit (150e) calculates a second map indicating a feature amount of a voxel value of the first map or a region of a part of the first map. The display control unit (150c) displays an image on a display unit on the basis of the second map.

Description

Magnetic resonance imaging device and image processing device

[0001] The embodiments disclosed in this specification and the drawings relate to a magnetic resonance imaging apparatus and an image processing apparatus.

[0002] T2 using GRE (Gradient Echo) in magnetic resonance imaging (hereinafter sometimes referred to as "MRI"). ＊ In imaging, the phase component of the obtained MRI signal is proportional to the change in the local magnetic field, and this change in the local magnetic field is due to magnetic substances contained in the tissue, such as water and iron. Imaging methods that utilize this phase component include the QSM (Quantitative Susceptibility Map) method, which maps magnetic susceptibility, and the SWI (Susceptibility Weighted Imaging) method, and the PADRE (Phase-Difference-Enhanced) method, which selectively utilizes tissue-specific phase components to create images.

[0003] Although these technologies all utilize magnetic susceptibility obtained from phase information, from the perspective of clinical application, they cannot be said to have sufficient diagnostic ability for early dementia diagnosis, for example, in the asymptomatic stage of preclinical Alzheimer's disease.

[0004] U.S. Patent No. 6,501,272 U.S. Patent No. 6,658,280

[0005] Y. Wang, et al, “Quantitative Susceptibility Mapping (QSM): Decoding MRI Data for a Tissue Magnetic Biomarker”, Magnetic Resonance in Medicine 73:82―101 (2015)

[0006] One of the problems that the embodiments disclosed in this specification and the drawings aim to solve is to provide images with sufficient diagnostic capability for Alzheimer's disease and the like.

[0007] A magnetic resonance imaging apparatus according to an embodiment includes a sequence controller, a first calculator, a range selector, a second calculator, and a display controller. The sequence controller executes a pulse sequence to collect data from an imaging target. The first calculator calculates a first map indicating the magnitude of a quantity corresponding to volume magnetic susceptibility based on the data. The range selector determines a selection range of phase slope values ​​based on a phase change value and an echo time (TE) value in the first map. The second calculator calculates a second map indicating feature quantities of voxel values ​​of the first map or a partial region within the first map based on the first map and the selection range. The display controller causes a display unit to display an image based on the second map.

[0008] An image processing apparatus according to an embodiment includes a first calculator, a range selector, a second calculator, and a display controller. The first calculator calculates a first map indicating the magnitude of a quantity corresponding to volume magnetic susceptibility based on data collected from an imaging target by executing a pulse sequence. The range selector determines a selection range of phase slope values ​​based on a phase change value and an echo time (TE) value in the first map. The second calculator calculates a second map indicating feature quantities of voxel values ​​of the first map or a partial region within the first map based on the first map and the selection range. The display controller analyzes the second map and displays an image on a display.

[0009] FIG. 1 is a diagram illustrating an example of the configuration of a magnetic resonance imaging apparatus 100 according to an embodiment. FIG. 2 is a flowchart illustrating the flow of processing performed by the magnetic resonance imaging apparatus 100 according to an embodiment. FIG. 3 is a flowchart illustrating the processing of steps S100 and S200 in FIG. 2 in more detail. FIG. 4 is a diagram illustrating an example of a phase image acquired by the magnetic resonance imaging apparatus 100 according to an embodiment. FIG. 5 is a diagram illustrating the processing of creating a slope map by the magnetic resonance imaging apparatus 100 according to an embodiment. FIG. 6 is a diagram illustrating an example of an image generated by the magnetic resonance imaging apparatus 100 according to an embodiment. FIG. 7 is a diagram illustrating the linear transformation performed by the magnetic resonance imaging apparatus 100 according to an embodiment. FIG. 8 is a diagram illustrating the linear transformation performed by the magnetic resonance imaging apparatus 100 according to an embodiment. FIG. 9 is a flowchart illustrating the processing of steps S300 and S400 in FIG. 2 in more detail. FIG. 10 is a diagram illustrating an example of a histogram generated by the magnetic resonance imaging apparatus 100 according to an embodiment. FIG. 11 is a diagram showing an example of a histogram of slope map values, with the vertical axis representing the number of voxels and the horizontal axis representing the phase slope value. FIG. 12 is a diagram showing an example of a color bar controlled by a display control unit 150c according to an embodiment. FIG. 13 is a diagram showing an example of an image displayed by the magnetic resonance imaging apparatus 100 according to an embodiment. FIG. 14 is a diagram showing an example of an image displayed by the magnetic resonance imaging apparatus 100 according to an embodiment. FIG. 15 is a diagram showing an example of a histogram created by the magnetic resonance imaging apparatus 100 according to an embodiment. FIG. 16 is a diagram showing an example of processing performed by the magnetic resonance imaging apparatus 100 according to an embodiment. FIG. 17 is a diagram showing an example of processing performed by the magnetic resonance imaging apparatus 100 according to an embodiment. FIG. 18 is a diagram showing an example of processing performed by the magnetic resonance imaging apparatus 100 according to an embodiment. FIG. 19 is a diagram showing an example of processing performed by the magnetic resonance imaging apparatus 100 according to an embodiment. FIG. 20 is a diagram showing an example of processing performed by the magnetic resonance imaging apparatus 100 according to an embodiment.Fig. 21 is a diagram showing an example of processing performed by the magnetic resonance imaging apparatus 100 according to the embodiment. Fig. 22 is a diagram showing an example of processing performed by the magnetic resonance imaging apparatus 100 according to the embodiment. Fig. 23 is a diagram showing an example of an image displayed by the magnetic resonance imaging apparatus 100 according to the embodiment. Fig. 24 is a diagram showing an example of an image displayed by the magnetic resonance imaging apparatus 100 according to the embodiment.

[0010] (Embodiments) Hereinafter, embodiments of a magnetic resonance imaging apparatus and an image processing apparatus will be described in detail with reference to the drawings.

[0011] FIG. 1 shows an example of the configuration of a magnetic resonance imaging apparatus 100 according to an embodiment. The magnetic resonance imaging apparatus 100 includes a static magnetic field magnet 101, a static magnetic field power supply (not shown), a gradient magnetic field coil 103, a gradient magnetic field power supply 104, a bed 105, a bed control circuit 106, a transmission coil 107, a transmission circuit 108, a reception coil 109, a reception circuit 110, a sequence control unit 120, and an image processing device 130. Note that the magnetic resonance imaging apparatus 100 does not include a subject P (e.g., a human body). The configuration shown in FIG. 1 is merely an example. For example, the components of the sequence control unit 120 and the image processing device 130 may be integrated or separated as appropriate.

[0012] The static magnetic field magnet 101 is a magnet formed in a hollow, approximately cylindrical shape, and generates a static magnetic field in the space inside the cylinder in the direction of its central axis (Z-axis). The static magnetic field magnet 101 is, for example, a superconducting magnet, and is excited by receiving a current from a static magnetic field power supply. The static magnetic field power supply supplies a current to the static magnetic field magnet 101. As another example, the static magnetic field magnet 101 may be a permanent magnet, in which case the magnetic resonance imaging apparatus 100 does not need to include a static magnetic field power supply. Furthermore, the static magnetic field power supply may be provided separately from the magnetic resonance imaging apparatus 100.

[0013] The gradient magnetic field coil 103 is a hollow, approximately cylindrical coil and is disposed inside the static magnetic field magnet 101. The gradient magnetic field coil 103 is formed by combining three coils corresponding to the mutually orthogonal X, Y, and Z axes. These three coils are individually supplied with current from a gradient magnetic field power supply 104 to generate gradient magnetic fields along the X, Y, and Z axes, whose magnetic field strength in the Z direction varies depending on the distance from the center of each axis. The gradient magnetic fields of the X, Y, and Z axes generated by the gradient magnetic field coil 103 are, for example, a slice gradient magnetic field Gs, a phase encoding gradient magnetic field Ge, and a readout gradient magnetic field Gr. The gradient magnetic field power supply 104 supplies current to the gradient magnetic field coil 103.

[0014] The bed 105 includes a top plate 105a on which the subject P is placed, and under the control of a bed control circuit 106, the top plate 105a is inserted into the cavity (imaging port) of the gradient magnetic field coil 103 with the subject P placed thereon. The bed 105 is usually installed so that its longitudinal direction is parallel to the central axis of the static magnetic field magnet 101. Under the control of the image processing device 130, the bed control circuit 106 drives the bed 105 to move the top plate 105a in the longitudinal direction and up and down.

[0015] The transmitting coil 107 is disposed inside the gradient magnetic field coil 103, and generates a radio frequency magnetic field upon receiving RF (Radio Frequency) pulses from a transmitting circuit 108. The transmitting circuit 108 supplies the transmitting coil 107 with RF pulses corresponding to a Larmor frequency determined by the type of atom of interest and the magnetic field strength.

[0016] The receiving coil 109 is disposed inside the gradient magnetic field coil 103, and receives magnetic resonance signals (hereinafter referred to as "MR signals" as necessary) emitted from the subject P by excitation with a radio frequency magnetic field. Upon receiving the magnetic resonance signals, the receiving coil 109 outputs the received magnetic resonance signals to the receiving circuitry 110.

[0017] The above-described transmitting coil 107 and receiving coil 109 are merely examples. The coil may be configured by combining one or more of a coil having only a transmitting function, a coil having only a receiving function, or a coil having both a transmitting and receiving function.

[0018] The receiving circuitry 110 detects magnetic resonance signals output from the receiving coil 109 and generates magnetic resonance data based on the detected magnetic resonance signals. Specifically, the receiving circuitry 110 generates magnetic resonance data by digitally converting the magnetic resonance signals output from the receiving coil 109. The receiving circuitry 110 also transmits the generated magnetic resonance data to the sequence control unit 120. The receiving circuitry 110 may be provided on the gantry side that includes the static magnetic field magnet 101, the gradient magnetic field coil 103, etc.

[0019] The sequence control unit 120 performs imaging of the subject P by driving the gradient magnetic field power supply 104, the transmission circuitry 108, and the reception circuitry 110 based on sequence information transmitted from the image processing device 130. Here, the sequence information is information that defines a procedure for performing imaging. The sequence information defines the strength of the current that the gradient magnetic field power supply 104 supplies to the gradient magnetic field coil 103 and the timing of supplying the current, the strength of the RF pulse that the transmission circuitry 108 supplies to the transmission coil 107 and the timing of applying the RF pulse, and the timing of detecting a magnetic resonance signal by the reception circuitry 110. For example, the sequence control unit 120 is an integrated circuit such as an ASIC (Application Specific Integrated Circuit) or an FPGA (Field Programmable Gate Array), or an electronic circuit such as a CPU (Central Processing Unit) or an MPU (Micro Processing Unit).

[0020] The image processing device 130 performs processes such as image generation. In addition, the image processing device 130 performs overall control of the magnetic resonance imaging apparatus 100. The image processing device 130 includes a memory 132, an input device 134, a display 135, and a processing circuitry 150. The processing circuitry 150 includes a control unit 150a, a generation unit 150b, a display control unit 150c, a first calculation unit 150d, a second calculation unit 150e, a third calculation unit 150f, a fourth calculation unit 150g, and a range selection unit 150h.

[0021] The control unit 150a, the generation unit 150b, the display control unit 150c, the first calculation unit 150d, the second calculation unit 150e, the third calculation unit 105f, the fourth calculation unit 150g, and the range selection unit 150h are realized by a processing circuit that executes the corresponding functions. That is, the control unit 150a, the generation unit 150b, the display control unit 150c, the first calculation unit 150d, the second calculation unit 150e, the third calculation unit 105f, the fourth calculation unit 150g, and the range selection unit 150h are realized by the processing circuit 150 having a processor that reads out from the memory 132 a program having a function corresponding to each unit and executes the program.

[0022] Here, the term "processor" refers to, for example, a CPU (Central Processing Unit), a GPU (Graphical Processing Unit), an Application Specific Integrated Circuit (ASIC), a programmable logic device (for example, a Simple Programmable Logic Device (SPLD), a Complex Programmable Logic Device (CPLD), and a Field Programmable Gate Array (FPGA)).

[0023] The control unit 150a transmits sequence information to the sequence control unit 120 and receives magnetic resonance data from the sequence control unit 120. The control unit 150a also stores the received magnetic resonance data in the memory 132.

[0024] The control unit 150a arranges the magnetic resonance data stored in the memory 132 in the k-space. As a result, the memory 132 stores the k-space data.

[0025] The generator 150b generates a magnetic resonance image by performing image reconstruction on the magnetic resonance data acquired from the sequence controller 120. The generator 150b also reads out k-space data from the memory 132 and performs reconstruction processing such as Fourier transform on the read out k-space data to generate a magnetic resonance image.

[0026] The functions of the display control unit 150c, the first calculation unit 150d, the second calculation unit 150e, the third calculation unit 150f, the fourth calculation unit 150g, and the range selection unit 150h will be described later.

[0027] The memory 132 stores magnetic resonance data received by the control unit 150 a, k-space data arranged in k-space, image data generated by the generation unit 150 b, etc. For example, the memory 132 is a semiconductor memory element such as a RAM (Random Access Memory), a flash memory, a hard disk, an optical disk, etc.

[0028] The input device 134 accepts various instructions and information input from an operator. The input device 134 is, for example, a pointing device such as a mouse or a trackball, a selection device such as a mode switch, or an input device such as a keyboard. The display 135, under the control of the control unit 150a, displays a GUI (Graphical User Interface) for accepting input of imaging conditions, images generated by the generation unit 150b, etc. The display 135 is, for example, a display device such as a liquid crystal display.

[0029] The control unit 150a performs overall control of the magnetic resonance imaging apparatus 100, and controls imaging, image generation, image display, etc. For example, the control unit 150a accepts input of imaging conditions (imaging parameters, etc.) on a GUI and generates sequence information according to the accepted imaging conditions. The control unit 150a also transmits the generated sequence information to the sequence control unit 120.

[0030] Next, the background of the embodiment will be briefly described.

[0031] T2 using gradient echo (GRE) in magnetic resonance imaging ＊ In imaging, the phase component of the obtained MRI signal is proportional to the change in the local magnetic field, and this change in the local magnetic field is due to magnetic substances contained in the tissue, such as water and iron. There are several imaging methods that utilize this phase component, including the Quantitative Susceptibility Map (QSM) method, which maps magnetic susceptibility, and Susceptibility Weighted Imaging (SWI), and the Phase-difference-enhanced (PADRE) method, which selectively utilizes tissue-specific phase components to create images.

[0032] Although these technologies all utilize magnetic susceptibility obtained from phase information, from the perspective of clinical application, they cannot be said to have sufficient diagnostic ability, particularly in the early diagnosis of dementia in asymptomatic preclinical Alzheimer's disease.

[0033] In view of this background, a magnetic resonance imaging apparatus 100 according to an embodiment includes a sequence control unit 120, a first calculation unit 150d, a second calculation unit 150e, and a display control unit 150c. The sequence control unit 120 executes a pulse sequence to collect data from an imaging target. The first calculation unit 150d calculates a first map indicating the magnitude of a quantity corresponding to volume magnetic susceptibility based on the data. The second calculation unit 150e calculates a second map indicating feature quantities of voxel values ​​of the first map or a partial region within the first map based on the first map. The display control unit 150c displays an image on a display 135 as a display unit based on the second map.

[0034] The image processing device 130 according to the embodiment includes a first calculator 150d, a second calculator 150e, and a display controller 150c. The first calculator 150d calculates a first map indicating the magnitude of a quantity corresponding to volume magnetic susceptibility based on data collected from an imaging target by executing a pulse sequence. The second calculator 150e calculates a second map indicating feature quantities of voxel values ​​of the first map or a partial region within the first map based on the first map. The display controller 150c analyzes the second map and displays an image on a display 135 serving as a display unit.

[0035] The processing performed by the magnetic resonance imaging apparatus 100 according to the embodiment will be described below with reference to Figures 2 to 23. Figure 2 is a flowchart showing the flow of processing performed by the magnetic resonance imaging apparatus 100 according to the embodiment. Figure 2 is a flowchart showing an overall image of the processing, and details of the processing in steps S100 and S200 are explained in more detail in Figure 3, and details of the processing in steps S300 and S400 are explained in more detail in Figure 9.

[0036] First, in step S100, the sequence control unit 120 executes a pulse sequence to collect data from the imaging target. Here, the pulse sequence executed by the sequence control unit 120 in step S100 is, for example, a sequence that generates multiple echoes by one excitation using a high-frequency magnetic field. The upper limit of the number of multiple echoes that can be generated is determined by the signal-to-noise ratio of the echo signals, and is determined by the static magnetic field strength, gradient magnetic field performance, or T2 ＊ The number of echoes also depends on the type of weighted imaging. In general clinical applications, it is about 16 echoes. As an example, as shown in FIG. 3, in step S110, the sequence control unit 120 executes a pulse sequence of the Multi-GRE method, which is a pulse sequence of the GRE (Gradient Echo) method executed for a plurality of TEs (Echo Times), and obtains T2 ＊Enhanced imaging is performed. Hereinafter, the imaging coordinate system executed in step S110 is referred to as a first imaging coordinate system. The sequence controller 120 executes a pulse sequence in which the echo times are four different echo times, for example, TE1, TE2, TE3, and TE4, and collects data of four echoes. Here, typically, the sequence controller 120 executes the pulse sequence so that the time interval (ΔTE) between each echo is constant.

[0037] Fig. 4 shows an example of a phase image obtained by the pulse sequence executed by the sequence controller 120. Fig. 4 shows each phase image in a predetermined slice cross section in a two-dimensional multi-slice image or a three-dimensional phase image obtained when the sequence controller 120 executes a pulse sequence for four echoes, TE1, TE2, TE3, and TE4, in step S100.

[0038] Subsequently, in step S200, the first calculator 150d calculates a first map indicating the magnitude of a quantity corresponding to the volume magnetic susceptibility based on the magnetic resonance data collected in step S100. Details of the process of step S200 are shown in FIG.

[0039] First, in step S210, the first calculator 150d calculates a plurality of T2 ＊ A phase unwrapping process is performed on the phase image.

[0040] Next, in step S215, the first calculator 150d applies a high-pass filter to each phase image subjected to the phase unwrapping process in step S210. This makes it possible to remove inhomogeneity in the static magnetic field caused by biological tissue in the static magnetic field magnet. Instead of (or in addition to) the high-pass filter, a background removal filter may be applied. This makes it possible to remove low-frequency noise occurring in the background of the phase image. This step may be omitted if the inhomogeneity of the static magnetic field is sufficiently guaranteed or if there is no noise occurring in the background of the phase image.

[0041] Next, in step S220, to reduce phase noise, the first calculation unit 150d applies a Wiener filter to the phase image for each echo to which the high-pass filter was applied in step S220. Here, the Wiener filter is an adaptive filter that has a low filtering effect in areas where signal values ​​change significantly, i.e., areas representing structure, but has a high filtering effect on small noise components. When the diagnostic performance was compared between the case where a Gaussian filter, which performs noise reduction uniformly over the entire image, and the case where a Wiener filter was used, better results were obtained when using the Wiener filter than when using the Gaussian filter.

[0042] Subsequently, in step S225, the first calculator 150d calculates, for each voxel, a value (Δθ / ΔTE) obtained by dividing the value of the phase change (Δθ) for each echo by the value of the change in TE (ΔTE) for each echo. The first calculator generates a phase slope value (Δθ / ΔTE) map based on the value Δθ / ΔTE calculated in step S225. The phase slope value is also synonymous with the terms "phase gradient," "slope value," and "slope value."

[0043] Details of this process will be described assuming that the sequence control unit 120 executes a pulse sequence using the Multi-GRE method in step S110 and collects data for four echoes TE1, TE2, TE3, and TE4.

[0044] In this case, the phase value θ, which is the voxel value for each TE, is expressed by the following equation (1).

[0045]

[0046] where γ is the gyromagnetic ratio, B is the magnitude of the static magnetic field, and χ is the volume magnetic susceptibility. Strictly speaking, χ is a three-dimensional tensor quantity, but for clinical application, the expression is simplified to volume magnetic susceptibility.

[0047] Here, from equation (1), the phase difference Δθ between each echo is given by the following equation (2) when the difference in TE between each echo is ΔTE.

[0048]

[0049] From equation (2), the phase slope value (Δθ / ΔTE) is the value obtained by dividing the phase change value (Δθ) between each echo by the TE change value (ΔTE) between each echo, and is expressed by the following equation (3).

[0050]

[0051] Here, the phase slope value (Δθ / ΔTE) for TE has units of rad / sec and is a physical quantity (units [rad / sec]) having the dimension of angular frequency ω. Since the equation ω = 2πf (f is frequency) holds, the phase slope value can also be expressed in units of [Hz]. From equation (3), the phase slope value (Δθ / ΔTE) for TE is proportional to the volume magnetic susceptibility χ, and therefore the phase slope value (Δθ / ΔTE) for TE can be interpreted as a quantity that roughly represents the volume magnetic susceptibility χ.

[0052] Fig. 5 shows a plot of the phase value θ of the voxel indicated by the arrow in Fig. 4 against the echo time TE. Points 10a, 10b, 10c, and 10d are data points plotting the phase value θ of a given voxel, for example, the given voxel indicated by the arrow in Fig. 4, in the acquisitions of TE1, TE2, TE3, and TE4, respectively.

[0053] The first calculator 150d calculates a value (Δθ / ΔTE) by dividing the value of the phase change (Δθ) between each echo by the value of the change in TE (ΔTE) between each echo. As an example, in the case of two echoes in the Multi-GRE method, the first calculator 150d calculates the value of the phase change Δθ based on the difference between the phase of point 10a and the phase of point 10b, calculates the value of the change in TE ΔTE between each echo based on the difference between the TE of point 10a and the TE of point 10b, and calculates Δθ / ΔTE based on these. As another example, in the case of three or more echoes in the Multi-GRE method, the first calculator 150d calculates the value of Δθ / ΔTE based on the data at points 10a and 10b, calculates the value of Δθ / ΔTE based on the data at points 10b and 10c, calculates the value of Δθ / ΔTE based on the data at points 10c and 10d, and calculates the phase slope value (Δθ / ΔTE) by averaging these values.As another example, the first calculator 150d calculates a straight line 11 that best fits the data points at points 10a, 10b, 10c, and 10d using, for example, the least squares method, and calculates the phase slope value (Δθ / ΔTE) from the slope of the straight line 11. 6 shows a phase slope value (Δθ / ΔTE) map 31 in which the phase slope value (Δθ / ΔTE) calculated for each voxel is mapped as a voxel based on images 30a, 30b, 30c, and 30d, which are phase images having TE values ​​of TE1, TE2, TE3, and TE4, respectively. The first calculator generates the first map based on the phase slope value Δθ / ΔTE calculated in step S225, through steps described below.

[0054] The unit of the phase slope value (Δθ / ΔTE) obtained by dividing the value of the phase change between each echo (Δθ) by the value of the change in TE between each echo (ΔTE) will be explained below. When the units are converted so that Equation (3) can be expressed using the values ​​of the data collected by the magnetic resonance imaging apparatus 100, the phase value (rad) in the MRI phase image takes a value in the range of 2π between -π and π, and the discretization is performed, for example, with 12 bits. Therefore, the change in pixel value of the discretized phase (Δθ bｉｔ ) takes integer values ​​from 0 to 4095, with a width of 4096. Therefore, in equation (3), the volume magnetic susceptibility per voxel is expressed as χvoxel Then, the following equation (4) holds.

[0055]

[0056] Here, in the case of hydrogen (i.e., normal MRI), γ / 2π=42.57×10 6 [Hz / T], the following equation (5) holds true.

[0057]

[0058] 1ml (1000mm 3 ) the volume magnetic susceptibility χ corresponding to the amount of magnetic material per mｌ is the volume of a voxel, V vｏｌ [mm 3 ], it is expressed by the following equation (6).

[0059]

[0060] When equation (6) is substituted into equation (5), the following equation (7) is established.

[0061]

[0062] That is, when the first calculator 150d evaluates the value of Δθ / ΔTE in units of the left side of equation (7), the first calculator 150d uses equation (7) and the volume V of the voxel as ｖｏｘｅｌ By using the above, the volume magnetic susceptibility χ corresponding to the amount of magnetic material per 1 ml can be calculated. ｍｌ The value of can be calculated.

[0063] Subsequently, in step S230, the first calculation unit 150d performs a predetermined linear transformation on the phase slope value (Δθ / ΔTE) for TE calculated in step S225.

[0064] Since the soft tissue of a living body is mainly composed of water and protein, it exhibits diamagnetism, and its volume magnetic susceptibility χ is approximately -9 ppm. Therefore, as shown in curve 20 of Figure 7, the histogram of the phase slope value (Δθ / ΔTE) of the voxels of a three-dimensional structure, with the horizontal axis representing the phase slope value (Δθ / ΔTE) and the vertical axis representing the number of voxels that have these values, generally exhibits a normal distribution, and the peak value of this normal distribution takes a positive value on the horizontal axis. Depending on the MRI device used, after the living body is placed in the static magnetic field and gradient magnetic field coil, the magnitude of the static magnetic field is adjusted to match the magnetic susceptibility of the soft tissue of the living body, and the peak value of the histogram of the phase slope value of curve 20 of Figure 7 may actually be approximately 0.

[0065] Here, for example, iron deposition associated with the accumulation of amyloid β protein (hereinafter referred to as "Aβ protein") in Alzheimer's disease and iron ions that increase with aging exhibit paramagnetic properties, so that when a lesion is present, the curve 20 shifts in the direction of the arrow 23, which is the paramagnetic side, due to the increase in iron deposition and iron ions. In other words, when a lesion is present, the pixel value shifts in the negative direction.

[0066] However, when displaying a phase slope value (Δθ / ΔTE) map, it is easier for the user to understand if the image shows pixel values ​​shifting in the positive direction in the presence of a lesion, as in a PET image that visualizes Aβ protein. Therefore, in the magnetic resonance imaging apparatus 100 according to the embodiment, the first calculator 150d performs a linear transformation on the Δθ / ΔTE value calculated in step S225 to generate data that is easier for the user to understand.

[0067] Specifically, the first calculator 150d performs a linear transformation on Δθ / ΔTE calculated in step S225 by multiplying the value on the horizontal axis by −1 and inverting the horizontal axis, thereby creating the data shown by curve 21 in FIG. 8. Because the sign is inverted at this stage, when a lesion occurs and the curve 21 shifts toward the paramagnetic side, the pixel values ​​of curve 21 shift in the positive direction. That is, when a lesion occurs, the pixel values ​​of curve 21 shift in the positive direction.

[0068] Here, curve 21 peaks when the signal value is negative under normal conditions where no lesion is present. However, it is easier for the user to understand if the signal value peaks when the signal value is positive. Therefore, the first calculator 150d further performs a linear transformation on curve 21, adding a constant value to the horizontal axis value so that the horizontal axis value at the peak of the signal value becomes positive, thereby generating curve 22. That is, the first calculator 150d inverts the horizontal axis of the histogram of the phase slope value (Δθ / ΔTE) calculated in step S225, and then performs a linear transformation by adding a constant value to the horizontal axis value to generate curve 22 and a phase slope value (Δθ / ΔTE) map. In other words, the first calculator 150d performs a transformation to convert the phase slope value (Δθ / ΔTE) calculated in step S225 to −Δθ / ΔTE+α. As a result, the peak value of the normally distributed histogram will be a positive value when normal, and the peak value of the histogram will be a large positive value when diseased, resulting in pixel values ​​being large positive values, making it possible to generate data that is easy for the user to understand.

[0069] In this way, the first calculator 150d calculates a phase slope (Δθ / ΔTE) map by linearly transforming the value obtained by dividing the phase change between echoes by the change in TE between echoes so that the converted value becomes positive when the imaging target contains a paramagnetic component. This allows the creation of data that is intuitively easy for the user to understand.

[0070] 3 , in step S240, the first calculator 150d generates an echo intensity image in the first imaging coordinate system based on the data obtained from the pulse sequence executed in step S110. Subsequently, in step S245, the sequence controller 120 performs T1-weighted imaging in a second imaging coordinate system different from the first imaging coordinate system in which imaging was executed in step S110.

[0071] Next, in step S250, the first calculator 150d generates a T1-weighted image after deformation and registration based on the T1-weighted imaging performed in the second imaging coordinate system in step S245. This deformation and registration makes the T1-weighted image equivalent to an image captured in the first imaging coordinate system. The T1-weighted image generated in step S250 is used to generate a vascular component image mask in step S255, a cerebrospinal fluid (CSF) mask in step S260, and create a standard parcel atlas in step S420.

[0072] Subsequently, in step S255, the first calculator 150d generates a vascular component image mask based on the phase image that has been subjected to the unwrapping process and that has been generated in step S210, and the T1-weighted image that has been subjected to the process of step S250. Furthermore, in step S260, the first calculator 150d generates a CSF mask based on the T1-weighted image that has been subjected to the process of step S250.

[0073] In step S270, the first calculator 150d generates a first map, a slope map, based on the phase slope value (Δθ / ΔTE) map obtained after linear transformation in step S230, the blood vessel component image mask generated in step S255, and the CSF mask generated in step S270.

[0074] That is, the first calculation unit 150d removes the vascular components and CSF components from the phase slope value (Δθ / ΔTE) map after linear transformation using the vascular component image mask and CSF mask, and creates a first map as the brain parenchyma to be the target of feature analysis in S300.

[0075] As described above, in step S200, the first calculator 150d creates a phase slope value (Δθ / ΔTE) map based on the data collected in step S100, and calculates a first map indicating the magnitude of a quantity corresponding to the volume magnetic susceptibility of only the brain parenchyma by removing the vascular components and CSF components. However, the embodiment is not limited to the case where the first map indicating the magnitude of a quantity corresponding to the volume magnetic susceptibility is created by creating a phase slope value (Δθ / ΔTE) map. The first calculator 150d may create the first map indicating the magnitude of a quantity corresponding to the volume magnetic susceptibility of only the brain parenchyma by removing the vascular components and CSF components from an image, such as an image obtained by a quantitative susceptibility map (QSM) method or a susceptibility weighted imaging (SWI) method, as the phase slope value (Δθ / ΔTE) map.

[0076] Next, in step S300, the second calculator 150e calculates a second map indicating feature quantities of the signal values ​​of the first map based on the first map. The second map is used to generate a third map and a fourth map, for example, in steps S500 and S700, which will be described later. FIG. 9 shows the details of the processing of step S300.

[0077] Subsequently, in step S370, the second calculation unit 150e calculates texture features based on the first map calculated in step S270, and generates a second map representing the texture features. Here, the texture features may be, for example, the mean value m and variance σ of multiple voxel values ​​included in a specific three-dimensional structure of the first map. 2 , standard deviation σ, entropy e, Z-value, skewness or kurtosis, etc. Here, entropy e is, for example, Shannon's information entropy. Furthermore, the Z-value is calculated using a population of normal individuals, with its mean value μ and standard deviation σ n Then, Z = (m - μ) / σ n means the value given by

[0078] A particularly promising texture feature is m / e / σ, a combination of texture values, where m, e, and σ represent the mean value, entropy, and standard deviation of multiple voxel values ​​contained in a specific three-dimensional structure of the first map, respectively.

[0079] If the texture feature amount in step S370 includes entropy, kurtosis, skewness, or the like, the process in step S300 proceeds to step S350, in which the second calculator 150e calculates a histogram of the first map (hereinafter referred to as the slope map) from a plurality of voxel values ​​included in a specific three-dimensional region of the first map.

[0080] In step S355, the range selector 150h selects a range of phase slope values. That is, the range selector 150h determines the selection range of the phase slope values ​​based on the phase change value and the TE (Echo Time) value of the first map. Because the MRI signal (MRI magnetic field strength) output from the MRI apparatus is not fixed and has no standard specification, the histogram calculated in step S350 is also not fixed. Therefore, since the count numbers of the histogram calculated in step S350 generally have a shape similar to a normal distribution, the minimum and maximum values ​​at which the count numbers can be considered to be approximately 0 are selected. The range from the selected minimum value to the maximum value becomes the "selection range" of the phase slope value. The minimum and maximum values ​​for selecting the range of the phase slope value may be determined according to a predetermined standard or may be determined arbitrarily by manual input, etc. As will be described later, in order to capture the original spatial distribution, a range wider than that of the limited QSM method may be selected.

[0081] FIG. 10 shows a histogram 70 of the first map of voxels generated in step S320, excluding masked vessel or CSF regions.

[0082] Here, we explain the differences between the QSM method and our method for histograms. The QSM method analyzes the minute magnetic fields generated by magnetic materials in their vicinity based on electromagnetic field theory. It is a magnetic susceptibility mapping method that utilizes the principle that the effect of this minute magnetic field (vector field) on the phase component of an MRI signal (complex number) is proportional to the magnetic susceptibility. In this case, the magnetic susceptibility analysis targets a vector field, resulting in a directional tensor magnetic susceptibility. Because the QSM method performs rigorous magnetic susceptibility analysis based on minute magnetic fields, it can be said that the analysis is limited in terms of the MRI phase component. In other words, because the phase slope value is limited in the QSM method, the entropy and variance do not capture the true spatial distribution. Furthermore, because the QSM method optimizes dipole inversion according to the characteristics of each individual image, variations in the images, such as outliers, can result in significant variations in the final QSM images. In this respect, the phase slope value without dipole inversion is less likely to vary between individuals and can be said to be stable.

[0083] On the other hand, the analysis target of this method is a quantity corresponding to volume magnetic susceptibility. The essential difference from QSM is that the analysis can be performed by freely selecting a range of phase slope values ​​within the distribution of phase slope values ​​in the phase direction (frequency direction) (see Figure 11). Figure 11 shows an example in which the range of phase slope values ​​is selected from a minimum of 0 [rad / sec] to a maximum of 25 [rad / sec]. The selected range may be selected to facilitate diagnosis depending on the diagnostic objective, white matter, gray matter, etc. In other words, the entire range of phase slope values ​​may be selected, or a portion of the phase slope values ​​may be selected as long as it is wider than the QSM method. By selecting a wider range of phase slope values ​​than the QSM method, the original spatial distribution of texture features such as entropy and variance can be captured.

[0084] In step S360, the second calculator 150e calculates statistics of the first map, which is a slope map, based on the histogram 70 calculated in step S350 and the selected range of phase slope values ​​selected in step S355. Based on this, the second calculator 150e calculates texture features in step S370, which will be described later. That is, if the texture features include entropy e, kurtosis, skewness, etc., the features calculated in step S370 are calculated from the histogram of voxels created in step S350, excluding masked blood vessel or CSF regions.

[0085] In this way, in step S300, the second calculation unit 150e calculates the second map representing texture features such as the mean, standard deviation, and entropy, and combination features of texture features such as m / e / σ.

[0086] The effectiveness of these texture features in identifying lesions will be shown by the results of Parcel analysis performed in step S400, which will be described later. The effectiveness of the texture features listed here will be explained after explaining the processes in steps S400 and S500.

[0087] Next, in step S400, the third calculation unit 150f calculates a third map by aggregating the second map for each functional unit of the imaging target, for example, known as a parcel. A more detailed processing flow of step S400 will be described below with reference to FIG. 9 .

[0088] First, the standard parcel atlas will be described. A parcel refers to a functional unit of an imaging target, and an atlas refers to a map that represents location information. When the imaging target is the brain, for example, there is a standard parcel atlas called AAL (Automated Anatomical Labeling), which is a global standard that divides brain tissue into 116 functional regions. In addition to AAL, there are standard parcel atlases for gray matter and white matter. In step S420, at least one standard parcel atlas is selected according to the purpose of analysis. Furthermore, as necessary, characteristic parcels are selected from each standard parcel atlas. That is, the third calculation unit 150f acquires each parcel, i.e., functional unit, in association with its location information.

[0089] Next, in step S430, the third calculation unit 150f performs deformation and registration of the standard parcel atlas selected in step S420 on the T1-weighted image after deformation and registration generated in step S250, and creates a parcel region for analysis and display.

[0090] Subsequently, in step S450, the third calculation unit 150f calculates a third map by aggregating the second map for each parcel based on the second map generated in step S370 and the parcel region for analysis and display created in step S430. In this manner, in step S400, the third calculation unit 150f calculates a third map by aggregating the second map for each functional unit of the imaging object. As an example, the third calculation unit 150f calculates an average value of the second map values, which are texture features, for each functional unit of the imaging object, and sets the calculated average value as the value in the third map for all voxels in that functional unit. As an example, if m / e / σ is selected as the texture feature in the second map in step S370, the third calculation unit 150f calculates a value of m / e / σ, which is a combination of texture features, for each parcel, and sets the calculated value as the value in the third map for the voxels included in that parcel. In other words, the third calculator 150f generates the third map so that, for example, all voxel values ​​belonging to the same parcel are the same.

[0091] Returning to FIG. 2, subsequently, in step S500, the display control unit 150c causes the display 135 serving as the display unit to display the third map generated in step S400.

[0092] As an example, the display control unit 150c refers to the color table 75 shown in FIG. 12 and causes the display 135 serving as a display unit to display the third map generated in step S400 in color for each functional unit, as shown in FIG. 13 . As an example, if m / e / σ is selected as the texture feature in the second map in step S370, the display control unit 150c adjusts the color tone of the color table 75 so that, for example, a high m / e / σ value results in a reddish color and a low m / e / σ value results in a bluish color. Then, the display control unit 150c causes the display 135 serving as a display unit to display the third map generated in step S400 in color. Here, the third map is a map in which a texture feature is assigned to each functional unit. In FIG. 13 , for example, regions 41, 42, and 43 belong to different parcels, and therefore have different values ​​according to the texture feature of each parcel and are displayed in different colors.

[0093] As another example, the display control unit 150c may superimpose the third map generated in step S400 on the display 135 as a display unit, along with a T1-weighted image or a vascular image. For example, FIG. 14 shows an example of such a superimposed image. In FIG. 14, for example, regions 80 and 81 are regions of the third map, in which texture features are assigned to each functional unit and displayed in color. In contrast, region 82 is a region in which T1-weighted images or vascular images are primarily displayed and displayed in monochrome. That is, the display control unit 150c displays the texture features assigned to each functional unit in color based on the third map generated in step S400, while displaying regions for which texture features have not been calculated in monochrome using T1-weighted images or vascular images. This allows the user to understand the positional relationship between the brain tissue structure, vascular structure, and parcels, thereby generating an image that is easy for the user to understand and enabling accurate diagnosis.

[0094] When performing the superimposed display, the display control unit 150c may use a toggle display function, i.e., a function for switching between two screens with a single button, to switch between the portion displaying the third map and the portion displaying the T1-weighted image or vascular image, making it easier to compare the two images. Furthermore, a function for adjusting the transparency of the parcel side may be used to increase the transparency, making it easier to compare the positional relationship between the brain tissue structure or vascular structure and the parcel.

[0095] The display control unit 150c may also receive a mouse operation from the user via the input device 134 and display the name of the parcel at the position designated by the user on the display 135. In the standard parcel atlas in step S420, a number is assigned to each parcel, and the name of the parcel can be displayed based on the number.

[0096] Furthermore, when the display control unit 150c receives a mouse operation from the user via the input device 134, it may display the texture feature quantity at the position designated by the user on the display 135 as a display unit. As an example, the display control unit 150c may display the m / e / σ values, average value m, etc. at the position pointed to by the mouse cursor or the like on the display 135 as a display unit. Furthermore, the display control unit 150c may also display the standard values ​​and Z values ​​of these numerical values ​​on the display 135.

[0097] Next, the selection of texture features in step S370 will be explained again. In order to consider what texture features should be adopted, a first map was created for 45 cases, and slope map analysis was performed. These 45 cases were classified into cases without Aβ protein (PET-negative) and cases with Aβ protein (PET-positive) based on the image findings using amyloid PET images. As a result, there were 18 cases that were PET-negative and 27 cases that were PET-positive.

[0098] In this study, the DK atlas, which is a type of standard parcel atlas for the cortex, was used in step S420. The number of parcels in the DK atlas is 24.

[0099] Furthermore, the texture feature values ​​calculated in step S370 are the mean value m, the standard deviation σ, and the entropy e. These texture feature values ​​are calculated in step S360 from the histogram of all voxels excluding the blood vessels and CSF regions masked by the tissue within the parcel, calculated in step S350. The texture feature values ​​are not limited to the mean value m, the standard deviation σ, and the entropy e described above. Other possible values ​​include a fractal dimension (value), a gray-level co-occurrence matrix (GLCM), a gray-level size zone matrix (GLSZM), a neighborhood gray-tone difference matrix (NGTDM), local binary patterns (LBP), a Laplacian distribution, a run length matrix, and frequency filtering (such as a wavelet transform or a Fourier transform).

[0100] Regarding the texture feature quantities of the mean, standard deviation, and entropy described above, we analyzed how the texture feature quantities changed between PET-negative cases and PET-positive cases.

[0101] As a result, the following results were obtained. Note that the comparison between the two groups of PET-positive cases and PET-negative cases was performed by t-test for each of the 24 parcels. The following p-values ​​are the average p-values ​​for each parcel.

[0102] When the mean slope map m was used as the texture feature, a comparison of the values ​​between the two groups showed that the PET-negative cases were lower than the PET-positive cases, with a p-value of p<0.15. The mean value of the slope map histogram tended to be higher on the PET-positive side. This is thought to be due to an increase in paramagnetic components caused by an increase in iron components in the PET-positive cases.

[0103] Furthermore, when the standard deviation σ of the slope map was used as the texture feature, a comparison of the values ​​between the two groups showed that PET-negative cases were greater than PET-positive cases, with a p-value of p<0.1. A tendency was observed for the spread of the slope map histogram to be lower on the PET-positive side. This is thought to be because, in the absence of proteins such as Aβ to which iron components are bound, the slope map values ​​are approximately normally distributed, but as the amount of iron-bound proteins increases, the number of voxels with corresponding slope map values ​​increases, and the area near the peak of the normal distribution of the histogram of slope map values ​​increases, resulting in a relatively low standard deviation representing the spread of the normal distribution.

[0104] Furthermore, when the entropy e of the slope map was used as a texture feature, a comparison of the values ​​between the two groups showed that the PET-negative cases were lower than the PET-positive cases, with a p-value of p<0.15. A tendency for the randomness of the slope map histogram to be higher on the PET-positive side was observed. This is because, in the absence of highly specific proteins such as Aβ bound to iron components, the randomness of the spatial distribution of the slope map values ​​is large. However, as the amount of proteins such as Aβ bound to iron components increases, the number of voxels with corresponding slope map values ​​increases, and the specific granularity increases, which is thought to be due to the reduced randomness of the spatial distribution of the slope map values.

[0105] Based on the above considerations, for example, in step S370, the second calculation unit 150e employs m / e / σ, obtained by dividing the mean m, which increases in PET-positive cases, by the entropy e, which decreases in PET-positive cases, and then dividing this by the standard deviation σ, which decreases in PET-positive cases, as the texture feature. This results in PET-negative cases being less than PET-positive cases, with a p-value of p<0.05. In other words, it can be seen that a value obtained by dividing the mean m of the voxel values ​​in the three-dimensional structure of the first map by the entropy e and then by the standard deviation σ is a promising option. A high value of this feature can be considered to indicate an increase in the amount of iron components due to amyloid accumulation or aging.

[0106] Returning to FIG. 2 , the processing of steps S600 and S700 will be described. Steps S400 and S500 described a case where a third map is generated and displayed by aggregating the second map by functional unit. Steps S600 and S700 described a case where a fourth map is generated by aggregating the second map by voxel. However, since simply calculating the value for each voxel can result in discontinuities, the fourth map is generated by performing processing such as averaging or smoothing filter processing.

[0107] In step S600, the fourth calculator 150g applies an image filter to the second map to generate a fourth map for each voxel. That is, since the second map, which is a slope map, tends to have discrete values ​​as is, the fourth calculator 150g applies a low-pass image filter, such as a Gaussian filter, to the second map to remove graininess from the second map and generate an image similar to a PET image, thereby generating a fourth map for each voxel.

[0108] When a gray level (also called gray scale or shading) image, monochrome color display, or color display is performed, the fourth calculation unit 150g creates a histogram 90 of a specific region of the slope map, as shown in Fig. 15. Based on the histogram 90, the fourth calculation unit 150g performs settings such that the average value of the slope map is positioned in the center of a display table, such as a color table, used for display. Note that the specific region in this case is a standard region or a region where Aβ protein accumulation is likely to occur.

[0109] In addition, when setting the above-mentioned display table, it may be changed appropriately depending on the purpose of the test, such as asymptomatic preclinical Alzheimer's disease (preclinical AD) mainly tested in medical checkups, or mild dementia (MCI) or Alzheimer's disease (AD) patients mainly tested in hospitals.

[0110] An example of such a situation is shown in FIG. 16. In FIG. 16, consider the case where the histogram of a specific region of the slope map is given as histogram 90. FIG. 16 shows the case where the range indicated by region 95b is selected and the display table indicated by region 95a is matched, and the case where the range indicated by region 96b is selected and the display table indicated by region 96a is matched. Normally, the selected range and the display table are matched, as in the relationship between region 96b and region 96a, but it is not necessary for the selected range and the display table to match. If they do not match, the colors in the mismatched range will not be displayed. This is useful when you want to narrow the range of colors to display.

[0111] First, in step S360, the second calculation unit 150e sets a range of voxel values ​​to be used to calculate texture features according to the diagnostic purpose of the imaging. For example, if the diagnostic purpose is a medical checkup or other examination, the second calculation unit 150e determines that voxel values ​​in the histogram 90 whose voxel value range is in region 95b are to be used to calculate texture features. Medical checkups and other examinations include many normal patients who do not have Alzheimer's disease. For screening tests, it is optimal to use data from region 95b, which is a wide region, for the histogram created to calculate texture feature values. Furthermore, if the diagnostic purpose is a hospital examination or other examination, the second calculation unit 150e determines that voxel values ​​in the histogram 90 whose voxel value range is in region 96b are to be used to calculate texture features. For hospital examinations of Alzheimer's disease patients, it is optimal to use data from region 96b, which is a relatively narrow region, for texture feature values.

[0112] In step S370, the second calculator 150e determines texture features based on a range of voxel values ​​determined according to each diagnostic purpose. That is, the texture features calculated in step S370 may be determined based on all voxel values ​​included in a partial region (spatial target range) in the first map, as already described in step S350, or may be calculated based on texture features of voxel values ​​within a certain range of values, i.e., a specific numerical range on the slope map.

[0113] Furthermore, independently of changing the range of values ​​used to calculate the texture feature amount according to the diagnostic purpose, in step S500 or step S700, the display control unit 150c may change the color table used when displaying the third map or the fourth map on the display 135 according to the diagnostic purpose of imaging, in order to perform effective diagnosis using the texture feature amount map.

[0114] As an example, for an examination in a medical checkup or the like, the display control unit 150c displays the third map or the fourth map using the color table shown in area 95a. Furthermore, for use in a diagnosis or examination in a hospital, the display control unit 150c displays the third map or the fourth map using the color table shown in area 96a. In this way, by changing the range of the color map in the display depending on the diagnostic purpose, it is possible to visualize an image that is optimal for diagnosis.

[0115] In this case, the color table set by the display control unit 150c in step S500 or step S700 may be determined based on the range of the histogram 90 used to determine the texture feature in step S360 or step S370. For example, in the case of an examination in a medical checkup or the like, the display control unit 150c may determine the color table 95a based on region 95b. In addition, in the case of a diagnosis or examination in a hospital, the display control unit 150c may determine the color table shown in region 96a based on region 96b.

[0116] In the above example, the case where the color table range is calculated based on the texture feature calculation range has been described, but the embodiment is not limited to this. As an example, the texture feature calculation range may not be limited to a specific region, and all voxels may be subject to texture feature calculation, and different color tables may be used to visualize the image depending on the diagnostic purpose, etc.

[0117] As another example of setting the display table, the display control unit 150c may display the fourth map on the display 135 as a display unit by setting the optimal display table for each region, such as the gray matter region or the white matter region.

[0118] An example of such a situation is shown in Fig. 17. In Fig. 17, consider a case where a histogram of a specific region of a slope map is given as histogram 90. Fig. 17 shows a case where the range shown in region 97b is selected and the display table shown in region 97a is fitted, and a case where the range shown in region 98b is selected and the display table shown in region 98a is fitted.

[0119] First, in step S360, the second calculator 150e sets a range of voxel values ​​to be used to calculate texture features according to the imaging region. For example, if the imaging region is a white matter region, the second calculator 150e determines that voxel values ​​in the range of voxel values ​​in region 97b in the histogram 90 are to be used to calculate texture features. Furthermore, if the imaging region is a gray matter region, the second calculator 150e determines that voxel values ​​in the range of voxel values ​​in region 98b in the histogram 90 are to be used to calculate texture features.

[0120] In step S370, the second calculation unit 150e determines texture feature amounts based on the range of voxel values ​​determined for each imaging region.

[0121] That is, the texture feature calculated in step S370 may be determined based on all voxel values ​​included in a partial region (spatial target range) in the first map, as already described in step S350, or may be calculated based on the texture feature of voxel values ​​whose voxel values ​​are within a certain range of values, i.e., a specific numerical range on the slope map.

[0122] Furthermore, independently of changing the range of values ​​used to calculate the texture feature depending on the imaging area, in step S500 or step S700, the display control unit 150c may change the color table used when displaying the third map or the fourth map on the display 135 depending on the diagnostic purpose of imaging, in order to perform effective diagnosis using the texture feature map.

[0123] As an example, when the imaging region is white matter, the display controller 150c displays the third map or the fourth map using the color table shown in region 97a. Furthermore, when the imaging region is gray matter, the display controller 150c displays the third map or the fourth map using the color table shown in region 98a. In this way, by changing the range of the color map displayed depending on the diagnostic purpose, it is possible to visualize an image optimal for diagnosis.

[0124] In this case, the color table set by the display control unit 150c in step S500 or step S700 may be determined based on the range of the histogram 90 used to determine the texture feature in step S360 or step S370. For example, when the imaging region is white matter, the display control unit 150c may determine the color table indicated by region 97a based on region 97b. Furthermore, when the imaging region is gray matter, the display control unit 150c may determine the color table indicated by region 98a based on region 98b.

[0125] In the above example, the case where the color table range is calculated based on the texture feature calculation range has been described, but the embodiment is not limited to this. As an example, the texture feature calculation range may not be limited to a specific region, and all voxels may be subject to texture feature calculation, and the image may be visualized using different color tables depending on the imaging region, etc.

[0126] Returning to the description of the example of the method of aggregating the second map on a voxel basis to correct the fourth map, the method of aggregating the second map on a voxel basis to correct the fourth map is not limited to the above-described method. As an example, the fourth calculation unit 150g may calculate a central voxel value based on texture features of voxel values ​​of the first map within a predetermined range centered on the central voxel, and then move the central voxel by one voxel to generate a fourth map for each voxel. As an example, the fourth calculation unit 150g may divide the average value m of the voxel values ​​of the first map within a predetermined range centered on the central voxel by the entropy e within the predetermined range of voxel values ​​of the second map, and then divide the result by the standard deviation σ of the voxel values ​​of the second map within the predetermined range to calculate the fourth map at the central voxel. Here, the predetermined range may be, for example, a rectangular parallelepiped or cubic shape.

[0127] 18 , the fourth calculator 150g calculates the average value m of the voxel values ​​of the second map, the entropy e of the voxel values ​​of the first map, and the standard deviation σ of the voxel values ​​of the first map within a predetermined range for a 5×5×5 cube or cubic VOI (Volume of Interest) region 41 of 125 voxels, and then divides the average value m of the voxel values ​​of the first map for the 125 voxel region 41 by the entropy e of the voxel values ​​of the first map for the 125 voxel region 41, and calculates the value m / e / σ by dividing the average value m of the voxel values ​​of the first map for the 125 voxel region 41 by the standard deviation σ of the voxel values ​​of the first map for the 125 voxel region 41 a. The fourth calculator 150g acquires the calculated value of m / e / σ as the value of the fourth map for the central voxel 40a.

[0128] The fourth calculation unit 150g calculates the fourth map by shifting the central voxel 40a by one voxel at a time over the entire three-dimensional region, i.e., the entire first map. The fourth calculation unit 150g may display the fourth map calculated in step S700 on the display 135 as a display unit, or may further apply a Gaussian filter to the fourth map calculated in this step.

[0129] Note that when the fourth map is generated from the first map by the above-described process, the region in which the fourth map is generated is smaller than the region in which the texture feature values ​​are calculated and the predetermined range described above. For example, as shown in FIG. 19 , consider a case in which texture feature values ​​are calculated in a rectangle P'QRS' indicated by a boundary 45, and the fourth calculator 150g calculates the fourth map using the texture feature data of this region. Here, the central voxel 40a and the integral range of the texture feature for the central voxel 40a are indicated as region 41a. In this case, the range in which the fourth calculator 150g can generate the fourth map is the range of the rectangle PQRS indicated by a boundary 46a, in which texture feature data exists in the predetermined range corresponding to the central voxel.

[0130] In the above example, the predetermined range VOI over which the fourth calculator 150g integrates the texture feature when generating the fourth map is described as a rectangular range, but the embodiment is not limited to this. The predetermined range VOI over which the fourth calculator 150g integrates the texture feature when generating the fourth map may be a spherical range, for example, as shown in FIG. 20 . In other words, the predetermined range is a range with a three-dimensional spherical shape. This allows the integration range of the texture feature to be within a certain distance from the central voxel, thereby improving the accuracy of the generated fourth map.

[0131] The upper part of Fig. 20 shows a three-dimensional representation of a specific central voxel 40b and a region 41b, which is a predetermined range for integrating texture features for the central voxel 40b. The lower part of Fig. 20 shows a two-dimensional representation of the central voxel 40b and a region 41b, which is a predetermined range for integrating texture features for the central voxel 40b, by displaying a central cross section of the VOI of the region 41b. The fourth calculator 150g calculates the mean value m of the voxel values ​​of the first map of the region 41b, which is a spherical VOI region, the entropy e of the voxel values ​​of the first map, and the standard deviation σ of the voxel values ​​of the first map, and then calculates the value m / e / σ by dividing the mean value m of the voxel values ​​of the first map of the voxels located in the region 41b, which is a spherical VOI region, by the entropy e of the voxel values ​​of the first map of the voxels located in the region 41b, which is a spherical VOI region, and dividing the mean value m by the standard deviation σ of the voxel values ​​of the first map of the voxels located in the region 41b, which is a spherical VOI region. The fourth calculator 150g acquires the calculated value of m / e / σ as the value of the fourth map for the central voxel 40b.

[0132] The fourth calculation unit 150g performs the above-described process by shifting the central voxel by one voxel at a time, as shown in FIG. 21 , to calculate the fourth map for the entire three-dimensional region, i.e., the entire first map. Here, the central voxel 40b indicates a certain central voxel, and the region 41b indicates the shape of a predetermined range within which the central voxel is integrated. The region 47 is the region from which data on the texture feature values ​​is obtained, and the boundary 45b is its boundary. Furthermore, the curve 46b indicates the trajectory of the center of a sphere, i.e., the central voxel 40b, when a sphere of the size indicated by the region 41b is inscribed in the curve 25b. In this case, the shape of the fourth map calculated by the fourth calculation unit 150g is the shape of the trajectory of the central voxel when a spherical VOI is inscribed in the first map, i.e., the shape of the curve 46b. The fourth calculation unit 150g may display the fourth map calculated in step S700 on the display 135 as a display unit, or may further apply a Gaussian filter to the fourth map calculated in this step.

[0133] 22, the fourth calculator 150g may set a boundary 45b of a region 47 for which data on the value of the texture feature has been obtained as the range through which the central voxel 40b passes when calculating the fourth map. This makes the boundary of the region for which the fourth map is obtained equal to the boundary 45b of the region 47 for which data on the value of the texture feature has been obtained.

[0134] 22 , when integrating the texture feature for each central voxel 40b, the range of integration may fall outside the region 47 for which texture feature value data has been obtained, resulting in missing texture feature values. In such cases, the fourth calculation unit 150g removes the range in which the texture feature value is missing from the integration range, calculates the average value of the texture feature values ​​for the range excluding the range in which the texture feature value is missing, and generates the fourth map. That is, the fourth calculation unit 150g calculates the fourth map for the same region 47 as the first map by calculating the central voxel value based on the texture feature in a region excluding the region in which the texture feature has not been calculated from the predetermined range in which the texture feature is integrated.

[0135] 20 to 22 , the first map is a map indicating the magnitude of a quantity corresponding to volume magnetic susceptibility calculated based on data collected from an imaging target by executing a pulse sequence, and the application to a magnetic resonance imaging apparatus has been described. However, the embodiment is not limited to this. The first map may be another type of map, and the embodiment is also applicable to image processing apparatuses that handle medical images other than those of magnetic resonance imaging apparatuses. That is, the embodiment is also applicable to an image processing apparatus that includes a calculation unit that calculates a central voxel value based on texture feature values ​​of a predetermined spherical range centered on the central voxel of the voxel values ​​of the first map, and moves the predetermined range and the central voxel by one voxel to generate a map for each voxel.

[0136] Subsequently, in step S700, the display control unit 150c displays the fourth map for each voxel generated in step S600 on the display 135 as a display unit. An example of the fourth map displayed in this manner is shown in Fig. 23. An image 91 is the fourth map generated for each voxel.

[0137] However, the embodiment is not limited to this, and the display control unit 150c may display the third map calculated for each functional unit in step S400 and the fourth map calculated for each voxel in step S600 side by side on the display 135 as a display unit.

[0138] For example, an example of such a display is shown in Figure 24. The upper part of Figure 24 shows an example of a case in which a third map generated by performing parcel analysis using parcels of the cortical system and a fourth map generated by performing voxel-based data processing are displayed side by side. That is, the display control unit 150c displays the third map calculated for each functional unit of the cortical system parcels in step S400 as image 50 and the fourth map calculated for each voxel in step S600 as image 51 side by side on the display 135 as a display unit.

[0139] 24 shows an example of a case in which a third map generated by parcel analysis using a standard parcel atlas of the white matter system and a fourth map generated by voxel-based data processing are displayed side by side. That is, the display controller 150c displays the third map calculated for each functional unit of the parcels of the white matter system in step S400 as image 52 and the fourth map calculated for each voxel in step S600 as image 53 side by side on the display 135.

[0140] The display control unit 150c may accept an input from the user and switch between different standard parcel atlases for display, or may regularly switch between different standard parcel atlases for display.

[0141] The method of the embodiment can also be used, for example, to observe changes over time and use the results for diagnosis. For example, by tracking changes over time in texture features using the method of the embodiment, it is possible to observe changes over time in the accumulation of Aβ protein in a patient, which can be useful for diagnosis.

[0142] In the above embodiment, the formula (m / e / σ) is used as the texture feature. However, the present invention is not limited to this. For example, the formula (m / e / σ) may be used, which is a first-order wide-area texture feature based on a histogram, for example, to express the texture feature by using a mean value m and a variance σ. 2 The texture feature may be any one of m / e / σ, standard deviation σ, entropy e, Z value, and skewness. Similarly, the texture feature may be one or more of fractal dimension (value), GLCM (Gray-Level Co-occurrence Matrix), GLSZM (Gray Level Size Zone Matrix), NGTDM (Neighborhood Gray-Tone-Difference Matrix), LBP (Local Binary Patterns), Laplacian distribution, Run Length matrix, and frequency filtering (such as wavelet transform or Fourier transform). Furthermore, the texture feature may be any one of m / e / σ, mean value m, variance σ, and so on. 2 The texture feature may be an expression (e.g., m / e or m×e / σ) using two or more of the following: σ, standard deviation σ, entropy e, Z value, skewness, and fractal dimension value, and an operator. Furthermore, second-order local or higher-order texture features that are not based on histograms may also be used. In other words, any texture feature that can characteristically capture and display brain lesions may be used.

[0143] According to at least one of the embodiments described above, it is possible to provide images with sufficient diagnostic capability for brain lesions such as Alzheimer's disease and demyelination.

[0144] Although several embodiments have been described, these embodiments are presented as examples and are not intended to limit the scope of the invention. These embodiments can be implemented in various other forms and can be variously replaced.

[0145] REFERENCE SIGNS LIST 100 Magnetic resonance imaging apparatus 108 Transmission circuit 110 Reception circuit 120 Sequence control unit 130 Image processing device 132 Memory 134 Input device 135 Display 150 Processing circuit 150a Control unit 150b Generation unit 150c Display control unit 150d First calculation unit 150e Second calculation unit 150f Third calculation unit 150g Fourth calculation unit 150h Range selection unit

Claims

1. A magnetic resonance imaging apparatus comprising: a sequence control unit that executes a pulse sequence to collect data from an imaging subject; a first calculation unit that calculates a first map indicating the magnitude of a quantity corresponding to volume magnetic susceptibility based on the data; a range selection unit that determines a selection range of phase slope values ​​based on a phase change value and an echo time (TE) value of the first map; a second calculation unit that calculates a second map indicating a feature amount of voxel values ​​of the first map or a portion of an area within the first map based on the first map and the selection range; and a display control unit that displays an image on a display unit based on the second map.

2. A magnetic resonance imaging apparatus according to claim 1, wherein the feature amount is calculated based on a texture feature amount of the first map or a partial area within the first map.

3. A magnetic resonance imaging apparatus as described in claim 2, wherein the feature is calculated based on texture features of all voxel values ​​contained in the first map or a portion of the first map, or voxel values ​​within a certain range of values.

4. The magnetic resonance imaging apparatus of claim 2, wherein the feature value is the average value of the voxel values ​​of the first map divided by the entropy of the voxel values ​​of the first map divided by the standard deviation of the voxel values ​​of the first map.

5. The magnetic resonance imaging apparatus of claim 1, further comprising a fourth calculation unit that calculates a central voxel value based on texture features of a predetermined range of voxel values ​​of the first map centered on a central voxel, and moves the central voxel by one voxel to generate a fourth map for each voxel, and the display control unit displays the fourth map for each voxel on the display unit.

6. The magnetic resonance imaging apparatus of claim 5, wherein the fourth calculation unit calculates the central voxel value by dividing the average value of a predetermined range centered on the central voxel by the entropy of the voxel values ​​in the predetermined range and dividing the result by the standard deviation of the voxel values ​​in the predetermined range.

7. A magnetic resonance imaging apparatus according to claim 5, wherein the predetermined range is a range having a rectangular parallelepiped or cubic shape.

8. A magnetic resonance imaging apparatus according to claim 5, wherein the predetermined range is a range having a three-dimensional spherical shape.

9. A magnetic resonance imaging apparatus as described in claim 5, wherein the fourth calculation unit calculates the fourth map for the same area as the first map by calculating the central voxel value based on the texture feature in an area of ​​the specified range excluding areas in which the texture feature has not been calculated.

10. The magnetic resonance imaging apparatus according to claim 1, wherein the pulse sequence is a GRE (Gradient-Echo) pulse sequence executed for a plurality of TEs (Echo Times).

11. The magnetic resonance imaging apparatus of claim 1, wherein the pulse sequence is a sequence that generates a plurality of echoes, and the first calculation unit calculates the first map based on a value obtained by dividing a value of a phase change between the plurality of echoes by a value of a change in TE between the echoes.

12. The magnetic resonance imaging apparatus of claim 1, wherein the pulse sequence is a sequence that generates a plurality of echoes, and the first calculation unit calculates, as the first map, a value obtained by dividing a value of a phase change between each echo by a value of a change in TE between each echo and linearly transforming the value so that the converted value is positive when the imaging subject contains a paramagnetic component.

13. The magnetic resonance imaging apparatus according to claim 1, wherein the pulse sequence is a sequence that generates a plurality of echoes, and the first calculation unit calculates the first map after applying a Wiener filter to a phase image for each of the echoes.

14. The magnetic resonance imaging apparatus according to claim 1, further comprising a third calculation unit that calculates a third map by aggregating the second maps for each functional unit of the imaging subject.

15. The magnetic resonance imaging apparatus according to claim 14, wherein said display control unit refers to a color table and causes said display unit to display said third map in color for each of said functional units.

16. The magnetic resonance imaging apparatus according to claim 14, wherein the display control unit causes the third map to be superimposed on a T1 weighted image or a blood vessel image on the display unit.

17. The magnetic resonance imaging apparatus of claim 1, further comprising a fourth calculation unit that applies an image filter to the second map to generate a fourth map for each voxel, and the display control unit causes the display unit to display the fourth map for each voxel.

18. The magnetic resonance imaging apparatus according to claim 17, wherein the display control unit causes the display unit to display in parallel the third map calculated for each functional unit and the fourth map calculated for each voxel.

19. An image processing device comprising: a first calculation unit that calculates a first map indicating the magnitude of a quantity corresponding to volume magnetic susceptibility based on data collected from an imaging subject by executing a pulse sequence; a range selection unit that determines a selection range of phase slope values ​​based on a phase change value and an echo time (TE) value of the first map; a second calculation unit that calculates a second map indicating a feature amount of voxel values ​​of the first map or a portion of an area within the first map based on the first map and the selection range; and a display control unit that analyzes the second map and displays an image on a display unit.

20. An image processing device having a calculation unit that calculates a central voxel value based on the texture features of a predetermined spherical range centered on a central voxel of the voxel values ​​of a first map, and moves the predetermined range and the central voxel by one voxel at a time to generate a map for each voxel.

Images