Patent Application
20180059198
The present invention relates to a magnetic resonance imaging (MRI) technique. More particularly, it relates to an image processing technique that utilizes an acquired image to calculate a magnetic susceptibility distribution within a living body.
A magnetic resonance imaging apparatus is a medical-use diagnostic noninvasive imaging apparatus, utilizing nuclear magnetic resonance phenomenon, which is a phenomenon that hydrogen nucleus (protons) placed in a static magnetic field are resonant with an RF magnetic field at a specific frequency. Since nuclear magnetic resonance signals are changed depending on various physical properties, such as proton density and relaxation time, an image obtained by the MRI can depict various biological information, such as a structure or a composition of living tissues and cell properties.
In recent years, a magnetic susceptibility difference between living tissues receives attention, as one of the physical properties being measurable by the MRI. The magnetic susceptibility is a physical property that represents a degree of magnetic polarization (magnetization) of materials in the static magnetic field. In a living body, there are contained paramagnetic substances such as deoxyhemoglobin in venous blood and iron protein, and diamagnetic substances such as water constituting a large part of the living tissues and calcium serving as a basis of calcification. Creating an image quantitatively from a difference in magnetic susceptibility, i.e., the magnetic susceptibility difference between the living tissues, may be applicable to diagnosis of cerebral ischemia disease, prediction of radiation treatment effect against cancer, and identification of neurodegenerative disease.
A method for creating an image of the magnetic susceptibility difference between living tissues by utilizing the MRI is referred to as quantitative susceptibility mapping (QSM). The QSM is a method of calculating spatial variation in magnetic field caused by the magnetic susceptibility difference between living tissues, from phase information of an MR image being measured, and obtaining a magnetic susceptibility distribution according to a relational expression between the magnetic field and the magnetic susceptibility.
However, a distribution of the magnetic field is obtained by subjecting the magnetic susceptibility distribution to spatial convolution integral. Therefore, obtaining the magnetic susceptibility distribution from the magnetic field distribution is an inverse problem of the convolution integral, and thus a unique solution cannot be obtained.
In general, the least-square method is employed to obtain the magnetic susceptibility distribution from the magnetic field distribution. In this case, an error function is introduced, and a value that minimizes this error function is obtained as a solution. Representative examples of this method are, for example, TKD (Truncated-based K-space Division) (e.g., see Non Patent Document 1), Iterative SWIM (Susceptibility Weighted Imaging and Mapping) (e.g., see Non Patent Document 2 and Patent Document 1), and MEDI (Morphology enabled dipole inversion) method (e.g., see Non Patent Document 3 and Patent Document 2). Here, the TKD is a method of calculating a magnetic susceptibility distribution according to operations on k-space of the magnetic field distribution and of the point dipole magnetic field; the iterative SWIM is a method of merging through an iterative operation, the magnetic susceptibility distribution calculated according to the TKD method, with the magnetic susceptibility distribution where a fine structure is extracted according to thresholding; and the MEDI is a method that utilizes the regularized least-square method.
In the TKD method, a region for which inverse calculations cannot be performed is truncated at zero or at an optional value. Therefore, the TKD method has an advantage being simple and short in calculation time, but there are problems that artifacts may occur, and accuracy in the calculated magnetic susceptibility may be reduced, due to the truncation in the vicinity of a magic angle region.
The iterative SWIM method applies a masking process allowing the fine structure to be maintained, to the magnetic susceptibility distribution obtained by the TKD method, and iterates the process for updating the magnetic susceptibility distribution. This method allows the magnetic susceptibility distribution to be calculated with keeping the fine structure, but there is a problem that artifacts cannot be removed completely. In order to extract the fine structure, various thresholding processes are required. Therefore, there is a problem that if there are not provided thresholds appropriate for subjects respectively, the fine structure cannot be extracted properly.
The MEDI method calculates the magnetic susceptibility distribution, by minimizing through the iterative operation, a term that provides a smoothing effect, referred to as a regularization term, assuming a relational expression between the magnetic field and the magnetic susceptibility as a constraint. This regularized least-square method is a process that truncates a value of the signal with the magnitude corresponding to a value of a regularization parameter. Therefore, it is necessary to set an adequate regularization term, according to noise that has mixed into the magnetic field distribution. Generally, since it is not possible to obtain in advance noise information mixing into the magnetic field distribution, an adequate regularization parameter is calculated from a change curve of residual which is obtained by varying the regularization parameter more than once. Therefore, in the MEDI method, enormous calculation time is required for setting the adequate regularization parameter.
The present invention has been made in view of the situations as described above, and an object of the present invention is to provide a technique that is used in calculating the magnetic susceptibility distribution with the use of MRI, for reducing artifacts and noise through a simple process, so as to improve accuracy and precision of the magnetic susceptibility distribution being calculated.
According to the present invention, a magnetic field distribution is obtained via a phase image, from a complex image acquired by MRI. Smoothing process is iteratively applied to a magnetic susceptibility distribution that is calculated from the magnetic field distribution, under the constraint based on a relational expression between the magnetic field and the magnetic susceptibility. Specifically, the magnetic susceptibility distribution is calculated according to a minimization process for minimizing an error function through the iterative operation, the error function being defined by a predetermined initial magnetic susceptibility distribution and the magnetic field distribution. On this occasion, the updated magnetic susceptibility distribution calculated according to the error function during the minimizing process, is smoothed every time the updated magnetic susceptibility distribution is calculated.
In calculating the magnetic susceptibility distribution by using the MRI, artifacts and noise are reduced through the simple process, thereby improving accuracy and precision of the magnetic susceptibility distribution being calculated.
An embodiment to which the present invention is applied will now be described. Hereinafter, in all the figures illustrating the embodiment of the present invention, elements with an identical function are labeled with the same reference numeral, unless otherwise noted, and they will not be redundantly explained. It should also be noted that the descriptions below will not limit the present invention.
In the present embodiment, through the use of an MRI apparatus, a complex image I at a predetermined echo time is obtained, a magnetic field distribution δ that captures changes in magnetic field caused by a magnetic susceptibility difference between living tissues, is calculated from the complex image I, and a magnetic susceptibility distribution χ is calculated from the magnetic field distribution δ. The magnetic susceptibility distribution is calculated according to so-called QSM method. Prior to describing the MRI apparatus implementing the steps above, a method of calculating the magnetic susceptibility distribution χ from the complex image I according to the QSM method will be described.
In the QSM method, a local change in magnetic field caused by the magnetic susceptibility difference between living tissues is calculated from a phase image that is obtained by the Gradient Echo(GrE) method.
When a position vector is represented as r, a relative magnetic field distribution δ(r), caused by the magnetic susceptibility difference between tissues, is expressed by the following formula 1.
where P(r) is the phase image, γ is the nuclear gyromagnetic ratio of proton, B0 is static magnetic field strength, and TE is the echo time.
The magnetic field distribution δ(r) is expressed by the following formula 2 according to Maxwell's equations concerning the static magnetic field, by using the magnetic susceptibility distribution χ(r) within a living body.
where α is an angle formed by the vector(r′−r) with the direction of static magnetic field, and d(r) is the point-dipole magnetic field.
As shown by the formula 2, the magnetic field distribution δ(r) is represented by the convolution integral between the magnetic susceptibility distribution χ(r) and the point-dipole magnetic field d(r). Therefore, Fourier transform is applied to both sides of the formula 2, whereby the formula 2 is transformed to the following formula 3.
where k=(kx, ky, kz) represents position vectors in the k-space, Δ(k), X(k), and D(k) are Fourier components, respectively, of the magnetic field distribution δ(r), the magnetic susceptibility distribution χ(r), and the point-dipole magnetic field d(r).
As indicated by the formula 3, the Fourier component X(k) of the magnetic susceptibility distribution can be calculated by dividing the Fourier component Δ(k) of the magnetic field distribution, by the Fourier component D(k) of the point-dipole magnetic field. However, according to the formula 3, since the reciprocals diverge in a domain in the vicinity of D(k)=0, it is not possible to calculate X(k) straightforwardly.
The domain where D(k)=0 is referred to as the magic angle, forming a region of reverse bicone with the apex nearly double the angle of 54.7° with respect to the magnetic field direction. Due to this magic angle, the QSM method which calculates the magnetic susceptibility distribution from the magnetic field distribution comes down to an ill-conditioned inverse problem.
As a representative method for solving this problem, the TKD method, the Iterative SWIM method (hereinafter, simply referred to as “SWIM method”), and the MEDI method are provided as described above.
In the TKD method, the reciprocal of Fourier component D(k) of the point-dipole magnetic field d(r) is obtained, and the vicinity of the magic angle region thereof is truncated at zero or at any value. Thereafter, the obtained value is multiplied by the Fourier component Δ(k) of the measured magnetic field distribution δ(r), and the product is subjected to the inverse Fourier transform to the real space data, thereby calculating the magnetic susceptibility distribution χ(r).
In the Iterative SWIM method, thresholding is applied to the magnetic susceptibility distribution calculated by a method such as the TKD method, and only a fine structure is extracted. Next, the magnetic susceptibility distribution of the fine structure is calculated, which is obtained by masking a part of the magnetic susceptibility other than the fine structure. The original magnetic susceptibility distribution and the magnetic susceptibility distribution of the fine structure are mixed in the k-space, thereby calculating a modified magnetic susceptibility distribution. Thresholding is applied to the modified magnetic susceptibility distribution, and then, the fine structure is extracted again. By repeating a series of those processes above, the magnetic susceptibility distribution is calculated.
In the MEDI method, as shown by the following formula 4, assuming the relational expression between the magnetic field and the magnetic susceptibility, as a constraint, the term referred to as regularization term which provides a smoothing effect, is minimized through the iterative operation, thereby calculating the magnetic susceptibility distribution χ(r).
where χ is the magnetic susceptibility distribution vector to be calculated, δ is the magnetic field distribution vector calculated from the phase image, W is the matrix having a weight vector as a diagonal component which is calculated based on an absolute value image, M is the matrix having a binary edge mask vector as the diagonal component, assuming the edge region as 0 and the other area as 1 on the basis of space gradient of the absolute value image. FD is an operator representing the convolution integral, ∇(nabla) is a three-dimensional space-gradient operator, ∥a∥n is Ln norm of any vector a. The second term λ in the formula 4 is a regularization parameter for striking a balance between the smoothing effect and accuracy of the magnetic susceptibility.
As shown in the formula 4, the MEDI method uses as the regularization term, L1 norm of the vector obtained by multiplying Total Variation (TV) ∇χ of the magnetic susceptibility distribution in the course of calculation, by the edge mask M which is based on the absolute value image. In general, a least-square solution with L1 norm regularization of any vector a, is likely to be a sparse solution. Therefore, according to the MEDI method, a region with less variation in the absolute value image is smoothed, resulting in that the magnetic susceptibility distribution that retains the edge region with a major change is obtained.
It should be understood that the regularized least-square method including the MEDI method is a process for truncating a signal value with a magnitude that corresponds to a value of the regularization parameter. Therefore, if the regularization parameter is too large, the fine structure is also likely to be eliminated along with noise, and further reducing the magnetic susceptibility of the entire image, resulting in that the calculated magnetic susceptibility is deteriorated. On the other hand, if the regularization parameter is too small, noise is likely to be increased, even though the fine structure and the magnetic susceptibility of the entire image are maintained.
Therefore, in the MEDI method, it is necessary to provide an adequate regularization parameter. In the case where the adequate regularization parameter is provided, a standard deviation of a residual of the first term in the formula 4 becomes equal to the standard deviation of the noise that is mixed into the measured magnetic field distribution. In general, however, since it is not possible to obtain in advance, information of the noise mixed into the magnetic field distribution, the regularization parameter is made to vary more than once to calculate residuals, and the adequate regularization parameter is obtained according to a change curve of the residuals thus calculated. Accordingly, the MEDI method requires enormous calculation time for providing the adequate regularization parameter.
In the present embodiment, considering the aforementioned publicly-known methods, an MRI apparatus and QSM method are provided in calculating the magnetic susceptibility distribution by using MRI, so as to obtain highly accurate and highly precise magnetic susceptibility distribution, according to a simple process along with reducing artifacts and noise.
Firstly, an magnetic resonance imaging apparatus (MRI apparatus) of the present embodiment will be described.
The present embodiment may employ any types of the MRI apparatus having those external views as shown above.
It should be noted that these are a few examples, and the MRI apparatus of the present embodiment is not limited to those examples shown here. In the present embodiment, various publicly-known MRI apparatuses may be employed, including any mode or any type thereof.
The MRI apparatus 100 will now be described as a representative example, unless otherwise distinguished. The present embodiment employs, for example, a coordinate system assuming a direction of the static magnetic field of the MRI apparatus 100 as z-direction, and two other directions vertical to the z-direction as x-direction and y-direction; i.e., the direction parallel to a bed surface where a subject to be measured is placed is the x-direction and the other direction is the y-direction. It should also be noted that the static magnetic field will be simply referred to as magnetic field in the following.
Various types of the static magnetic field coil 102 may be employed, depending on the structures of the MRI apparatuses 100, 120, and 130 as shown in
The sequence controller 114 controls operations of the gradient magnetic field power supply 112 being a power source for driving the gradient coil 103, the shim power supply 113 being a power source for driving the shim coil 104, the transmitter 107, and the receiver 108, thereby controlling application of the gradient magnetic field and the RF magnetic field, and receiving timing of the nuclear magnetic resonance signal. A time chart of the control is referred to as a pulse sequence, whose settings are previously configured depending on the measurement, and it is stored, for example, in a storage device being provided in the computer 109 as described below.
The computer 109 controls the entire operation of the MRI apparatus 100, along with performing various computations on the nuclear magnetic resonance signals being received. In the present embodiment, the computer generates a complex image of at least one echo time, a phase image, a magnetic field distribution, a magnetic susceptibility distribution, and the like. The computer 109 is an information processor including a CPU, a memory, and a storage device, and the computer 109 is connected to a display 110, an external storage device 111, and an input device 115.
The display 110 is an interface for an operator, configured to display information such as a result obtained by computations. The input device 115 is an interface configured to allow the operator to input conditions, parameters, and the like, necessary for the computations performed in the present embodiment. The input device 115 of the present embodiment allows the user to input measurement parameters, such as the number of echoes to be measured, a reference echo time, and an echo interval. The external storage device 111, along with the storage device, holds data such as the data used in various computations executed by the computer 109, data obtained through the computations, conditions and parameters being inputted.
As described above, in the present embodiment, under the constraint based on the relational expression between the magnetic field magnetic field and the magnetic susceptibility, a smoothing process is performed through the iterative operation, thereby calculating the magnetic susceptibility distribution, easily with a high degree of precision. Next, a description will be provided regarding a functional configuration of the computer 109 for implementing the aforementioned functions.
The computer 109 is provided with a CPU, a memory, and a storage device. The CPU loads programs held in advance in the storage device into the memory and executes those programs, whereby those various functions are implemented by the computer 109. The storage device or the external storage device 111 may store various data used for processing of the functions, and various data generated in process.
At least one of the various functions implemented by the computer 109 may also be implemented by an information processor that is independent of the MRI apparatus 100 and data transmittable and receivable thereto and therefrom. Furthermore, all or a part of the functions may be implemented by hardware such as ASIC (Application Specific Integrated Circuit) and FPGA (Field-Programmable Gate Array).
There will be described a flow of the process according to the QSM method (QSM process) implemented by the aforementioned functions of the computer 109 of the present embodiment.
Firstly, the measurement controller 310 controls the sequence controller 114 according to a predetermined pulse sequence, thereby measuring an echo signal at a preset echo time (step S1001).
Then, the image reconstructor 320 places thus obtained echo signals in the k-space, applies Fourier transform thereto, thereby reconstructing a complex image I (step S1002).
The magnetic field distribution calculator 330 performs a magnetic field distribution calculating process for calculating the magnetic field distribution δ from the complex image I (step S1003).
Then, the magnetic susceptibility distribution calculator 340 performs a magnetic susceptibility distribution calculating process for calculating the magnetic susceptibility distribution χ from the magnetic field distribution δ being obtained (step S1004), and displays thus obtained magnetic susceptibility distribution χ on the display 110 (step S1005).
When the magnetic susceptibility distribution χ is displayed on the display 110, it is also possible to display an image obtained in the course of the QSM process on the display 110, if required, in addition to the magnetic field distribution δ, and the like, calculated in step S1003.
Details of each step of the QSM process according to the present embodiment will now be described.
The measurement controller 310 activates the sequence controller 114, according to the pulse sequence that is configured on the basis of the parameters inputted by a user via the input device 115, and performs measurement so as to obtain a nuclear magnetic resonance signal (echo signal) at a preset echo time (TE). The sequence controller 114 controls each part of the MRI apparatus 100 according to the instructions from the measurement controller 310, so as to perform the measurement. In the present embodiment, an echo signal of at least one echo time is obtained.
An example of the pulse sequence used by the measurement controller 310 for the measurement will be described. By way of example, in the present embodiment, a pulse sequence of GrE (Gradient Echo) system may be employed. The pulse sequence of this GrE system is sensitive to local magnetic field variation generated by magnetic susceptibility differences between living tissues.
In the RSSG sequence 510, high-frequency (RF) magnetic field pulse 502 is applied along with application of the slice gradient magnetic field pulse 501, thereby exciting magnetization of a predetermined slice within the subject 101. Then, a slice encoding gradient magnetic field pulse 503 and a phase encoding gradient magnetic field pulse 504 are applied, so as to add positional information in the slice-direction and in the phase-encoding direction, to a phase of the magnetization.
After applying a readout gradient magnetic field pulse 505 for dephasing, to disperse a phase of nuclear magnetization within a pixel, one nuclear magnetic resonance signal (echo) 507 is measured, along with applying a readout gradient magnetic field pulse 506 for adding positional information in the readout direction. And finally, a slice encoding gradient magnetic field pulse 508 and a phase encoding gradient magnetic field pulse 509 for rephasing are applied, for convergence of the phase of nuclear magnetization that has been dephased by the slice encoding gradient magnetic field pulse 503 and the phase encoding gradient magnetic field pulse 504.
The measurement controller 310 executes the procedures above, repeatedly every repetition time TR, while varying strength of the slice encoding gradient magnetic field pulses 503 and 508 (slice-encoding count ks) and of the phase encoding gradient magnetic field pulses 504 and 509 (phase-encoding count kp), and the phase of the RF pulse 502, whereby echoes necessary for acquiring one image are measured. In this situation, the phase of the RF pulse 502 may be incremented by 117 degrees, for instance. In
The measured echoes are arranged in three-dimensional k-space having kr, kp, and ks as coordinate axes. In this situation, one echo occupies one line which is parallel to the kr axis in the k-space. As an absolute value image obtained by this RSSG sequence 510, T1 (longitudinal relaxation time) weighted image is obtained when TE is set short, whereas T2* weighted image in which the phase dispersion within the pixel is reflected is obtained when TE is set long.
The RSSG sequence 510 is taken as an example here, which is one method of the Cartesian imaging for acquiring data in parallel to the coordinate axes of the k-space. In the present embodiment, however, any sequence may be employed to acquire data in any k-space domain.
By way of example, another sequence may be employed for acquiring plural echoes at a single TR, and arranging those echoes in various domains of k-space. Alternatively, another sequence may be employed which performs the measurement more than once, and data is obtained in different k-space domains respectively by the plural measurements, and integrates those data items to obtain one image. Further alternatively, non-Cartesian imaging may be employed, such as radial scanning to acquire data in k-space in rotational manner.
It is also possible to employ a sequence for creating one complex image and one phase image, from different complex images and phase images acquired at plural TEs. In the case where a pulse sequence that implements Dixon method for acquiring an image from plural signals that are obtained at plural TEs, the magnetic field distribution δ can be obtained in the process of separating water from fat, for instance. Therefore, in this case, a process for generating the phase image P as described in the following becomes unnecessary.
Next, the image reconstruction process in the aforementioned step S1002 will be described. The image reconstructor 320 places in the k-space the echo signals at the echo time measured in step S1001, and applies Fourier transform to the echo signals. Then, image reconstructor 320 calculates the complex image I where each pixel is a complex number.
Next, there will be described the magnetic field distribution calculating process in the aforementioned step S1003, according to the magnetic field distribution calculator 330. The magnetic field distribution calculator 330 of the present embodiment calculates the magnetic field distribution δ from the complex image I that is reconstructed by the image reconstructor 320. The magnetic field distribution δ being calculated captures local variation of the magnetic field, caused by the magnetic susceptibility differences between living tissues.
The magnetic field distribution calculator 330 firstly calculates the phase image P from the complex image I, and applies various phase image processing to thus calculated phase image P, so as to remove background phase. Then, the magnetic field distribution δ is obtained from the phase image P after performing such phase image processing.
In the phase image processing, a phase component originated from the magnetic susceptibility is extracted, by removing the phase component originated from any element other than the magnetic susceptibility. Specifically, in the calculated phase image P, there are performed a phase unwrapping process for removing a phase exceeding the range from −n to n and folded back, and a global phase change eliminating process for eliminating global phase change caused by a biological form.
In order to implement the processes above, the magnetic field distribution calculator 330 of the present embodiment is provided with, as shown in FIG. 3, a phase image calculator 331 configured to calculate the phase image P from the complex image I, a phase unwrapping processor 332 configured to perform the phase unwrapping process, a global phase change eliminator 333 configured to perform the global phase change eliminating process, and a magnetic field distribution converter 334 configured to convert the phase image P after the phase processing into the magnetic field distribution δ.
The phase image calculator 331 obtains as a pixel value, an argument of each pixel value of the complex image I, so as to calculate the phase image P. In other words, the phase image P is an image having pixel values, which are phase components in the complex numbers of the pixels in the complex image I. Here, the pixel value P(r) of the phase image P at the pixel r is calculated according to the formula 5, by using the pixel value I(r) of the complex image I:
P(r)=arg{I(r)} (5)
The phase unwrapping processor 332 executes the phase unwrapping process on the phase image P, for eliminating the phase exceeding the range from −π to π and folded back. In a partial region of the phase image P, the phase value exceeding the range from −π to π is folded back into the range from −π to π. Therefore, in order to obtain an accurate phase in the entire portion to be imaged (e.g., a head), it is necessary to correct this folding back. In the present embodiment, by using an region growing method, for instance, a phase value folded back into the range from −π to π is corrected.
The global phase change eliminator 333 eliminates a global phase change that occurs due to the biological form, in the phase image P after the phase unwrapping process is performed.
Each pixel value of the phase image P is a sum of the global phase change caused by static magnetic field inhomogeneity that may occur depending on forms of the portion to be imaged, or the like, and a local phase change caused by the magnetic susceptibility change between tissues. The global phase change and the local phase change are associated with, respectively, a low-frequency component and a high-frequency component in the k-space of the phase image P. According to this process, only the global phase change is eliminated, and the local phase change caused by the magnetic susceptibility change between tissues can be extracted.
The global phase change eliminator 333 of the present embodiment performs as the global phase change eliminating process, firstly applies a low-pass filter process to every two-dimensional image of the three-dimensional image (original image) being obtained, and calculates a low resolution image. Then, the original image is subjected to complex division by the low resolution image. Accordingly, the global phase change included in the low resolution image is eliminated from the original image.
There are various methods for eliminating the global phase change. There is one method, by way of example, that extracts the global phase change by performing a fitting on a three-dimensional image according to a low-order polynomial expression, and subtracts the global phase change from the original image. In eliminating the global phase change, this kind of different methods may be applicable.
It should be noted that the phase unwrapping process and the global phase change eliminating process of the phase image processing are just examples, and they are not limited to those examples. Any one of those two processes may be omitted. Those processes may be performed in any order. Neither of those two processes may be performed.
The magnetic field distribution converter 334 utilizes the formula 1, to convert the phase image P after the phase image processing is performed, into the magnetic field distribution δ.
There will be described one example of the magnetic field distribution calculating process flow, according to the functions of the magnetic field distribution calculator 330 of the present embodiment.
The phase image calculator 331 calculates the phase image P from the complex image I (step S1101).
Then, the magnetic field distribution calculator 330 performs various phase image processing on the phase image P, thereby obtaining a phase image after the phase processing. In the present embodiment, for example, the phase unwrapping processor 332 performs the phase unwrapping process on the phase image P being calculated (step S1102), and the global phase change eliminator 333 performs the global phase change eliminating process on the phase image P after the phase unwrapping process is performed (step S1103). As described above, the phase image processing may be performed in any order. Only one of the phase image processes may be performed, or neither of them may be performed. Furthermore, the present embodiment is not limited to those example as discussed above.
Lastly, the magnetic field distribution converter 334 converts the phase image after the phase image processing, into the magnetic field distribution (step S1104), and the process ends.
Next, there will be described the magnetic susceptibility distribution calculating process of the aforementioned step S1004, according to the magnetic susceptibility distribution calculator 340. The magnetic susceptibility distribution calculator 340 of the present embodiment utilizes the magnetic field distribution δ calculated by the magnetic field distribution calculator 330, so as to obtain the magnetic susceptibility distribution χ.
The magnetic susceptibility distribution calculator 340 performs a minimization process that minimizes an error function through an iterative operation, the error function being defined by a predetermined initial magnetic susceptibility distribution χ0 and the magnetic field distribution δ, thereby calculating the magnetic susceptibility distribution χ. In this situation, the updated magnetic susceptibility distribution χUPD calculated according to the error function, during the minimization process, is smoothed every time the updated magnetic susceptibility distribution is calculated. It should be noted that the updated magnetic susceptibility distribution χUPD is obtained by updating the initial magnetic susceptibility distribution χ0, on the basis of the error function.
In order to implement those as described above, the magnetic susceptibility distribution calculator 340 of the present embodiment is provided with an error function setter 341 configured to set the error function, an updated magnetic susceptibility distribution calculator 342 configured to calculate the updated magnetic susceptibility distribution χUPD from the magnetic field distribution δ and initial magnetic susceptibility distribution χ0, by using the error function being set, a smoothing processor 343 configured to smooth the updated magnetic susceptibility distribution χUPD, and to generate a combined magnetic susceptibility distribution χADD, and a determiner 346 configured to determine the end of the iterative operation, along with updating the initial magnetic susceptibility distribution χ0 to the combined magnetic susceptibility distribution χADD being calculated.
The smoothing processor 343 is provided with, a filter processor 344 configured to generate the magnetic susceptibility distribution after the filter process is performed (processed magnetic susceptibility distribution) χSMT, by applying a predetermined smoothing filter to the updated magnetic susceptibility distribution χUPD, and a combiner 345 configured to combine the updated magnetic susceptibility distribution χUPD and the processed magnetic susceptibility distribution χSMT, so as to generate the combined magnetic susceptibility distribution χADD.
As shown in the formula 2, the magnetic field distribution δ is obtained by applying the spatial convolution integral to the magnetic susceptibility distribution χ. If the formula 2 is expressed using vectors and matrix, intended for all the pixels within the phase image P, the magnetic field distribution is expressed as the following formula 6.
δ=FDχ (6)
where δ is the column vector of the magnetic field distribution having the magnitude of the entire pixel count N, and χ is the column vector of the magnetic susceptibility distribution to be calculated, and FD is the matrix having the magnitude of N×N, corresponding to the convolution operation applied to the magnetic susceptibility distribution χ.
The updated magnetic susceptibility distribution calculator 342 of the present embodiment obtains the magnetic susceptibility distribution χ from the formula 6, according to the least-square method. In this situation, the error function is introduced and minimized, whereby a solution is obtained.
The error function setter 341 of the present embodiment sets the error function that is used when the updated magnetic susceptibility distribution calculator 342 calculates the updated magnetic susceptibility distribution. The error function e(χ0) being provided may be expressed by the following formula 7, for instance.
[Formula 7]
e(χ0)=∥W(δ−FDχ0)∥22 (7)
where W is the matrix having the weight vector (hereinafter, referred to as “weight W”) as the diagonal component calculated on the basis of the phase value component of the complex image I, and FD is the operator representing the convolution integral, and ∥a∥22 represents squares of L2 norm (Euclidean norm) of optional vector a.
In other words, the error function is provided by assigning the weight W to the function defined by the magnetic field distribution δ and a result of the convolution operation performed on the initial magnetic susceptibility distribution χ0. This weight W is calculated on the basis of the phase value component (phase image P) of the complex image as discussed above.
Specifically, the error function setter 341 calculates phase fluctuation of each pixel from the phase image P, and further calculates the weight W according to the magnitude of thus obtained phase fluctuation. By way of example, the calculation is made in such a manner that the larger is the phase fluctuation, the smaller becomes the weight W.
A statistical indicator such as standard deviation may be employed to represent the phase fluctuation, for example, calculated by using plural pixel values (phase values) in the region around the pixel targeted for calculating the phase fluctuation.
It should be noted that the phase image P to be used may be the one obtained from the complex image I, or the one obtained after performing the aforementioned various phase image processing on the phase image P obtained from the complex image I.
As discussed above, the updated magnetic susceptibility distribution calculator 342 uses the error function e (χ0) to calculate the updated magnetic susceptibility distribution χUPD from the initial magnetic susceptibility distribution χ0. In this example here, the updated magnetic susceptibility distribution calculator 342 calculates the gradient ∇e(χ0) of the error function e(χ0), and calculates the updated magnetic susceptibility distribution χUPD according to the following formula 8, by using thus calculated gradient ∇e(χ0):
χUPD=χ0+∇e(χ0) (8)
In the present embodiment, the initial magnetic susceptibility distribution χ0 which is initially provided may be optional. By way of example, a distribution with all the values set to zero may be the initial magnetic susceptibility distribution χ0.
The smoothing processor 343 performs smoothing by applying the smoothing filter to the updated magnetic susceptibility distribution χUPD, thereby obtaining the processed magnetic susceptibility distribution χSMT. In this process, the filter processor 344 employs an edge preserving smoothing filter as the smoothing filter, and performs a linear smoothing process adaptively to local magnetic susceptibility fluctuation of each pixel value in the updated magnetic susceptibility distribution χUPD.
Specifically, assuming an optional pixel position in the updated magnetic susceptibility distribution χUPD is r, an average value in the local region n×n (n is an integer at least one) around the position is μ(r), and dispersion is σ(r)2, the pixel value χSMT(r) of the processed magnetic susceptibility distribution χSMT being the magnetic susceptibility distribution after applying the edge preserving smoothing filter at the position r, is expressed by the following formula 9:
where ν2 is the parameter corresponding to noise dispersion mixed into the magnetic susceptibility distribution χ. In the present embodiment, there is employed as the parameter, an average value of the local dispersion σ(r)2 at each pixel in the updated magnetic susceptibility distribution χUPD.
According to the edge preserving smoothing filter which is applied in the present embodiment, as indicated by the formula 9, smoothing acts strongly on the region where local dispersion is small, yielding a value close to the local average value, and thus a quantitative value of the magnetic susceptibility can be maintained. On the other hand, the region where the local dispersion is large, the smoothing acts weakly, preserving an edge portion that includes a fine structure. Therefore, by applying this edge preserving smoothing filter to the updated magnetic susceptibility distribution χUPD, noise in the magnetic susceptibility distribution is reduced, and an edge between fine veins and between tissues with large magnetic susceptibility fluctuation is automatically extracted and preserved.
In the present embodiment, the aforementioned edge preserving smoothing filter is employed, but a filter used in the smoothing process is not limited to this example. By way of example, a publicly-known edge preserving smoothing filter, such as a median filter and a bilateral filter, may be employed.
The combiner 345 merges the updated magnetic susceptibility distribution χUPD with the processed magnetic susceptibility distribution χSMT, thereby calculating the combined magnetic susceptibility distribution χADD. Merging is performed so that increase of noise is prevented along with maintaining the fine structure of the magnetic susceptibility distribution. Byway of example, merging is performed in such a manner that the Fourier component of the processed magnetic susceptibility distribution χSMT to which the edge preserve smoothing has been applied, is inserted into the vicinity of the magic angle in the k-space that is susceptible to noise.
The combiner 345 converts the updated magnetic susceptibility distribution χUPD and the processed magnetic susceptibility distribution χSMT into the k-space data, the k-space data items obtained after the conversion are respectively weighted, and are added together, thereby obtaining the combined magnetic susceptibility distribution χADD.
Specifically, the combiner 345 firstly converts the updated magnetic susceptibility distribution χUPD and the processed magnetic susceptibility distribution χSMT, respectively, into k-space data, according to Fourier transform. Then, each k-space data item (Fourier component) after conversion is weighted so as to achieve the objective as described above, and thus weighted k-space data items are added together. After the addition is performed, the combined magnetic susceptibility distribution χADD is obtained according to inverse Fourier transform.
Weighting in the k-space may be performed by creating a weighted image G and multiplying the k-space data by the weighted image G. At this time, in the weighted image GUPD by which the k-space data obtained from the updated magnetic susceptibility distribution χUPD is multiplied, for example, the weight assigned to the data in the vicinity of the magic angle region may be set smaller than the weight assigned to another region. On the other hand, in the weighted image GSMT by which k-space data obtained from the processed magnetic susceptibility distribution χSMT is multiplied, the weight assigned to the data in the vicinity of the magic angle region is set larger than the weight assigned to the data in another region.
This is because the updated magnetic susceptibility distribution χUPD has relatively large amount of noise, and it is necessary to reduce contribution of data in the vicinity of the magic angle region, so as to prevent increase of noise. On the other hand, noise has been reduced in the processed magnetic susceptibility distribution χSMT, and therefore, data contribution in the vicinity of the magic angle region is increased so as to maintain the fine structure.
As the weighted images GUPD(k) and GSMT(k) that satisfy those conditions above, weighted images defined by the following formula 10 may be employed, for instance, by using the Fourier component D(k) of the point-dipole magnetic field and an adjustment parameter b for adjusting the magnitude of the magic angle region. It should be noted that k indicates data position in the k-space.
In the aforementioned formula 10, a region where the Fourier component value is smaller than the predetermined value (b) is assumed as the magic angle region in the updated magnetic susceptibility distribution χUPD.
In the case where the adjustment parameter b is 0, GUPD(k) becomes 1 and GSMT(k) becomes 0. Accordingly, this is a simple least-square method which minimizes the error function of the formula 6 without performing the edge preserving smoothing process.
The determiner 346 determines the end of the iterative operation in the minimizing process for minimizing the error function e(χ0) through the iterative operation. When it is determined that the process continues, that is, iteration of the operation is determined, the latest combined magnetic susceptibility distribution χADD at that timing is set as the initial magnetic susceptibility distribution χ0, thereby updating the initial magnetic susceptibility distribution χ0. When the end of the operation is determined, the latest combined magnetic susceptibility distribution χADD at that timing is outputted as the magnetic susceptibility distribution.
In the present embodiment, when the updated magnetic susceptibility distribution χUPD is calculated, the determiner 346 uses this updated magnetic susceptibility distribution χUPD and the initial magnetic susceptibility distribution χ0 at that point of time, to determine if it is the end of the operation. The determiner 346 utilizes an evaluation function f(χUPD) expressed by the following formula 11, for instance, and when its value is smaller than a predetermined threshold ε, the end of the operation is determined.
Next, a flow of the magnetic susceptibility distribution calculating process according to the magnetic susceptibility distribution calculator 340 of the present embodiment will be described.
In the present embodiment, the minimization process for minimizing the error function defined by the predetermined initial magnetic susceptibility distribution χ0 and the magnetic field distribution δ, is performed through the iterative operation, thereby calculating the magnetic susceptibility distribution χ. At this time, during the minimization process, the smoothing process is performed on the updated magnetic susceptibility distribution calculated from the error function, whereby the combined magnetic susceptibility distribution is obtained, updating the initial magnetic susceptibility distribution with the combined magnetic susceptibility distribution, and then the process is repeated. Specifically, the following process is performed.
The magnetic susceptibility distribution calculator 340 provides the initial magnetic susceptibility distribution χ0 410 (step S1201). Then, the error function setter 341 utilizes the initial magnetic susceptibility distribution χ0 410 and the magnetic field distribution δ 400, so as to provide the error function e(χ0) (step S1202).
Then, the updated magnetic susceptibility distribution calculator 342 calculates the gradient ∇e(χ0) of the error function (step S1203) and according to the aforementioned formula 8, the initial magnetic susceptibility distribution χ0 is updated, whereby the updated magnetic susceptibility distribution χUPD 420 is calculated (step S1204).
The determiner 346 utilizes the evaluation function f(χUPD) according to the aforementioned formula 11, and determines if the operation ends (step S1205).
When it is determined that the operation continues in step S1205, the filter processor 344 applies the predetermined smoothing filter to the updated magnetic susceptibility distribution χUPD 420 (Filtering), thereby calculating the processed magnetic susceptibility distribution χSMT 430 (step S1206).
The combiner 345 merges the updated magnetic susceptibility distribution χUPD 420 with the processed magnetic susceptibility distribution χSMT 430, and the combined magnetic susceptibility distribution χADD 440 is obtained (step S1207). Merging is performed by multiplying k-space data items 421 and 431, obtained through Fourier transform of the updated magnetic susceptibility distribution χUPD 420 and the processed magnetic susceptibility distribution χSMT 430, respectively, by the weighted images 422 and 432 according to the aforementioned formula 10, adding the products together, followed by inverse Fourier transform (IFT) on the result of the weighted addition.
The determiner 346 sets the combined magnetic susceptibility distribution χADD 440 as the initial magnetic susceptibility distribution χ0 410 (step S1208), returning to step S1202, and the process is repeated.
When the end of the operation is determined in step S1205, the determiner 346 assumes the latest combined magnetic susceptibility distribution χADD obtained at that point of time, as the magnetic susceptibility distribution χ outputted by the magnetic susceptibility distribution calculator 340 (step S1209), and finishes the process.
If the combined magnetic susceptibility distribution χADD is not obtained in step S1205, i.e., when the end of the operation is determined in the first operation, the combined magnetic susceptibility distribution χADD is calculated from the updated magnetic susceptibility distribution χUPD at that point of time according to the aforementioned method, and this calculated combined magnetic susceptibility distribution χADD is assumed as the magnetic susceptibility distribution χ to be outputted.
The determiner 346 may utilize the combined magnetic susceptibility distribution χADD to determine if it is the end of the operation. In other words, it is also possible that when the combined magnetic susceptibility distribution χADD is calculated, the combined magnetic susceptibility distribution χADD and the initial magnetic susceptibility distribution χ0 at this point of time are used to determine the end of the operation.
In this case, as the evaluation function for determining the end of the operation, f(χADD) as the following is employed.
The magnetic susceptibility distribution calculator 340 provides the initial magnetic susceptibility distribution χ0 410 (step S1301). Then, the error function setter 341 utilizes the initial magnetic susceptibility distribution χ0 410 and the magnetic field distribution δ 400, so as to provide the error function e(χ0) (step S1302).
Next, the updated magnetic susceptibility distribution calculator 342 calculates the gradient ∇e(χ0) of the error function (step S1303), and calculates the updated magnetic susceptibility distribution χUPD 420 according to the formula 8 (step S1304).
The filter processor 344 applies the predetermined smoothing filter to the updated magnetic susceptibility distribution χUPD 420 (Filtering), thereby calculating the processed magnetic susceptibility distribution χSMT 430 (step S1305).
The combiner 345 merges the updated magnetic susceptibility distribution χUPD 420 with the processed magnetic susceptibility distribution χSMT 430, whereby the combined magnetic susceptibility distribution χADD 440 is obtained (step S1306). Merging is performed by multiplying the k-space data items 421 and 431, obtained through Fourier transform of the updated magnetic susceptibility distribution χUPD 420 and the processed magnetic susceptibility distribution χSMT 430, respectively, by the weighted images 422 and 432 according to the aforementioned formula 10, adding the products together, followed by inverse Fourier transform (IFT) on the result of the weighted addition.
Next, the determiner 346 utilizes the evaluation function f(χADD) of the formula 12, so as to determine the end of the operation (step S1307).
When it is determined the operation continues in step S1307, the magnetic susceptibility distribution calculator 340 sets the combined magnetic susceptibility distribution χADD 440 as the initial magnetic susceptibility distribution χ0 410 (step S1308), returning to step S1302, and the process is repeated.
When the end of the operation is determined in step S1307, the latest combined magnetic susceptibility distribution χADD obtained at that point of time is assumed as the magnetic susceptibility distribution χ outputted by the magnetic susceptibility distribution calculator 340 (step S1309), and the process is finished.
If the determination is made according to the combined magnetic susceptibility distribution χADD, it is possible to obtain an image with a large effect of the filtering process, relative to the case where the determination is made according to the updated magnetic susceptibility distribution χUPD. On the other hand, if the determination is made according to the updated magnetic susceptibility distribution χUPD, it is determined by data closer to an actual magnetic susceptibility distribution, relative to the case where the determination is made by the combined magnetic susceptibility distribution χADD, and therefore this allows the determination on the basis of more appropriate scope of data.
As described above, the MRI apparatus of the present embodiment is provided with, the measurement controller 310 configured to measure the complex signals of the nuclear magnetic resonance signals generated from the subject 101, in response to irradiation with RF magnetic field pulses, the image reconstructor 320 configured to reconstruct from the complex signals, the complex image I including the pixel value which is a complex number, the magnetic field distribution calculator 330 configured to generate the phase image P from the complex image I and to calculate the magnetic field distribution δ from the phase image P, and the magnetic susceptibility distribution calculator 340 configured to calculate the magnetic susceptibility distribution χ from the magnetic field distribution δ, wherein, the magnetic susceptibility distribution calculator 340 performs the minimization process for minimizing through the iterative operation, the error function e(χ0) defined by the predetermined initial magnetic susceptibility distribution χ0 and the magnetic field distribution δ, thereby calculating the magnetic susceptibility distribution χ, and performs smoothing on the updated magnetic susceptibility distribution χUPD that is calculated based on the error function e(χ0) during the minimization process, every time the updated magnetic susceptibility distribution χUPD is calculated.
The magnetic susceptibility distribution calculator 340 is provided with the updated magnetic susceptibility distribution calculator 342 configured to use the error function e(χ0) so as to calculate the updated magnetic susceptibility distribution χUPD from the initial magnetic susceptibility distribution χ0, the smoothing processor 343 configured to smooth the updated magnetic susceptibility distribution χUPD, so as to generate the combined magnetic susceptibility distribution χADD, and the determiner 346 configured to determine the end of the iterative operation, and to set the combined magnetic susceptibility distribution χADD as the initial magnetic susceptibility distribution χ0. The smoothing processor 343 is provided with the filter processor 344 configured to apply the predetermined smoothing filter to the updated magnetic susceptibility distribution χUPD, so as to obtain the processed magnetic susceptibility distribution χSMT, and the combiner 345 configured to merge the updated magnetic susceptibility distribution χUPD with the processed magnetic susceptibility distribution χSMT, and to generate the combined magnetic susceptibility distribution χADD. The determiner 346 may set the combined magnetic susceptibility distribution χADD as the magnetic susceptibility distribution χ, which is obtained when the end of the operation is determined.
The error function e(χ0) is set by assigning the weight W to the function defined by the magnetic field distribution δ and the result of the convolution operation performed on the initial magnetic susceptibility distribution χ0. The weight W may be calculated based on the phase image P obtained from the complex image I.
As thus described, in the present embodiment, the iterative operation of the smoothing process under the constraint of the relational expression between the magnetic field and the magnetic susceptibility, enables calculation of the highly precise magnetic susceptibility distribution.
In other words, in the present embodiment, the updated magnetic susceptibility distribution χUPD is smoothed, which is obtained in the course of the minimization process for minimizing the error function, and this process is repeated. In the smoothing, the processed magnetic susceptibility distribution χSMT smoothed by the smoothing filter and the updated magnetic susceptibility distribution χUPDcalculated by using the gradient of the error function, are weighted in k-space and merged together, so that noise is reduced and the fine structure is maintained.
As described above, in the present embodiment, data of noise part is averaged by filtering according to the smoothing filter, allowing the noise in the calculated magnetic susceptibility distribution χ to be reduced. Further, according to the merging process, the fine structure can be maintained in the magnetic susceptibility distribution χ.
Unlike the conventional methods, a regularization term is not used in the present embodiment, and thus it is possible to obtain the magnetic susceptibility distribution χ where noise is reduced along with maintaining the fine structure, through the simple process. In other words, according to the present embodiment, enormous calculations for setting adequate regularization parameters are not required.
According to the present embodiment, values in the vicinity of the magic angle region are not truncated in solving the ill-conditioned problem. Therefore, artifacts caused by this truncation may be reduced, thereby enhancing accuracy of thus calculated magnetic susceptibility distribution.
It should be further noted that the processed magnetic susceptibility distribution χSMT used for reducing noise and maintaining the fine structure in the updated magnetic susceptibility distribution χUPD, and the combined magnetic susceptibility distribution χADD for updating the initial magnetic susceptibility distribution χ0 in the iterative operation, are both calculated from the updated magnetic susceptibility distribution χUPD. In other words, these are data items that depend on phase information.
In the present embodiment, such data depending on the phase information is used to update the magnetic susceptibility distribution calculated through the iterative operation of the minimization process, and thus this yields an advantage that inconsistencies may be reduced between the measured phase information and the magnetic susceptibility distribution being calculated.
According to the present embodiment, the weight W in the error function is calculated from the phase image. The phase image is higher in precision in reflecting the phase fluctuation, relative to an absolute value image, and thus, it is possible to obtain the weight W on which this phase fluctuation is reflected precisely. Since the magnetic susceptibility distribution χ is calculated according to the error function using this weight, the magnetic susceptibility distribution can be calculated with a higher degree of precision.
As described above, according to the present embodiment, smoothing process is repeated through the simple process, preserving edge information such as the fine structure along with maintaining the relationship between the magnetic field and the magnetic susceptibility, thereby reducing artifacts and noise, and enhancing accuracy and precision of the magnetic susceptibility distribution being calculated.
The weighted image G which is used when the combiner 345 merges the updated magnetic susceptibility distribution χUPD with the processed magnetic susceptibility distribution χSMT, is not limited to the weighted images as described above. By way of example, GUPD(k) and GSMT(k) may be constant values, respectively. That is, the entire data within each k-space may be weighted homogeneously.
Assignment of the unique weight may allow easy control of SNR (signal-to-noise) of the image that is ultimately obtained.
In the aforementioned embodiment, the weight W in the error function e(χ0) that is provided by the error function setter 341 is calculated from the phase value component of the complex image I, but this is not the only example. The weight W may be calculated from an absolute value component of the complex image I. Calculation from the absolute value component may allow the weight to be obtained through an easy and simple process.
| Number | Date | Country | Kind |
|---|---|---|---|
| 2014-228843 | Nov 2014 | JP | national |
| Filing Document | Filing Date | Country | Kind |
|---|---|---|---|
| PCT/JP2015/079477 | 10/19/2015 | WO | 00 |
