Magnetic resonance numerical simulation device, method, and program

Abstract

To improve the calculation efficiency of magnetic resonance numerical value simulation.SOLUTION: A magnetic resonance numerical value simulation device according to an embodiment comprises an input unit and a numerical value calculation unit. The input unit accepts the input of numerical values before update of magnetization and numerical values before update of partial differentiation in each of the spatial direction and / or angular frequency direction with respect to the magnetization for an isochromat. The numerical value calculation unit calculates the update formula that represents the time evolution of magnetic resonance by using some of the numerical values before the update of magnetization and the numerical values before the update of partial differentiation, and calculates the numerical values after the update of magnetization and the numerical values after the update of spatial partial differentiation in the calculation target direction in the spatial direction and / or angular frequency direction with respect to the magnetization for the isochromat.SELECTED DRAWING: Figure 1

Description

[Technical field]

[0001] The embodiments disclosed in this specification and the drawings relate to a magnetic resonance numerical simulation apparatus, method, and program. [Background technology]

[0002] In magnetic resonance numerical simulation, the Bloch equations are used to simulate magnetic resonance numerical imaging. The Bloch equations are equations that describe the time evolution of magnetization in a magnetic field, and are the fundamental equations of magnetic resonance imaging. MAGSI (Magnetic Gradient Simulation of Intra-voxel dephasing) is a numerical solution method for the Bloch equations. MAGSI sets a virtual molecule at the center of the voxel, and calculates the magnetization at any position within the voxel using the spatial partial derivative of the magnetization under the assumption that the magnetic field within the voxel is linear. However, since MAGSI naively calculates the spatial partial derivative of the magnetization, there is room for improvement in the calculation efficiency. [Prior art documents] [Patent documents]

[0003] [Patent Document 1] JP 2023-003389 A [Non-patent literature]

[0004] [Non-Patent Document 1] Hidenori Takeshima. “A Fast and Practical Computation Method for Magnetic Resonance Simulators” Magnetic Resonance in Medicine, DOI: 10.1002 / mrm.29646 [Non-Patent Document 2] Thies H. Jochimsen, Andreas Schafer, Roland Bammer, Michael E. Moseley, “Efficient simulation of magnetic resonance imaging with Bloch-Torrey equations using intra-voxel magnetization gradients” Journal of Magnetic Resonance, 2006; 180:29-38. [Non-Patent Document 3] Graf C, Rund A, Aigner CS, Stollberger R. “Accuracy and Performance Analysis for Bloch and Bloch-McConnell Simulation Methods” Journal of Magnetic Resonance. 2021; 329: 107011. [Non-Patent Document 4] Latta P, Gruwel MLH, Jellus V, Tomanek B. “Bloch simulations with intra-voxel spin dephasing” Journal of Magnetic Resonance. 2010; 203:44-51. Summary of the Invention [Problem to be solved by the invention]

[0005] One of the problems to be solved by the embodiments disclosed in this specification and the drawings is to improve the calculation efficiency of magnetic resonance numerical simulation. However, the problems to be solved by the embodiments disclosed in this specification and the drawings are not limited to the above problem. Problems corresponding to the effects of each configuration shown in the embodiments described later can also be positioned as other problems. [Means for solving the problem]

[0006] A magnetic resonance numerical simulation apparatus according to an embodiment includes an input unit and a numerical calculation unit. The input unit inputs, for an isochromat, a pre-update value of magnetization and a pre-update value of a partial differential in each of the spatial direction and / or the angular frequency direction with respect to the magnetization. The numerical calculation unit uses some of the pre-update value of the magnetization and the pre-update value of the partial differential to calculate an update equation that represents the time evolution of magnetic resonance, and calculates, for the isochromat, a post-update value of the magnetization and a post-update value of the partial differential in a calculation target direction among the spatial direction and / or the angular frequency direction with respect to the magnetization. [Brief description of the drawings]

[0007] [Figure 1] FIG. 1 is a diagram showing an example of the configuration of a magnetic resonance numerical simulation apparatus according to the first embodiment. [Diagram 2] FIG. 2 is a diagram illustrating components constituting a magnetization vector M(t) that is the object of calculation in a magnetic resonance numerical simulation. [Diagram 3] FIG. 3 is a diagram showing a schematic diagram of a calculation process of the magnetization vector M(t+Δt) according to a comparative example. [Figure 4] FIG. 4 is a diagram showing calculation targets and contributing terms related to the magnetization vector M(t). [Diagram 5] FIG. 5 is a diagram showing a schematic diagram of the calculation process of the update equation for the magnetization vector. [Figure 6] FIG. 6 is a diagram illustrating a processing procedure of a magnetic resonance numerical simulation according to Example 1 of the first embodiment. [Figure 7] FIG. 7 is a diagram showing an example of a pulse sequence according to the spin echo method. [Figure 8] FIG. 8 is a diagram illustrating a processing procedure of a magnetic resonance numerical simulation according to Example 2 of the first embodiment. [Figure 9] FIG. 9 is a diagram illustrating an example of the configuration of a magnetic resonance numerical simulation apparatus according to the second embodiment. [Figure 10]FIG. 10 is a diagram illustrating a processing procedure of a magnetic resonance numerical simulation according to Example 1 of the second embodiment. [Figure 11] FIG. 11 is a diagram illustrating a processing procedure of a magnetic resonance numerical simulation according to Example 2 of the second embodiment. DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

[0008] Hereinafter, a magnetic resonance numerical simulation apparatus, method, and program according to the present embodiment will be described in detail with reference to the drawings.

[0009] The magnetic resonance numerical simulation device according to the present embodiment is a computer that executes a magnetic resonance numerical simulation. The magnetic resonance numerical simulation according to the present embodiment simulates the time evolution of magnetic resonance by a given pulse sequence for an isochromat in each voxel set in an inspection object model. Specifically, MAGSI is used to calculate magnetization at any position in a voxel using spatial partial derivatives of magnetization under the assumption that the magnetic field in the voxel is linear. The magnetic resonance numerical simulation device, method, and program according to the present embodiment aim to improve the calculation efficiency of magnetization and / or partial derivatives to improve calculation efficiency.

[0010] (First embodiment) Fig. 1 is a diagram showing an example of the configuration of a magnetic resonance numerical simulation device 1 according to the first embodiment. As shown in Fig. 1, the magnetic resonance numerical simulation device 1 is a computer having a processing circuit 11, a communication device 12, a display device 13, an input device 14, and a storage device 15. Data communication between the processing circuit 11, the communication device 12, the display device 13, the input device 14, and the storage device 15 is performed via a bus.

[0011] The processing circuit 11 has a processor such as a CPU (Central Processing Unit). The processor starts a magnetic resonance numerical simulation program installed in the storage device 15 or the like, thereby realizing an input function 111, a numerical calculation function 112, an output control function 113, and the like. Each of the functions 111 to 113 is not limited to being realized by a single processing circuit. A processing circuit may be configured by combining a plurality of independent processors, and each processor may execute a magnetic resonance numerical simulation program to realize each of the functions 111 to 113. In addition, each of the functions 111 to 113 may be implemented as a module that constitutes the magnetic resonance numerical simulation program, or may be implemented as individual hardware.

[0012] By the input function 111, the processing circuit 11 inputs various information used in the magnetic resonance numerical simulation. As an example, the processing circuit 11 inputs a sequence fragment of a pulse sequence to be calculated. The numerical calculation function 112 performs a magnetic resonance numerical calculation for each sequence fragment. By the input function 111, the processing circuit 11 inputs, for an isochromat, a numerical value before the update of magnetization and a numerical value before the update of partial differentials in each of the spatial direction and the angular frequency direction for the magnetization. An isochromat means a set of nuclear magnetic moments having the same frequency set in a voxel. As the numerical value before the update, a numerical value manually input by the input device 14, a numerical value after the update calculated by the numerical calculation function 112, or an arbitrary value set in advance may be input.

[0013] By using the numerical calculation function 112, the processing circuit 11 executes an update calculation of the magnetization vector for each of the sequence fragments input by the input function 111. As an example, the processing circuit 11 calculates an update formula (hereinafter simply referred to as an update formula) expressing the time evolution of magnetic resonance using some of the pre-update values ​​of the magnetization and the pre-update values ​​of the partial differentials input by the input function 111, and calculates, for the isochromat, an updated value of the magnetization and an updated value of the spatial partial differential in the calculation target direction among the spatial directions and / or angular frequency directions for the magnetization. Specifically, the processing circuit 11 calculates the update formula using the pre-update values ​​of the magnetization without using the pre-update values ​​of the partial differentials, and calculates the updated value of the magnetization. In addition, the processing circuit 11 calculates the update formula using the pre-update values ​​of the partial differentials in the calculation target direction without using the pre-update values ​​of the partial differentials in the non-calculation target directions among the spatial directions and / or angular frequency directions, and calculates the updated value of the partial differential in the calculation target direction. Hereinafter, the numerical value of magnetization in the magnetization vector is called magnetization value, and the numerical value of the partial differential of magnetization is called partial differential value. Bloch equations or their developed forms are used as update equations. The Bloch equations according to this embodiment include not only the original Bloch equations, but also developed forms that add a diffusion term by Torrey or a magnetization transfer term by McConnell.

[0014] The processing circuit 11 may generate data (hereinafter, simulation data) corresponding to the purpose of the magnetic resonance numerical simulation based on the updated magnetization value and the updated partial differential value in the calculation target direction by using the numerical calculation function 112. Examples of the simulation data include a predicted output value of an A / D converter and an MR image based on the predicted output value. The magnetization value and the spatial partial differential value may be set in the simulation data.

[0015] By using the output control function 113, the processing circuit 11 outputs simulation data such as updated magnetization values ​​and partial differential values ​​calculated by the numerical calculation function 112. As an example, the processing circuit 11 displays the simulation data on the display device 13, stores it in the storage device 15, or transmits it to another computer via the communication device 12.

[0016] The communication device 12 is an interface that connects to a magnetic resonance imaging apparatus, a workstation, a PACS (Picture Archiving and Communication System), a HIS (Hospital Information System), a RIS (Radiology Information System), etc. via a LAN (Local Area Network) etc. The communication device 12 transmits and receives various data to and from the connected devices.

[0017] The display device 13 displays various data in accordance with the output control function 113 of the processing circuit 11. As the display device 13, a liquid crystal display (LCD), a cathode ray tube (CRT) display, an organic electroluminescence display (OELD), a plasma display, or any other display may be appropriately used. The display device 13 may also be a projector.

[0018] The input device 14 accepts various input operations from a user, converts the accepted input operations into electrical signals, and outputs the electrical signals to the processing circuit 11. Specifically, the input device 14 is connected to input devices such as a mouse, a keyboard, a trackball, a switch, a button, a joystick, a touchpad, and a touch panel display. The input device 14 outputs electrical signals corresponding to the input operations to the input device 14 to the processing circuit 11. The input device 14 may also be an input device provided in another computer connected via a network or the like. The input device 14 may also be a voice recognition device that converts a voice signal collected by a microphone into an instruction signal.

[0019] The storage device 15 is a storage device such as a ROM (Read Only Memory), a RAM (Random Access Memory), a HDD (Hard Disk Drive), an SSD (Solid State Drive), or an integrated circuit storage device that stores various data. The storage device 15 stores, for example, a magnetic resonance numerical simulation program. In addition to the above-mentioned storage devices, the storage device 15 may be a portable storage medium such as a CD (Compact Disc), a DVD (Digital Versatile Disc), or a flash memory, or a drive device that reads and writes various information between a semiconductor memory element, etc. The storage device 15 may be in another computer connected to the magnetic resonance numerical simulation device 1 via a network.

[0020] As described above, the magnetic resonance numerical simulation device according to this embodiment includes a processing circuit 11. The processing circuit 11 inputs, for an isochromat, a pre-update value of magnetization and a pre-update value of a partial differential in each of the spatial direction and / or the angular frequency direction with respect to the magnetization. The processing circuit 11 uses some of the pre-update values ​​of the magnetization and the partial differential to calculate an update equation that represents the time evolution of magnetic resonance, and calculates, for the isochromat, a post-update value of magnetization and a post-update value of a partial differential in a calculation target direction of the spatial direction and / or the angular frequency direction.

[0021] Specifically, when calculating the magnetization value before the update, the processing circuit 11 uses the magnetization value before the update as a part of the numerical value. In other words, the processing circuit 11 may not use the partial differential value before the update as a part of the numerical value. When calculating the partial differential value in the calculation target direction before the update, the processing circuit 11 uses the magnetization value before the update and the partial differential value in the calculation target direction as a part of the numerical value. In other words, the processing circuit 11 may not use the partial differential value in the non-calculation target direction before the update as a part of the numerical value. Note that each of the partial differential value in the spatial direction and / or the partial differential value in the angular frequency direction may have one or more components. For example, the partial differential value in the spatial direction may have multiple components such as a partial differential value in the readout gradient magnetic field direction, a partial differential value in the phase encoding gradient magnetic field direction, and a partial differential value in the slice selection gradient magnetic field direction. The partial differential value in the calculation target direction refers to a part of these multiple components. Similarly, the partial differential value in the angular frequency direction may have one or more components. In other words, when calculating the partial differential value in the direction to be calculated before the update, the processing circuit 11 does not need to use all the numerical values ​​of the magnetization value before the update, the partial differential value in the direction to be calculated, and the partial differential value in the direction not to be calculated in order to calculate the updated numerical value of the partial differential in the direction to be calculated among the spatial directions and / or angular frequency directions, but only uses the magnetization value before the update and the partial differential value in the direction to be calculated.

[0022] According to the above configuration, this embodiment can reduce the amount of calculation of the update equation compared to the case where all the numerical values ​​of the updated magnetization value and partial differential value are used. As a result, the calculation time and memory consumption can be reduced, and the calculation efficiency of the magnetic resonance numerical simulation can be improved.

[0023] Here, this embodiment will be compared with Non-Patent Document 4. Non-Patent Document 4 states that for each of the x, y, and z directions in the Bloch equations, -6The magnetization when a small deviation of about is given is calculated, and the magnetization difference with respect to each axis is regarded as a differential to calculate the phase, thereby simulating intra-voxel spin dephasing. Compared to Non-Patent Document 4, this embodiment does not replace the differential with a differential, but calculates a spatial partial differential value. In addition, this embodiment utilizes the fact that there is no interaction between spatial partial differentials in different partial differential directions, and applies the update formula only to the components that contribute to the calculation of the updated numerical value. This makes it possible for this embodiment to calculate the updated numerical value in a short time without neglecting the accuracy of the result.

[0024] An example of the magnetic resonance numerical simulation apparatus 1 according to the first embodiment will be described below. Note that the magnetization vectors and update formulas shown below are merely examples, and the present embodiment is not limited to these.

[0025] FIG. 2 is a diagram illustrating components constituting a magnetization vector M(t) that is a calculation target of the magnetic resonance numerical simulation according to the present embodiment. The magnetization vector M(t) shown in FIG. 2 represents a vector sum of magnetization moments of isochromats in a voxel at time t. Specifically, the magnetization vector M(t) is M x , M y , M z , ∂ x M x , ∂ x M y , ∂ x M z , ∂ y M x , ∂ y M y , ∂ y M z , ∂ z M x , ∂ z M y , ∂ z M z , ∂ ω M x , ∂ ω M y , ∂ ω M z and M0. x , M y and Mz represents the spatial component of magnetization. Specifically, M x , M y and M z represent the magnetization components in the x, y, and z directions, respectively. x M x , ∂ x M y , ∂ x M z , ∂ y M x , ∂ y M y , ∂ y M z , ∂ z M x , ∂ z M y and ∂ z M z represents the spatial partial differential of the magnetization in each of multiple spatial directions. Specifically, ∂ x M x , ∂ x M y and ∂ x M z are the magnetization components M x , M y and M z represents the partial differential component in the x direction with respect to ∂ y M x , ∂ y M y , ∂ y M z are the magnetization components M x , M y and M z represents the partial differential component in the y direction with respect to ∂ z M x , ∂ z M y and ∂ z M z are the magnetization components M x , M y and M z represents the partial differential component in the z direction with respect to ∂ ω M x , ∂ ω M y , ∂ ω M z are the magnetization components M x , My and M z represents the partial differential component in the angular frequency direction with respect to . M0 represents the magnetization of the isochromat in thermal equilibrium. T means the transposition.

[0026] In this embodiment, the x-axis, y-axis, and z-axis can be set in any direction, but for convenience of calculation, the readout gradient magnetic field direction is set to the x-axis, the phase encoding gradient magnetic field direction is set to the y-axis, and the slice selection gradient magnetic field direction is set to the z-axis. Also, partial differential components in the angular frequency direction may or may not be taken into consideration.

[0027] The processing circuit 11 applies an update equation to the magnetization vector M(t) before the update for each of the multiple voxels set in the inspection target model, and calculates the magnetization vector M(t+Δt) after a given time step Δt as the updated magnetization vector. At this time, the processing circuit 11 calculates the magnetization at any position in the voxel using the spatial partial derivative of the magnetization under the assumption that the magnetic field in each voxel is linear. The updated magnetization at any position in the voxel can be calculated using the time derivative and the spatial partial derivative.

[0028] 3 is a diagram showing a schematic diagram of a calculation process of a magnetization vector M(t+Δt) according to a comparative example. As shown in FIG. 3, in the comparative example, one update formula F for the magnetization vector M(t) is used. ALL Calculate the update formula F ALL is each component M of the magnetization vector M(t). x , M y , M z , ∂ x M x , ∂ x M y , ∂ x M z , ∂ y M x , ∂ y M y , ∂ y M z , ∂ z M x , ∂ z M y , ∂ z Mz , ∂ ω M x , ∂ ω M y , ∂ ω M z In this case, all of these components M x , M y , M z , ∂ x M x , ∂ x M y , ∂ x M z , ∂ y M x , ∂ y M y , ∂ y M z , ∂ z M x , ∂ z M y , ∂ z M z , ∂ ω M x , ∂ ω M y , ∂ ω M z Using and M0, the update formula F ALL Calculate the updated magnetization value M x,y,z (t+Δt) and the updated spatial partial differential value ∂ x M x,y,z (t+Δt), ∂ y M x,y,z (t+Δt) and ∂ z M x,y,z (t+Δt) is calculated. Note that M x,y,z (t+Δt) is M x (t+Δt), M y (t+Δt), M z This is an abbreviated notation of (t+Δt).

[0029] FIG. 4 is a diagram showing a calculation target and a contribution term related to a magnetization vector M(t). A calculation target means a component of interest among the components constituting a magnetization vector. In terms of spatial partial differential values, a calculation target is determined by a spatial direction in which partial differentiation is performed (hereinafter, referred to as a partial differential direction), such as the x direction, y direction, and z direction. A contribution term is a component used to calculate a calculation target, in other words, an update formula F ALL The term "non-contributing term" refers to a component that is used in the term used to calculate the calculation target among the terms in the update formula F. ALL In the case of a spatial partial differential value, a non-contributing term means a component that is not included in a term used in the calculation of the target value. In the case of a spatial partial differential value, a non-contributing term means a partial differential direction other than the target value of the calculation.

[0030] As shown in FIG. 4(1), the calculation target is the magnetization value M x , M y and M z If so, the contribution term is M x , M y , M z , M0, and the non-contributing terms are the partial differentials in all spatial and angular frequency directions, ∂ x M x , ∂ x M y , ∂ x M z , ∂ y M x , ∂ y M y , ∂ y M z , ∂ z M x , ∂ z M y , ∂ z M z , ∂ ω M x , ∂ ω M y , ∂ ω M z As shown in Figure 4 (2), the calculation target is the spatial partial differential value ∂ x M x , ∂ x M y , ∂ xM z If so, the contribution term is M x , M y , M z , ∂ x M x , ∂ x M y , ∂ x M z and the non-contributing terms are the partial differentials ∂ in the y, z, and ω directions. y M x , ∂ y M y , ∂ y M z , ∂ z M x , ∂ z M y , ∂ z M z , ∂ ω M x , ∂ ω M y , ∂ ω M z As shown in FIG. 4(3), the calculation target is the spatial partial differential value ∂ y M x , ∂ y M y , ∂ y M z If so, the contribution term is M x , M y , M z , ∂ y M x , ∂ y M y , ∂ y M z and the non-contributing terms are the partial differentials ∂ in the x, z, and ω directions. x M x , ∂ x M y , ∂ x M z , ∂ z M x , ∂ z M y , ∂ z M z , ∂ ω M x , ∂ ω M y , ∂ ω M zAs shown in FIG. 4(4), the calculation target is the spatial partial differential value ∂ z M x , ∂ z M y and ∂ z M z If so, the contribution term is M x , M y , M z , ∂ z M x , ∂ z M y , ∂ z M z and the non-contributing terms are the partial differentials ∂ in the x, y, and ω directions. y M x , ∂ y M y , ∂ y M z , ∂ z M x , ∂ z M y , ∂ z M z , ∂ ω M x , ∂ ω M y , ∂ ω M z As shown in FIG. 4(5), the calculation target is the partial differential value ∂ ω M x , ∂ ω M y and ∂ ω M z If so, the contribution term is M x , M y , M z , ∂ ω M x , ∂ ω M y , ∂ ω M z and the non-contributing terms are the spatial partial derivatives in all spatial directions, ∂ x M x , ∂ x M y , ∂ x M z , ∂ y M x , ∂ y M y , ∂y M z , ∂ z M x , ∂ z M y , ∂ z M z and M0.

[0031] For example, when the Bloch equation includes the Torrey diffusion term, the spatial partial differential value ∂ x M x , ∂ x M y , ∂ x M z , ∂ y M x , ∂ y M y , ∂ y M z , ∂ z M x , ∂ z M y , ∂ z M z is added for all calculations (even though it is a non-contributing term in the Bloch equations).

[0032] As shown in (1) to (5) of FIG. 4, it can be seen that the magnetization component of the magnetization vector, the partial differential component in the x direction, the partial differential component in the y direction, the partial differential component in the z direction, and the partial differential component in the ω direction each have contributing terms and non-contributing terms. For example, it is not necessary to use partial differential values ​​in all directions to calculate the magnetization value. In another example, it is not necessary to use spatial partial differential values ​​in other spatial directions to calculate the spatial partial differential component value in a certain direction. This can also be expressed as there being no interaction between partial differential components in different spatial directions. Therefore, the magnetization value, the partial differential value in the x direction, the partial differential value in the y direction, the partial differential value in the z direction, and the partial differential value in the ω direction are calculated separately.

[0033] 5 is a diagram showing a schematic diagram of a calculation process of an update formula of a magnetization vector according to the first embodiment. As shown in FIG. 5, in the first embodiment, the update formula of a magnetization vector is an update formula of a magnetization component (hereinafter, magnetization update formula) F Mand the spatial partial differential component update formula (hereinafter referred to as the spatial partial differential update formula) F ∂p and the update formula for the angular frequency partial differential component (hereinafter referred to as the angular frequency partial differential update formula or the ω direction partial differential update formula) F ∂ω The spatial partial differential update formula F ∂p is further expressed as the update equation for the spatial partial differential component in the x direction (hereinafter, the x direction partial differential update equation) F ∂x , the update formula for the spatial partial differential component in the y direction (hereinafter, the y direction partial differential update formula) F ∂y , the update equation for the spatial partial differential component in the z direction (hereinafter, the z direction partial differential update equation) F ∂z It can be divided into:

[0034] As shown in FIG. 5, the processing circuit 11 calculates the magnetization value M x、y、z When calculating (t+Δt), the contribution term M x , M y , M z , M0, and the non-contributing term ∂ x M x , ∂ x M y , ∂ x M z , ∂ y M x , ∂ y M y , ∂ y M z , ∂ z M x , ∂ z M y , ∂ z M z , ∂ ω M x , ∂ ω M y , ∂ ω M z Without using the magnetization update formula F M The processing circuit 11 calculates the x-direction partial differential value ∂ x M x、y、z When calculating (t+Δt), the contribution term M x , M y , M z , ∂ x M x , ∂ x M y , ∂ x M z , M0, and the non-contributing term ∂y M x , ∂ y M y , ∂ y M z , ∂ z M x , ∂ z M y , ∂ z M z , ∂ ω M x , ∂ ω M y , ∂ ω M z Instead of using the x-direction partial differential update formula F ∂x The processing circuit 11 calculates the y-direction partial differential value ∂ y M x、y、z When calculating (t+Δt), the contribution term M x , M y , M z , ∂ y M x , ∂ y M y , ∂ y M z , M0, and the non-contributing term ∂ x M x , ∂ x M y , ∂ x M z , ∂ z M x , ∂ z M y , ∂ z M z , ∂ ω M x , ∂ ω M y , ∂ ω M z Instead of using ∂y The processing circuit 11 calculates the z-direction partial differential value ∂ z M x、y、z When calculating (t+Δt), the contribution term M x , M y , M z , ∂ z M x , ∂ z M y , ∂ z M z , M0, and the non-contributing term ∂x M x , ∂ x M y , ∂ x M z , ∂ y M x , ∂ y M y , ∂ y M z , ∂ ω M x , ∂ ω M y , ∂ ω M z Instead of using the z-direction partial differential update formula F ∂z The processing circuit 11 calculates the ω-direction partial differential value ∂ ω M x、y、z When calculating (t+Δt), the contribution term M x , M y , M z , ∂ ω M x , ∂ ω M y , ∂ ω M z , M0, and the non-contributing term ∂ x M x , ∂ x M y , ∂ x M z , ∂ y M x , ∂ y M y , ∂ y M z , ∂ z M x , ∂ z M y , ∂ z M z Instead of using the ω-direction partial differential update formula F ∂ω Calculate the following.

[0035] As described above, the update formula F according to the first embodiment ALL are the mutually separated magnetization renewal equations F M and the spatial partial differential update equation F ∂ and the angular frequency partial differential update formula F ∂ω In order to calculate the updated magnetization value, the processing circuit 11 refers only to the magnetization value before the update and calculates the magnetization update formula F MHere, "referring to" means accessing and reading out the pre-update numerical value stored in the storage device 15. In order to calculate the updated spatial partial differential value in the calculation target direction, the processing circuit 11 refers only to the pre-update magnetization value and the spatial partial differential value in the calculation target direction, and calculates the spatial partial differential update formula F ∂p The spatial partial differential update formula F ∂p are the mutually separated x-direction partial differential update equations F ∂x , y-direction partial differential update equation F ∂y and the z-direction partial differential update equation F ∂z In order to calculate the updated x-direction partial differential value, the processing circuit 11 refers only to the magnetization value and the x-direction partial differential value before the update, and calculates the x-direction partial differential update formula F ∂x In addition, in order to calculate the updated y-direction partial differential value, the processing circuit 11 calculates the y-direction partial differential update formula F ∂y In order to calculate the updated z-direction partial differential value, the processing circuit 11 calculates the z-direction partial differential update formula F by referring only to the magnetization value and the z-direction partial differential value before the update. ∂z In order to calculate the updated angular frequency direction partial differential value, the processing circuit 11 calculates the angular frequency partial differential update formula F by referring only to the magnetization value and the angular frequency direction partial differential value before the update. ∂ω Calculate the following.

[0036] In this way, by separating the update equations related to the components that do not interact with each other, it becomes unnecessary to refer to the non-contributing terms when calculating the update equation related to the calculation target component. By calculating the update equation related to the calculation target component without referring to the non-contributing terms and referring to the contributing terms, the update equation F of the magnetization vector as shown in FIG. ALL In comparison with the case where the update formula is calculated, the calculation time and calculation load can be reduced.

[0037] Next, the processing procedure of the magnetic resonance numerical simulation according to the first embodiment will be described in two parts, Example 1 and Example 2.

[0038] [Example 1] The Bloch equation representing the time evolution of magnetic resonance is expressed by the following equation (1): Note that in this embodiment, the diffusion term is not taken into consideration.

[0039]

number

[0040] M is the magnetization vector, B is the effective magnetic field vector, and γ is the gyromagnetic ratio. T1 is the longitudinal relaxation time, and T2 is the transverse relaxation time. The x, y, and z components of the Bloch equation above can be expanded to the following equations (2), (3), and (4), respectively. Equations (2), (3), and (4) correspond to the magnetization renewal equations in the x, y, and z directions, respectively.

[0041]

number

[0042] As shown in the above formulas (2), (3), and (4), the calculation formula for the magnetization value includes a term related to magnetization and a term related to M0, but does not include a term related to the spatial partial differential. In this way, since it is not necessary to use the spatial partial differential value when calculating the magnetization value, the magnetization update formula and the spatial partial differential update formula can be separated. By separating them, it is possible to reduce the time required to calculate the updated magnetization value compared to the case where they are not separated.

[0043] Next, we derive the equation for the spatial partial differential of the Bloch equation. We assume that the differentials in the x, y, z, and t directions are all continuous within a small time. However, we ignore the positions where t becomes discontinuous as a result of the magnetic field being a piecewise constant. We also consider the static magnetic field, including the chemical shift, to be constant for a small voxel corresponding to one isochromat. Let p∈{x, y, z}, and then M x (t), M y (t) and M z The spatial partial derivatives of (t) are expressed by the following equations (5), (6), and (7), respectively.

[0044]

number

[0045] As in the above formulas (5), (6), and (7), the calculation formula for the spatial partial differential value in the calculation target direction p includes a partial differential term in the direction p, but does not include a partial differential term in other directions q (≠p). For example, when p=x in formulas (5), (6), and (7), only a term related to the calculation target direction x exists in formulas (5), (6), and (7), and no term including the non-calculation target directions y and / or z exists. In other words, since it is not necessary to use the spatial partial differential value in the non-calculation target direction when calculating the spatial partial differential value in the calculation target direction, the spatial partial differential update formula can be separated for each calculation target direction. Therefore, it is possible to reduce the time required to calculate the updated spatial partial differential value compared to the case where the formulas are not separated.

[0046] FIG. 6 is a diagram illustrating a processing procedure of a magnetic resonance numerical simulation according to Example 1. It is assumed that a pulse sequence to be calculated has already been set at the start of FIG. 6. The pulse sequence used in this embodiment is not particularly limited, and can be applied to any type of pulse sequence. Hereinafter, as an example, the pulse sequence used in this embodiment is a pulse sequence based on the spin echo method. It is assumed that the angular frequency directional partial differential is not taken into consideration in the example shown in FIG. 6.

[0047] FIG. 7 is a diagram showing an example of a pulse sequence by the spin echo method. In FIG. 7, the "RF" timeline shows a sequence of RF pulses applied by an RF transmission coil driven in response to a drive signal from a transmission circuit. The "Gx" timeline shows a sequence of readout gradient magnetic field pulses applied by a Gx coil driven in response to a drive signal from a gradient magnetic field power supply. The "Gy" timeline shows a sequence of phase encoding gradient magnetic field pulses applied by a Gy coil driven in response to a drive signal from a gradient magnetic field power supply. The "Gz" timeline shows a sequence of slice selection gradient magnetic field pulses applied by a Gz coil driven in response to a drive signal from a gradient magnetic field power supply. As shown in FIG. 7, the pulse sequence can be divided into a non-application period 41, which is a period during which an RF pulse is not applied, and an application period 42, which is a period during which an RF pulse is applied.

[0048] As shown in FIG. 7, in the spin echo method, a first slice selection gradient magnetic field pulse is applied together with a 90° pulse, a second slice selection gradient magnetic field pulse is applied together with a 180° pulse, and a phase encoding pulse of a given magnetic field gradient is applied. An echo signal (not shown) is generated so that the signal value peaks after the echo time has elapsed since the application of the 90° pulse. The echo signal is received by a receiving coil, and the received analog echo signal is converted into digital k-space data by an A / D converter (ADC: Analog digital converter). The k-space data is an example of a predicted output value of the A / D converter. When the echo signal is received, a readout gradient magnetic field pulse is applied. The time from the application of the 90° pulse to the application of the readout gradient magnetic field pulse is TR (time of repetition). The above pulse sequence is repeated for a given number of phase encoding steps.

[0049] 6, the processing circuitry 11 inputs a sequence fragment by implementing the input function 111 (step S1). The sequence fragment is a part of a pulse sequence divided according to an arbitrary rule. The pulse sequence may be divided every predetermined time, or may be divided by pulse type such as "RF", "Gx", "Gy", "Gz", or may be divided by a combination of these.

[0050] When step S1 is performed, the processing circuit 11 determines parameters of the magnetization update equation and the spatial partial differential update equation (hereinafter, update equation parameters) by implementing the numerical calculation function 112 (step S2). In step S2, the processing circuit 11 determines the update equation parameters by using the sequence fragment input in step S1. The update equation parameters are determined for each of all isochromats arranged on the model to be inspected. Examples of the update equation parameters include the T1 value, T2 value, and B value included in the update equation.

[0051] When step S2 is performed, the processing circuit 11, by implementing the input function 111, outputs the magnetization value M(t) before updating, the spatial partial differential value ∂ p M(t) and the magnetization value M0 in the thermal equilibrium state are input (step S3). In the first update calculation, M(t) before the update and ∂ p An arbitrary initial value is input for M(t). In the second and subsequent update operations, M(t) and ∂ p For M(t), M(t+Δt) and ∂ p M(t+Δt) is input. Any value may be input for the magnetization value M0. M(t) before the update and ∂ p It is assumed that M(t) and the magnetization value M0 are stored in the storage device 15 before step S2. The processing circuit 11 receives from the storage device 15 the pre-update M(t) and ∂ p M(t) and the magnetization value M0 are read out and input to the self.

[0052] When step S3 is performed, the processing circuit 11, by implementing the numerical calculation function 112, refers only to the magnetization values ​​M(t) and M0 to calculate the magnetization update equations (2), (3), and (4) and calculates the updated magnetization value M(t+Δt) (step S4). At this time, the processing circuit 11 calculates the spatial partial differential value ∂ p There is no need to refer to M(t). This improves the calculation speed of the magnetization value M(t+Δt). The calculated updated magnetization value M(t+Δt) is stored in the storage device 15. The magnetization update formula can be calculated using the Runge-Kutta method or the like.

[0053] When step S4 is performed, the processing circuit 11 realizes the numerical calculation function 112 to calculate the magnetization value M(t) and ∂ p By referring only to M(t), the spatial partial differential update equations (5), (6), and (7) are calculated, and the updated spatial partial differential value ∂ p In step S5, the processing circuit 11 calculates the spatial partial differential update equations (5), (6), and (7) for each direction p while switching the calculation target direction p between the x direction, the y direction, and the z direction in turn, to calculate the spatial partial differential value ∂ p At this time, the processing circuit 11 calculates the spatial partial differential value ∂ q There is no need to refer to M(t). This allows each spatial partial differential value ∂ p The calculation speed of M(t+Δt) is improved. The calculated updated spatial partial differential value ∂ p M(t+Δt) is stored in the storage device 15. The spatial partial differential update equation can be calculated using the Runge-Kutta method or the like.

[0054] Steps S3 to S5 are executed for all voxels of the inspection target model. In addition, for each voxel, the processing circuit 11 calculates the magnetization value at any position within the voxel using the spatial partial differential value under the assumption that the magnetic field within the voxel is linear.

[0055] When step S5 is performed, the processing circuit 11 determines whether or not to end the update calculation by implementing the numerical calculation function 112 (step S6). It is determined that the update calculation is to be ended when a predetermined update end condition is satisfied. The update end condition may be set, for example, to the number of updates reaching a predetermined value. If it is determined that the update calculation is not to be ended (step S6: NO), the processing circuit 11 repeats the processing from steps S3 to S6 for the next timing step until it is determined in step S6 that the update calculation is to be ended.

[0056] When it is determined that the update calculation is to be ended (step S6: YES), the processing circuit 11 outputs the magnetization value M(t+Δt) and the spatial partial differential value ∂ p M(t+Δt) is output (step S7).

[0057] When step S7 is performed, the processing circuit 11 determines whether or not to end the simulation by implementing the numerical calculation function 112 (step S8). It is determined that the simulation is to end when a predetermined simulation end condition is satisfied. The simulation end condition may be set, for example, to completion of processing for all sequence fragments included in the pulse sequence to be calculated. If it is determined that the simulation end condition is not satisfied (step S8: NO), the processing circuit 11 proceeds to step S1 and inputs the next sequence fragment by implementing the input function 111. Then, the processing circuit 11 executes steps S2 to S8 for the sequence fragment, and repeats steps S1 to S8 until it is determined that the simulation end condition is satisfied. As a result, the magnetization vector M(t) for each voxel and each isochromat included in the inspection target model by the pulse sequence at each predetermined time is obtained.

[0058] If it is determined that the simulation end condition is satisfied (step S8: YES), the processing circuitry 11 calculates simulation data corresponding to the simulation objective by implementing the numerical calculation function 112 (step S9). The simulation data may be set to the magnetization vector M(t), IQ signal, and MR image. The IQ signal means k-space data acquired by the ADC. The IQ signal is the x-component ΣM of the sum ΣM of multiple magnetization vectors M(T) corresponding to multiple isochromats. X and the y component ΣM Y MR images are generated by performing reconstruction operations such as Fourier transform on the IQ signals.

[0059] When step S9 is performed, the processing circuit 11 outputs the simulation data output in step S9 (step S10) by implementing the output function 113. The processing circuit 11 can output the simulation data to any output destination such as the communication device 12, the display device 13, the storage device 15, etc.

[0060] With the above, the magnetic resonance numerical simulation according to the first embodiment is completed.

[0061] [Example 2] Second embodiment Next, a magnetic resonance numerical simulation according to a second embodiment will be described.

[0062] The magnetization update equation, the spatial partial differential update equation, and the angular frequency partial differential update equation according to the second embodiment include a rotation term and a relaxation term. The processing circuit 11 according to the second embodiment, by implementing the numerical calculation function 112, refers to the magnetization value before the update, calculates the rotation term of the magnetization update equation to calculate the magnetization rotation term value, and applies the relaxation term of the magnetization update equation to the magnetization rotation term value to calculate the updated magnetization value. Similarly, the processing circuit 11, by implementing the numerical calculation function 112, refers to the magnetization value before the update and the partial differential value in the calculation target direction, calculates the rotation term of the spatial partial differential update equation and / or the angular frequency partial differential update equation to calculate the partial differential rotation term value, and applies the relaxation term of the spatial partial differential update equation to the partial differential rotation term value to calculate the partial differential value in the calculation target direction after the update.

[0063] The spatial partial differential update formula has a first term that depends on an effective magnetic field vector component in a direction perpendicular to the slice selection gradient magnetic field direction, a second term that depends on an effective magnetic field vector component in the slice selection gradient magnetic field direction, and a third term that does not depend on the effective magnetic field vector component, which are separated from each other. The processing circuit 11 according to the second embodiment calculates a rotation term of the first term by referring to the magnetization value before the update and the spatial partial differential value in the calculation target direction to calculate a first rotation term value, calculates a rotation term of the second term by referring to the magnetization value before the update and the spatial partial differential value in the calculation target direction to calculate a second rotation term value, and calculates a rotation term of the third term by referring to the magnetization value before the update and the spatial partial differential value in the calculation target direction to calculate a third rotation term value. The processing circuit 11 then applies the relaxation term of the spatial partial differential update formula to the sum of the first rotation term value, the second rotation term value, and the third rotation term value to calculate an updated spatial partial differential value in the calculation target direction.

[0064] First, an example of an update method according to the second embodiment will be specifically described. Note that all the solutions to the differential equations shown below can also be calculated by repeatedly updating in units of a small time using, for example, the Newton method or the Runge-Kutta method, and all the equations showing the update method based on the differential equations are examples of the solution methods.

[0065] The spatial partial differentials of the Bloch equations shown in (5), (6), and (7) above are approximately updated for small time periods using separation of variables.

[0066] B x (t) and B y (t)-dependent term, B z If we consider updating terms that depend on (t) and other terms in sequence, for example, M x Equation (5), which expresses the spatial partial derivative of the component, can be separated into the following equations (8), (9), and (10). Note that M y Ingredients, M z Formulas for the components can be generated similarly.

[0067]

number

[0068] Spatial RF pulse fluctuation ∂B x (t) / ∂p, ∂B y If (t) / ∂p is considered to be zero, it can be expressed as the following equation (11). This differential equation represents a rotation and can be updated using the Rodrigues rotation formula. M y Ingredients and M z The spatial partial differential update equation based on equation (8) or (11) is x (t) and B y (t), and is an example of the above-mentioned "first term that depends on the effective magnetic field vector component in the direction perpendicular to the slice selection gradient magnetic field direction." Hereinafter, the spatial differential update formula is expressed as B xy This is called the dependent spatial partial differential update equation.

[0069]

number

[0070] Equation (11) is based on the Rodrigues rotation formula, which ignores relaxation of equation (2), and is similar to the update formula described in Non-Patent Document 1, ∂M x , ∂M y , ∂M z For the vector that combines these, the magnetic field B x , B y , B z The rotation matrix can be updated by multiplying it by the corresponding 3x3 rotation matrix.

[0071] M x Ingredients and M y The spatial partial differential update equation based on equations (9) and (10) for the components is expressed by equation (12) for a short time. z (t) bar indicates B that changes within a short time Δt z It means the value when (t) is considered as a specific fixed value. M zThe spatial partial differential update equation based on equations (9) and (10) for the components is expressed by equation (13) for a small time.

number

[0072] Here, M xy =M x +iM y However, if there is no formula (8), B z If (t) is given by a linear expression for x, y, z, and ω, M xy The update formula for M z The update equation for can be applied even if time is not small.

[0073] The spatial partial differential update equation expressed by equation (12) is B z (t), and is an example of the above-mentioned "second term depending on the effective magnetic field vector component in the slice selection gradient magnetic field direction." Hereinafter, the spatial differential update formula is expressed as B z The spatial partial differential update equation expressed by equation (13) is called the spatial partial differential update equation of B x (t), B y (t), B z This is a term that does not depend on (t) and is an example of the above-mentioned "third term that does not depend on the effective magnetic field vector component." Hereinafter, this spatial differential update formula will be referred to as the B-independent spatial partial differential update formula.

[0074] The magnetization vector update formula can be expressed by the following formula (14) using the Rodrigues rotation formula. Formula (14) corresponds to formula (5) in Non-Patent Document 2.

[0075]

number

[0076] The magnetization vector M(t) on the right side of equation (14) represents the magnetization vector of the isochromat contained in the region at position r at time t among the multiple hydrogen atoms placed in the model to be inspected in the magnetic resonance numerical simulation. RF denotes the rotation terms about the z-axis and x-axis due to the application of an RF pulse that produces a flip angle α during a time interval Δt. That is, the rotation term R RF shows the effect of rotation on the magnetization vector M(t). R on the right hand side of Eq. (14) relax is the relaxation term that indicates the effect of magnetic relaxation on the magnetization vector M(t) during the time interval Δt. z (θ) indicates the rotation term about the z-axis with respect to the angle θ. Here, θ g is the angle related to the integral of the gradient field over the time interval Δt. i denotes the angle associated with the inhomogeneity of the static magnetic field over a time interval Δt.

[0077] By spatially differentiating equation (14) in the x, y, and z directions, respectively, it is possible to obtain the spatial partial differential update equations in the x, y, and z directions (Rodrigues rotation formula version of the spatial partial differential update equation) corresponding to equations (5), (6), and (7), respectively. By applying variable separation to the x-direction partial differential update equation, y-direction partial differential update equation, and z-direction partial differential update equation in the Rodrigues rotation formula version, B based on equation (8) or (11) can be obtained. xy The Rodrigues rotation formula version of the spatial partial differential update formula for dependence, expressed by equation (12), B z It is possible to obtain the Rodrigues rotation formula version of the B-dependent spatial partial differential update equation and the Rodrigues rotation formula version of the B-independent spatial partial differential update equation expressed by equation (13).

[0078] Fig. 8 is a diagram illustrating a processing procedure of a magnetic resonance numerical simulation according to Example 2 of the first embodiment. Fig. 8 shows the processing from step S3 to step S6 in Fig. 6. Each processing from step SA1 to SA6 is performed by the numerical calculation function 112. Note that in the example shown in Fig. 8, the angular frequency directional partial differential is not taken into consideration.

[0079] When step S3 is performed, the processing circuit 11 refers only to M(t) and M0 and calculates the rotation term Rot of the magnetization update equation. z (θ g ), Rot z (θ g ) and R RF The magnetization rotation term value is calculated (step SA1). By referring to only M(t) and M0, the rotation term Rot z (θ g ), Rot z (θ g ) and R RF The above calculation can be performed efficiently.

[0080] When step SA1 is performed, the processing circuit 11 adds the relaxation term R of the magnetization update formula to the magnetization rotation term value calculated in step SA1. relax The updated magnetization value M(t+Δt) is calculated by applying the above formula (step SA2). The updated magnetization value M(t+Δt) is stored in the storage device 15.

[0081] When step SA2 is performed, the processing circuit 11 calculates M(t) and ∂ p B with reference only to M(t) xy Rotation term Rot in the spatial partial differential update equation of dependence z (θ g ), Rot z (θ g ) and R RF Then, the first rotation term value is calculated (step SA3). p By referring only to M(t), B xy Rotation term Rot in the spatial partial differential update equation of dependence z (θ g ), Rot z (θ g) and R RF The above calculation can be performed efficiently.

[0082] When step SA3 is performed, the processing circuit 11 calculates M(t) and ∂ p B with reference only to M(t) z Rotation term Rot in the spatial partial differential update equation of dependence z (θ g ), Rot z (θ g ) and R RF Then, the second rotation term value is calculated (step SA4). p By referring only to M(t), B z Rotation term Rot in the spatial partial differential update equation of dependence z (θ g ), Rot z (θ g ) and R RF The calculation of B can be performed efficiently. z If the (t) bar has not yet been calculated, the B z (t) After calculating the bar, the rotation term Rot z (θ g ), Rot z (θ g ) and R RF It is advisable to perform the following calculation.

[0083] When step SA4 is performed, the processing circuit 11 calculates M(t) and ∂ p Rotation term Rot of the spatial partial differential update equation that is independent of B, referring only to M(t) z (θ g ), Rot z (θ g ) and R RF Then, the third rotation term value is calculated (step SA5). p By referring only to M(t), the rotation term Rot of the B-independent spatial partial differential update equation is z (θ g ), Rot z (θ g ) and R RF The above calculation can be performed efficiently.

[0084] After step SA5, the processing circuit 11 adds the first rotation term value calculated in step SA3, the second rotation term value calculated in step SA4, and the third rotation term value calculated in step SA4 (step SA6). Step SA6 calculates the sum of the first rotation term value, the second rotation term value, and the third rotation term value.

[0085] When step SA6 is performed, the processing circuit 11 adds the relaxation term R of the partial differential update formula to the added value obtained in step SA6. relax The spatial partial differential value ∂ after updating by applying p M(t+Δt) is calculated (step SA7). The updated spatial partial differential value ∂ p M(t+Δt) is stored in the storage device 15.

[0086] When step SA6 is performed, step S6 shown in FIG. 6 is executed.

[0087] As described above, according to the second embodiment, even in the update equation described using the Rodrigues rotation formula, the magnetization update equation, the spatial partial differential update equation, and each frequency partial differential update equation can be separated, and the magnetization value and the spatial partial differential value can be efficiently calculated using the contributing terms without using the non-contributing terms. Also, in the spatial partial differential update equation, the spatial partial differential update equation can be separated for each calculation target direction, so that the spatial partial differential value in each calculation target direction can be efficiently calculated using the contributing terms without using the non-contributing terms.

[0088] The processing order of steps SA3, SA4, and SA5 is not limited to this, and may be any order. xy Dependent spatial partial differential update equation (first term), B zIt is not necessary to calculate all of the B-dependent spatial partial differential update equation (second term) and the B-independent spatial partial differential update equation (third term), and the calculation targets may be appropriately selected as necessary. For example, the processing circuit 11 may calculate the first term, the second term, and the third term during the application period of the RF pulse, and calculate the second term and the third term during the non-application period of the RF pulse. As another example, the processing circuit 11 may calculate the first term, the second term, and the third term during both the application period and the non-application period of the RF pulse.

[0089] In addition, when simulating ADC sampling using the partial differentials in the x, y, and z directions obtained, for example, according to the method described in Non-Patent Document 2, a model in which molecules are uniformly placed in a unit rectangular parallelepiped is used, and M(t), ∂ x M(t), ∂ y Using M(t), the angle (∂ / ∂ p )φ is calculated by equation (16) in Non-Patent Document 2, and the spatial variation of the magnetic field can be reflected in the sampled values ​​by multiplying the integral value of the sinc function. When simulating T2 star damping using partial differential in the ω direction, the method described in Non-Patent Document 2 (where molecules are uniformly placed only in a specific frequency region) may be used, or, for example, as another method, a model in which molecules are placed with a probability according to a Lorentzian distribution for frequency may be used, and the partial differential in the ω direction may be multiplied by exp{(-1 / T2')·(∂ / ∂ω)φ)}. T2' is the T2 star damping (T2 * ) and T2 decay, 1 / T2´=(1 / T2 * )-(1 / T2).

[0090] Second embodiment The processing circuit 11 according to the second embodiment uses a combined transition matrix as an update equation describing the time evolution of magnetic resonance. The combined transition matrix is ​​a transition matrix T(t+Δt k-1) (K≧2), where k is a natural number equal to or greater than 1. The K transition matrices T(t), T(t+Δt1), , T(t+Δt K-1 ) combined transition matrix T combined can be expressed by the following equation (15). k represents the time interval from kt to (k+1)t.

[0091]

number

[0092] joint transition matrix T combined The update equation for the magnetization vector M(t) using can be expressed as follows (16).

[0093]

number

[0094] As shown in Fig. 4, the bond transition matrix T combined contains many elements, the number of which is the square of the number of components of the magnetization vector. combined Each element of T a,b is expressed as the bond transition matrix T combined can be expressed by the following equation (17): In equation (17), the number of components of the magnetization vector is set to “n”.

[0095]

number

[0096] When using a combined transition matrix, K multiplications of a transition matrix, whose elements are equal to the square of the number of components of the magnetization vector, are repeated for each voxel. When the inspection target model is three-dimensional, the number of voxels often exceeds one million, and the calculation load, including the amount of calculation and memory usage, is significant.

[0097] The magnetic resonance numerical simulation device according to the second embodiment applies the method of referring to only the contribution terms according to the first embodiment to the coupling transition matrix to speed up the magnetic resonance numerical simulation. The magnetic resonance numerical simulation device according to the second embodiment will be described below. In the following description, components having substantially the same functions as those in the first embodiment will be given the same reference numerals and will be described repeatedly only when necessary.

[0098] FIG. 9 is a diagram showing a configuration example of a magnetic resonance numerical simulation device 1 according to the second embodiment. As shown in FIG. 9, the processing circuit 11 of the magnetic resonance numerical simulation device 1 realizes an input function 111, a numerical calculation function 112, an output control function 113, a determination function 114, and a matrix calculation function 115. Each of the functions 111 to 115 is not limited to being realized by a single processing circuit. A processing circuit may be configured by combining a plurality of independent processors, and each processor may execute a physical simulation program to realize each of the functions 111 to 115. In addition, each of the functions 111 to 115 may be implemented as a module that constitutes a magnetic resonance numerical simulation program, or may be implemented as individual hardware.

[0099] By implementing the numerical calculation function 112, the processing circuit 11 applies the magnetization value and spatial partial differential value before the update to the bonding transition matrix to calculate the updated magnetization value and spatial partial differential value. The bonding transition matrix is ​​assumed to be calculated in advance by the matrix calculation function 115 and stored in the storage device 15. Since the update formula parameters of each matrix element of the bonding transition matrix differ for each sequence fragment, the bonding transition matrix is ​​calculated and stored for each sequence fragment. As an example, the bonding transition matrix may be stored in association with an identifier representing the sequence fragment.

[0100] By implementing the determination function 114, the processing circuit 11 determines whether or not a joint transition matrix corresponding to the sequence fragment to be calculated is stored in the storage device 15. If it is determined that the joint transition matrix is ​​stored, the processing circuit 11 reads out the joint transition matrix corresponding to the sequence fragment to be calculated from the storage device 15 and uses it.

[0101] By implementing the matrix calculation function 115, the processing circuit 11 calculates each matrix element of the bonding transition matrix based on an output vector obtained by applying an input vector including a predetermined magnetization value and a spatial partial differential value in each of a plurality of spatial directions to an update formula. The bonding transition matrix is ​​stored in the storage device 15. As an example, when the determination function 114 determines that a bonding transition matrix corresponding to a sequence fragment to be calculated is not stored in the storage device 15, the processing circuit 11 may calculate and use a bonding transition matrix corresponding to the sequence fragment. Note that, if the sequence is known in advance, the bonding transition matrix may be calculated in advance before the simulation. In addition, the sequence fragment to be calculated for the bonding transition matrix may be cut out based on, for example, the presence or absence of RF irradiation as a criterion, but is not limited to this criterion and may be a sequence cut out based on any criterion. In addition, it is not necessary to calculate the bonding transition matrix for all sequence fragments, and the bonding transition matrix may be calculated only for some sequence fragments.

[0102] The magnetic resonance numerical simulation apparatus 1 according to the second embodiment will be described in detail below.

[0103] First, the matrix elements contained in the bond transition matrix will be described.

[0104] A first matrix element of the coupling transition matrix that contributes to the computation of the time evolution of magnetic resonance has a non-zero value, and a second matrix element of the coupling transition matrix that does not contribute to the computation of the time evolution of magnetic resonance has a zero value. Specifically, the first matrix element is a matrix element corresponding to a combination of a certain component of a magnetization vector and a component that interacts with the certain component, and the second matrix element is a matrix element corresponding to a combination of a certain component of a magnetization vector and a component that does not interact with the certain component.

[0105] More specifically, M x , M y , M z is M x , M y , M z and does not interact with any other components than M0. Therefore, Mx , M y , M z For each non-zero matrix element, M x , M y , M z and M0. ∂ p M is ∂ q M(p≠q) and M0 do not interact. Therefore, ∂ p M x , ∂ p M y , ∂ p M z For each non-zero matrix element, M x , M y , M z , ∂ p M x , ∂ p M y , ∂ p M z In addition, since M0 is a constant, all rows related to M0 can be excluded.

[0106] The processing circuit 11 according to the second embodiment does not perform the calculation of the second matrix elements having a zero value. That is, the second matrix elements having a zero value are not incorporated into the calculation formula of the coupling transition matrix, and the calculation formula is composed of the first matrix elements having a non-zero value. This makes it possible to omit the calculation of many second matrix elements having a zero value, thereby improving the calculation speed.

[0107] The combined transition matrix is ​​stored in the storage device 15, but the storage device 15 stores the first matrix elements having a non-zero value and does not store the second matrix elements having a zero value, thereby making it possible to reduce the memory capacity required to hold the combined transition matrix.

[0108] Next, the calculation of the bond transition matrix by the matrix calculation function 115 will be described.

[0109] In the absence of diffusion terms, each ∂ p M is ∂ qIt is independent of M (p ≠ q). The joint transition matrix can be calculated based on the output vector corresponding to the input vector. The input vector includes a given magnetization value and a spatial partial derivative value in each of multiple spatial directions. The output vector is the magnetization vector obtained by applying the joint transition matrix to the input vector. Note that M0 is invariant in each simulation.

[0110] The bonding transition matrix can be easily obtained by preparing all independent unit vectors for the input vectors and executing the calculations described in the first embodiment. Alternatively, since the bonding transition matrix has many zero elements, some of the calculations can be omitted. For example, if the ω direction is not taken into consideration, the bonding transition matrix can be obtained by calculating only the following seven items. Note that if the ω direction is taken into consideration, ∂ ω Just add a term corresponding to

[0111] V1.M x =1, remaining 0 V2.M y =1, remaining 0 V3.M z =1, remaining 0 V4.∂ x M x =∂ y M x =∂ z M x =1, remaining 0 V5.∂ x M y =∂ y M y =∂ z M y =1, remaining 0 V6.∂ x M z =∂ y M z =∂ z M z =1, remaining 0 V7.All 0

[0112] The calculations described in the first embodiment are performed for each of V1 to V7. The elements of the coupling transition matrix corresponding to the update formula for M0 are known. If zero elements are used, M0 can be updated from V1, V2, V3, and V7. x, M y , M z , we obtain the elements of the joint transition matrix corresponding to the update equation for M0. Similarly, we use the zero elements to obtain p By utilizing the fact that there is no mutual interference between V1, V2, V3, V4 to V6, p From the term corresponding to ∂ p The elements of the bond transition matrix corresponding to

[0113] When using the Rodrigues rotation matrix, it is possible to calculate the bond transition matrix by sharing a part of the intermediate calculation for each calculation of steps SA1, SA3, and SA4 shown in FIG. 8 of Example 2 of the first embodiment. Alternatively, it is possible to calculate by explicitly writing down the matrix. When implementing the software using the initial values, it is possible to shorten the time required to calculate the bond transition matrix by simultaneously executing simulations for seven initial vectors for each isochromat. When obtaining the bond transition matrix by using the matrix product of the transition matrices, it is possible to shorten the time required to calculate the bond transition matrix by dividing the multiplication by each submatrix. In either case, it is possible to shorten the time required to calculate the bond transition matrix by sharing the values ​​of sin, cos, and exp as much as possible.

[0114] [Example 1] Fig. 10 is a diagram illustrating a processing procedure of a magnetic resonance numerical simulation according to Example 1 of the second embodiment. Fig. 10 shows the processing from step S3 to step S5 in Fig. 6. Each processing from step SB1 to SB2 is performed by the numerical calculation function 112.

[0115] As shown in FIG. 10, when step S3 is performed, the processing circuit 11 calculates the joint transition matrix T combined (t) is read from the storage device 15 (step SB1).

[0116] When step SB1 is performed, the processing circuit 11 calculates M(t), ∂ pM(t) and M0 are the joint transition matrix T combined (t), the updated M(t+Δt), ∂ p M(t+Δt) is calculated (step SB2). In step SB2, the processing circuit 11 does not perform the calculation of the second matrix element having a zero value, but performs the calculation of the first matrix element having a non-zero value. This makes it possible to reduce the calculation time of the bond transition matrix.

[0117] After step SB2 is performed, step S6 shown in FIG. 6 is executed.

[0118] [Example 2] 11 is a diagram illustrating a processing procedure of a magnetic resonance numerical simulation according to Example 2 of the second embodiment. FIG 11 shows the processing from step S3 to step S5 in FIG 6.

[0119] As shown in FIG. 11, when step S3 is performed, the processing circuit 11 realizes the determination function 114 to determine the joint transition matrix T combined It is determined whether (t) is stored in the storage device 15 (step SC1).

[0120] In step SC1, the bond transition matrix T combined If it is determined that (t) is not stored (step SC1: NO), the processing circuit 11 performs the matrix calculation function 115 to calculate the joint transition matrix T combined (t) is calculated (step SC2).

[0121] When step SC2 is performed, the processing circuit 11 realizes the output control function 113 to convert the coupling transition matrix T combined (t) is cached in the storage device 15 (step SC3). combined (t) is stored in association with the identifier of the sequence fragment input in step S1.

[0122] When step SC3 is performed or when step SC1 is performed, the coupling transition matrix T combined If it is determined that (t) is stored (step SC1: YES), the processing circuit 11 stores the joint transition matrix T combined (t) is read out from the storage device 15 (step SC4).

[0123] When step SC4 is performed, the processing circuit 11 realizes the numerical calculation function 112 to calculate M(t), ∂ p M(t) and M0 are the joint transition matrix T combined (t), the updated M(t+Δt), ∂ p M(t+Δt) is calculated (step SC5).

[0124] After step SC5, step S6 shown in FIG. 6 is executed.

[0125] According to the second embodiment, each time an unknown sequence fragment is detected, the combined transition matrix T combined By calculating (t) and reusing the results while dynamically caching them, the magnetic resonance numerical simulation can be accelerated.

[0126] (Modification) The above-described processing by the magnetic resonance numerical simulation apparatus 1 is one example, and various additions, modifications and / or deletions are possible without changing the gist of this embodiment.

[0127] [Variation 1] The processing circuit 11 according to the above embodiment uses all pre-update values ​​belonging to the contribution terms among the pre-update values ​​of the magnetization and the pre-update values ​​of the partial differential to calculate the update equation for the calculation target, and calculates the post-update values ​​of the calculation target. However, this embodiment is not limited to this. The processing circuit 11 does not need to use all pre-update values ​​belonging to the contribution terms, and may use only some of the pre-update values ​​belonging to the contribution terms. As an example, the processing circuit 11 according to the first modification may not use values ​​that can be considered to be zero or approximately zero among the pre-update values ​​belonging to the contribution terms.

[0128] [Variation 2] When calculating a predicted output value of an A / D converter (ADC), the processing circuit 11 according to the second modification does not calculate an updated partial differential value in a calculation target direction for a magnetization component in a slice selection gradient magnetic field direction (z direction), but calculates an updated partial differential value in a calculation target direction for a magnetization component in a gradient magnetic field direction other than the slice selection gradient magnetic field direction. For example, when calculating a predicted output value of an A / D converter (ADC), the processing circuit 11 according to the second modification calculates ∂ p M z ∂ p M x and ∂ p M y Only ∂ p M z This makes it possible to reduce the calculation time for calculating the predicted output value of an A / D converter (ADC).

[0129] The processing circuit 11 may perform ADC calculations in a continuous system using spatial partial derivatives. This can be performed by a calculation similar to that of Non-Patent Document 2 using a diffusion term.

[0130] According to at least one of the embodiments described above, it is possible to improve the calculation efficiency of the magnetic resonance numerical simulation.

[0131] The term "processor" used in the above description means a circuit such as a CPU, a GPU, an application specific integrated circuit (ASIC), a programmable logic device (e.g., a simple programmable logic device (SPLD), a complex programmable logic device (CPLD), and a field programmable gate array (FPGA)). The processor realizes a function by reading and executing a program stored in a memory circuit. Instead of storing a program in a memory circuit, the program may be directly built into the circuit of the processor. In this case, the processor realizes a function by reading and executing a program built into the circuit. On the other hand, when the processor is, for example, an ASIC, instead of storing the program in a memory circuit, the function is directly built into the circuit of the processor as a logic circuit. Note that each processor in this embodiment is not limited to being configured as a single circuit for each processor, but may be configured as a single processor by combining multiple independent circuits. Furthermore, the multiple components in FIG. 1 and FIG. 9 may be integrated into a single processor to realize the function.

[0132] Although some 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 various omissions, substitutions, modifications, and combinations of embodiments can be made without departing from the spirit of the invention. These embodiments and their modifications are included in the scope of the invention and its equivalents as described in the claims, as well as in the scope and spirit of the invention. [Explanation of symbols]

[0133] 1. Magnetic resonance numerical simulation device 11 Processing circuit 12. Communications equipment 13 Display equipment 14 Input Devices 15 Storage device 41 Non-applied period 42 Application Period 111 Input Function 112 Numerical Calculation Function 113 Output control function 114 Judgment function 115 Matrix calculation function

Claims

1. an input unit for inputting a value of magnetization before updating and a value of partial differential of the magnetization before updating in each of a spatial direction and / or an angular frequency direction; a numerical calculation unit that uses some of the values ​​of the magnetization before the update and the partial differential before the update to calculate an update equation that represents the time evolution of magnetic resonance, and calculates, for the isochromat, the value of the magnetization after the update and the value of the partial differential after the update in a calculation target direction among the spatial direction and / or the angular frequency direction with respect to the magnetization; A magnetic resonance numerical simulation device comprising:

2. The numerical calculation unit is Calculating the update equation using the value of the magnetization before the update without using the value of the partial differential before the update, and calculating the value of the magnetization after the update; calculating the update formula using the value of the partial differential before the update in the direction to be calculated, without using a value of the partial differential before the update in a direction that is not a target of calculation among the spatial direction and / or the angular frequency direction, and calculating a value after the update of the partial differential in the direction to be calculated; The magnetic resonance numerical simulation apparatus according to claim 1.

3. the update equations include a magnetization update equation, a spatial partial differential update equation, and an angular frequency partial differential update equation that are separated from each other, The numerical calculation unit is Calculating the magnetization update equation by referring to the value of the magnetization before the update in order to calculate the updated value of the magnetization; calculating the spatial partial differential update formula by referring to a pre-update value of the magnetization and a pre-update value of the partial differential of the spatial direction in the calculation target direction in order to calculate an updated value of the partial differential of the spatial direction in the calculation target direction; calculating an angular frequency partial differential update formula by referring to a value before the update of the magnetization and a value before the update of the partial differential in the angular frequency direction in order to calculate an updated value of the partial differential in the angular frequency direction; The magnetic resonance numerical simulation device according to claim 2.

4. The numerical calculation unit is calculating a rotation term value of the magnetization update equation by referring to a value of the magnetization before the update; A relaxation term in the magnetization update equation is applied to the magnetization rotation term value to calculate an updated value of the magnetization. The magnetic resonance numerical simulation apparatus according to claim 3.

5. The numerical calculation unit is calculating a rotation term in the spatial partial differential update equation and / or the angular frequency partial differential update equation by referring to the pre-update value of the magnetization and the pre-update value of the partial differential in the calculation target direction to calculate a partial differential rotation term value; calculating an updated value of the partial differential in the calculation target direction by applying a relaxation term in the spatial partial differential update equation and / or the angular frequency partial differential update equation to the partial differential rotation term value; 5. The magnetic resonance numerical simulation apparatus according to claim 4.

6. 6. The magnetic resonance numerical simulation device according to claim 5, wherein the spatial partial differential update equation has a first term that depends on an effective magnetic field vector component in a direction perpendicular to a slice selection gradient magnetic field direction, a second term that depends on an effective magnetic field vector component in the slice selection gradient magnetic field direction, and a third term that does not depend on the effective magnetic field vector component, which are separated from each other.

7. The numerical calculation unit is Calculating the first term, the second term, and the third term during an application period of an RF pulse; The second term and the third term are calculated during a period in which an RF pulse is not applied.

7. The magnetic resonance numerical simulation apparatus according to claim 6.

8. 7. The magnetic resonance numerical simulation apparatus according to claim 6, wherein said numerical calculation unit calculates said first term, said second term and said third term during an application period and a non-application period of an RF pulse.

9. The numerical calculation unit is to calculate an updated value of the spatial partial differential in the direction to be calculated, a rotation term of the first term is calculated by referring to the value of the magnetization before the update and the value of the spatial partial differential in the direction to be calculated to calculate a first rotation term value, a rotation term of the second term is calculated by referring to the value of the magnetization before the update and the value of the spatial partial differential in the direction to be calculated to calculate a second rotation term value, and a rotation term of the third term is calculated by referring to the value of the magnetization before the update and the value of the spatial partial differential in the direction to be calculated to calculate a third rotation term value, calculating an updated value of the spatial partial differential in the calculation target direction by applying a relaxation term in the spatial partial differential update equation to a sum of the first rotation term value, the second rotation term value, and the third rotation term value; 7. The magnetic resonance numerical simulation apparatus according to claim 6.

10. 3. The magnetic resonance numerical simulation apparatus according to claim 2, wherein the update equation is a coupled transition matrix that combines a plurality of transition matrices that describe the time evolution of the magnetization and the time evolution of the partial derivative for each of a plurality of time steps.

11. A first matrix element of the coupling transition matrix that contributes to the calculation of the time evolution of magnetic resonance has a non-zero value, A second matrix element of the coupling transition matrix that does not contribute to the calculation of the time evolution of magnetic resonance has a zero value, The numerical calculation unit does not perform an operation on the second matrix element having the zero value. The magnetic resonance numerical simulation device according to claim 10.

12. A storage unit that stores the bond transition matrix, the storage unit stores the first matrix elements and does not store the second matrix elements; The magnetic resonance numerical simulation device according to claim 11.

13. 11. The magnetic resonance numerical simulation device according to claim 10, further comprising a matrix calculation unit that calculates each matrix element of the coupling transition matrix based on an output vector obtained by applying an input vector including a predetermined magnetization value and a spatial partial differential value of the magnetization in each of the spatial directions to the update equation.

14. A storage unit that stores a bond transition matrix for each sequence fragment of a pulse sequence; a determination unit that determines whether a joint transition matrix corresponding to a sequence fragment to be calculated is stored in the storage unit, The numerical calculation unit is When it is determined by the determination unit that the matrix is ​​stored, the combination transition matrix corresponding to the sequence fragment to be calculated is read out from the storage unit and used; If it is determined by the determination unit that the sequence fragment is not stored, a bond transition matrix corresponding to the sequence fragment to be calculated is calculated, and the calculated bond transition matrix is ​​used. The magnetic resonance numerical simulation device according to claim 10.

15. 2. The magnetic resonance numerical simulation device according to claim 1, wherein, when calculating a predicted output value of an A / D converter, the numerical calculation unit does not calculate an updated value of the partial differential in the calculation target direction for the magnetization component in a slice selection gradient magnetic field direction, but calculates an updated value of the partial differential in the calculation target direction for the magnetization component in a gradient magnetic field direction other than the slice selection gradient magnetic field direction.

16. For an isochromat, a value of magnetization before updating and a value of partial differential of the magnetization in each of a spatial direction and / or an angular frequency direction before updating are input; Calculate an update equation that represents the time evolution of magnetic resonance using some of the values ​​of the magnetization before the update and the partial differential before the update, and calculate, for the isochromat, the value of the magnetization after the update and the value of the partial differential after the update in a calculation target direction among the spatial direction and / or the angular frequency direction with respect to the magnetization. The magnetic resonance numerical simulation method includes:

17. On the computer, A function for inputting a value of magnetization before updating and a value of partial differential of the magnetization in each of a spatial direction and / or an angular frequency direction before updating, for an isochromat; a function of calculating an update equation expressing the time evolution of magnetic resonance by using some of the values ​​of the magnetization before the update and the partial differential before the update, and calculating, for the isochromat, the value of the magnetization after the update and the value of the partial differential after the update in the direction to be calculated among the spatial direction and / or the angular frequency direction with respect to the magnetization; A magnetic resonance numerical simulation program that realizes this.

Images