Method for reconstructing a planar tomogram from magnetic resonance signals in an inhomogeneous magnetic fields

[0002] 1. Field of the Invention

[0003] The present invention relates to a method for reconstructing a planar tomogram of an examination subject from magnetic resonance signals in an inhomogeneous magnetic field with a magnetic resonance device that has a known inhomogeneous main magnetic field and/or at least one known non-linear gradient magnetic field.

[0004] 2. Description of the Prior Art

[0005] Imaging by means of magnetic resonance techniques makes use of the frequency dependency of the magnetic resonance signal on the magnetic field strength for purposes of spatial resolution. Common methods for reconstructing tomograms presume a homogenous main magnetic field and strictly linear gradient magnetic fields. Inhomogeneities of the main magnetic field cause distortions or deformations in the frequency encoding direction. Given non-linearities of the gradient fields, the distortions are present not only in the tomographic image plane, but also perpendicular thereto, given slice excitations with a selection gradient. The distortions or deformations relate to the geometric position of the reconstructed spin density in the examination subject and likewise to the reconstructed image intensity.

[0006] Heretofore, distortions have been corrected in the image plane, exclusively. PCT Application WO 95/30908 teaches a method wherein a generalized Fresnel transformation is performed in the readout direction (GFT reconstruction). The GFT transformation takes into account a known location dependency of the main magnetic field in the readout direction.

[0007] The method described in German OS 195 40 837 uses two auxiliary data records which describe a shifting of the measured location relative to an actual location of a signal origin. From one of the auxiliary data records, a corrected auxiliary data record is created. In an image data record, a location correction then occurs in a first coordinate direction using the corrected auxiliary data record. A first intensity correction also occurs. Subsequently, a location correction in the second coordinate direction occurs, along with a second intensity correction. Alternatively, the location correction can occur by a Fresnel transformation of the raw data record including the corrected auxiliary data.

[0008] The article “Simulation of the influence of magnetic field inhomogeneity and distortion correction in MR imaging” (J. Weiss, L. Budinsky, Magnetic Resonance Imaging,8, 1990: 483-489) teaches the correction of image distortions by a postprocessing of a conventionally acquired image. The information about magnetic field inhomogeneities which is required for this is obtained from the phase of separately registered spin-echo images. Deformations in the slice direction, i.e. slice curvatures, cannot be corrected with this method.

[0009] The slice curvature problem with respect to the main magnetic field non-homogeneity conventionally has been avoided by a 3D imaging wherein a spatial resolution occurs in the slice direction by means of an additional phase encoding. The phase encoding is relatively insensitive to main magnetic field inhomogeneities relating to deformations. The longer measuring time compared to multi-slice pickups and the resulting higher susceptibility to movement artefacts are disadvantageous, however. There are consequently fundamental restrictions in the application of certain techniques which are based on a rapid imaging, such as contrast agent examination or dynamic studies. Deformations relating to non-linear gradient fields heretofore have remained uncorrected.

[0010] An object of the present invention is to provide a method for reconstructing a planar tomogram from magnetic resonance signals in inhomogeneous magnetic fields, with which method deformations in the direction of the main magnetic field also can be corrected.

[0011] The object is achieved in a method in accordance with the invention having the following steps. Original image elements are generated by means of multi-slice excitation, these original image elements being generated from magnetic resonance signals of original volume elements that are arranged in the examination subject in several curved layers that are disposed adjacent to one another succession. The spatial position of the original image elements is determined by means of the known magnetic fields. Image elements of a planar tomogram are generated, these image elements being representative of the volume elements of a planar slice that is situated within the curved slice, from the original image elements which are arranged in an environment of the respective image elements.

[0012] The term “original volume element” means the measured, i.e. distorted, volume element (voxel). “Volume element “means the real, undistorted voxel, which agrees with the measured one given ideal magnetic fields.

[0013] First, several neighboring slices which are deformed by non-ideal magnetic fields are measured. From this, at least one planar slice which is free of deformations is reconstructed. Altogether, the multi-slice excitation covers a 3D volume which, in contrast to the ideal case wherein the main magnetic field is homogeneous and the gradient magnetic fields are strictly linear, does not form a cuboid, but instead forms a region bounded by curvilinear surfaces. An advantage of the inventive method is a complete correction of deformations due to basic field inhomogeneities and gradient field non-linearities in all three spatial directions. This produces a functionality as in the ideal case. The position of the image elements indicates the actual anatomical position. Thus, for example, distance measurements in the entire image region are anatomically accurate. The slices of 3D-corrected measurements can then be graphically positioned onto 2D-corrected or 3D-corrected reference images with straight auxiliary lines, without further difficulty.

[0014] In one embodiment the image elements have respective intensity values, and the original image elements have respective original intensity values, and the intensity values are generated from surrounding original intensity values which are weighted with allocated weight factors. The preparation of the planar tomograms from the measured 3D volumes which are deformed due to non-ideal magnetic fields occurs using reconstruction methods which are derived from the known multiplanar reconstruction methods. The required precise spatial position of the original volume elements is obtained from coefficient tables of an allocated series expansion of the magnetic field, as are already used in the deformation correction in an image plane. Given main magnetic field inhomogeneities and gradient field non-linearities, the position of an image element is shifted. Together with the shift in the tomographical plane itself, which already can be calculated with the algorithms of the 2D deformation correction, the spatial position can be determined for each image element in a slice and for all neighboring slices using the coefficient tables. The spatial position of the original image elements enter into the method for multiplanar reconstruction as input values, in order to reconstruct a planar slice therefrom.

[0015] In another embodiment the weighting factors are determined by corresponding volume contents of the original volume element and the volume element. Given sharp deviations of the magnetic fields from the ideal distribution, the different volume contents of the original volume elements and the volume elements are taken into consideration in the methods for multiplanar reconstruction.

[0016] In a further embodiment, the volume contents are taken into account in the weighting factors by determining the weighting factors as a ratio of the volume contents of the respective volume elements relative to the volume contents of respective original volume elements.

[0017] In another embodiment, the weight factors are respectively formed by the volume contents of an overlap volume of the volume element and the original volume element relative to the original volume element. All original image elements that surround the image element to be reconstructed thus can be taken into account in the reconstruction.

[0018] The volume ratio of the original volume elements and the volume elements can also be determined in another embodiment by the weighting factors, obtained directly from the gradient magnetic field and main magnetic field at the location of the original volume element. Since the original volume elements are generally not limited to a plane in a straight line, the determination of the volume contents can require an increased computing outlay. This problem is avoided if the correction of the deformed intensity is calculated by the analytical term into which the gradient and main magnetic fields are entered.

[0019] The computing time for reconstruction can be minimized when conventional methods for multiplanar reconstruction are utilized. This is possible according to a further embodiment by transforming the original image elements into a curvilinear coordinate system in which the original image elements are arranged on a cubic raster, and by also transforming the image elements of the planar slice into the curvilinear coordinate system, and by forming the weighting factors as the ratio of the volume contents of the volume element relative to the volume contents of the original volume element.

[0020] FIG. 1 shows the isocontour lines of an inhomogeneous total magnetic field as arises in a magnetic resonance tomography apparatus.

[0021] FIG. 2 is a basic illustration of non-planar slices which are obtained in the presence of an inhomogeneous total field as shown in FIG. 1.

[0022] FIG. 3 illustrates image elements of a planar slice and non-planar slices.

[0023] FIG. 4 is a schematic illustration for explaining a first embodiment of the inventive method for modified multiplanar reconstruction.

[0024] FIG. 5 is a schematic illustration describing the approximation of a curvilinearly bounded volume element based on geometric base bodies.

[0025] FIG. 6 is a schematic illustration of the non-planar and planar slices in a cubic grid.

[0026] FIG. 7 is a schematic illustration of the non-planar and planar slices in a curvilinear grid subsequent to a transformation in accordance with the inventive method.

[0027] FIG. 8 is a schematic illustration for explaining a second embodiment of the inventive method for modified multiplanar reconstruction.

[0028] FIG. 9 is an illustration for explaining the determination of weighting factors which are used in the second embodiment of the inventive reconstruction method.

[0029] FIG. 10 is a schematic illustration for explaining a third embodiment of the inventive method for modified multiplanar reconstruction.

[0030] FIG. 1 schematically depicts isocontour lines 2 of the magnitude of the field strength of an inhomogeneous total magnetic field in an xy plane perpendicular to the direction of a main magnetic field in an imaging volume of a diagnostic magnetic resonance device. For clarification, a complete xyz coordinate system 4 is illustrated. The inhomogeneous total field is composed of the inhomogeneous main magnetic field and on-linear gradient fields. A plus sign 6 references an increased field strength with respect to a central region, while a minus sign 8 references a lower field strength. Since the frequency of the magnetic resonance signals is set proportional to the magnetic field strength, a location shift in the imaging arises in conventional reconstruction methods, via the frequency encoding, due to the inhomogeneous total field. The frequency encoding occurs in a slice selection and in the readout of the magnetic resonance signals. In the direction of the phase encoding, deviations from the ideal main magnetic field are generally not harmful, since only phase differences are evaluated. Typically, the magnetic field is characterized by a series expansion using Legendre polynomials. The corresponding measured or calculated coefficients are inserted here into the appertaining polynomial expansion. The magnetic field thus is obtained dependent on the spatial coordinates. The inhomogeneous field portions can be determined according to the following method, for example.

[0031] The inhomogeneous field portions of basic and gradient fields in a magnetic resonance system are measured or calculated.

[0032] The total field is represented as a series development using Legendre polynomials and the coefficients are stored in a file.

[0033] For the measured slices, the local, inhomogeneous field values are determined, using the Legendre series expansion, from the position vectors of the original image elements of the measured slices and from the coefficient file.

[0034] FIG. 2 shows, in a sectional view, that curved slices 10 which are excited in an exemplary inhomogeneous total magnetic field (FIG. 1), are not planar, but curved. In addition, the slice thickness can vary. Three adjoining slices 10 are illustrated as examples in cross-section in FIG. 2. The excitation of the neighboring slices 10 occurs via multi-slice sequences, nearly all known sequences being in principle multi-slice-capable, with the marginal condition that in inhomogeneous fields, spin-echo sequences deliver better results than do gradient echo sequences. Multi-slice spin-echo sequences and multi-slice turbo-spin-echo sequences, as well as multi-slice HASTE sequences, are preferably used. Altogether, a 3D volume is covered via the multi-slice excitation which, contrary to the homogenous case wherein the main magnetic field is constant and the gradient magnetic fields are strictly linear, does not form a cuboid, but instead forms a region bounded by curved surfaces. Images are now reconstructed of planar slices 12 which are situated within the deformed slice 10. The planar slices 12 are also illustrated in cross-section in FIG. 2. In the reconstruction, deformations due to main field inhomogeneities and gradient field inhomogeneities in all three spatial directions are corrected.

[0035] FIG. 3 depicts a cross-section of a section of the curved slice 10, in which the planar slice 12 to be reconstructed is arranged. The curved slice 10 is formed by adjoining original volume elements 14 which are bounded by curved surfaces. By means of multi-slice excitation, original image elements for a tomogram representation are generated, which render the spin density of the original volume elements 14. The spatial positions of the original volume elements, and thus the positions of the original image elements in the tomogram, are obtained from the known magnetic field course. Besides the spatial position, the original image elements include an intensity value as a measure of the spin density in the form of a grayscale value. From the original image elements, image elements of a planar tomogram are constructed, which are representative of volume elements 16 of the planar slice 12 within the curved slice 10. for this purpose, image data are evaluated of the original volume elements 14, which are in an environment of the volume elements 16 to be represented. Four adjoining original volume elements 14 are shown in hatched fashion as examples in FIG. 3, from whose image information the image element that is representative of the volume element 16 of the planar slice 12 is reconstructed. The intensity value or grayscale value of the reconstructed image element is derived from weighted grayscale values of the original volume elements 14.

[0036] For the reconstruction, the methods detailed below are used, these being derived from the multiplanar reconstruction (MPR). These methods can be generally represented as follows: If i(p) is the intensity value of an original image element for a point p in space with the coordinates x, y, z, then the intensity value i′ of an image element at the point p′ with the coordinates x′,y′, z′, which does not correspond to the point p, is defined by a function f of the intensity values of the original image elements at n environmental points p₁to P_n:

[0037] i′(p′)=f(i(p₁), . . . ,i(p_n))

[0038] The function f is often a linear combination of the intensity values of the original image elements with suitable weight factors a_k:

[0039] i′(p′)=a₁i(p₁)+ . . . +a_ki(p_k)+ . . . +a_ni(p_n).

[0040] In the conventional methods of MRP, the factors a_kare derived from a distance value

[0041] d_k=|p′−p_k|.

[0042] Thus, in the case of the 2D reconstruction in a conventional correction method, which is referred to as a nearest neighbor method, only the intensity value of the original image element situated nearest the image element to be reconstructed is used. The intensity values of the other neighboring image elements are not taken into account in the reconstruction.

[0043] FIG. 4 details a nearest neighbor method which has been modified for 3D reconstruction. The volume elements 16 belonging to a planar slice 12 are situated within two curved slices 10. The center of gravity is taken as the representative p_kof the original volume elements 14, for example. Likewise, the center of gravity p′ is applied as the representative of the volume element 16. The original volume element 14 having the shortest distance d_min=min{d_i, . . . ,d_n} to the volume element to be reconstructed is then determined on the basis of the centers of gravity. The intensity value of the volume element 16 to be reconstructed is obtained by multiplying the intensity value of the nearest original volume element 14 by the ratio of the volume contents |V′| of the volume element 16 and the volume content |V_k| of the original volume element 14, as a weighting factor according to

[0044] i(p′)=i(p_min)|V′|/|V_k|.

[0045] Since, in contrast to the volume element 16, the original volume elements 14 are generally not limited to a plane in a straight line, the determination of the volume contents V_kcan be more complicated than in the case of simple geometric bodies. To save computing power, given low precision requirements, the conventional nearest neighbor method can be used, wherein corrections using weighting factors dependent on the volume contents of the participating volume elements 14, 16 do not occur.

[0046] If, however, higher precision requirements must be satisfied, then approximating methods can be employed only given volumes V_kthat are extremely difficult to calculate analytically. One of these is the approximation of the volume contents with simple geometric base bodies, wherein the elementary size is determined by the specified precision. The simplest case is a volume determination by an approximation of the volume element by an elementary cuboid 18; for clarity, FIG. 5 outlines this for the two-dimensional case. The 3D case is analogous.

[0047] Another way to correct the deformed intensity is the use of an analytical term into which the gradient and main magnetic fields are entered. An intensity correction factor I_c, that is obtained directly from the magnetic field quantities, corresponds to the volume ratio of the reconstructed volume element 16 to the original volume element 14, and is calculated according to the following relation:

[0048] I_c(x,y,z)=A·B·C,

[0049] wherein $A=\left|1+\frac{1}{G_x}·\frac{\delta }{\mathrm{\delta x}}\left(\mathrm{\Delta B}\left(x,y,z\right)+B_x\left(x,y,z\right)\right)\right|B=\left|1+\frac{1}{G_y}·\frac{\delta }{\mathrm{\delta y}}{B}_y\left(x,y,z\right)\right|C=\left|1+\frac{1}{G_z}·\frac{\delta }{\mathrm{\delta z}}{B}_z\left(x,y,z\right)\right|$

[0050] with the terms defined as follows:

[0051] x, y, z: the actual coordinates of the volume element, not the measured ones

[0052] ΔB: basic field inhomogeneity

[0053] B_x: gradient non-linearity in the frequency coding direction

[0054] B_y: gradient non-linearity in the phase coding direction

[0055] B_z: gradient non-linearity in the slice selection direction

[0056] G_x: field gradient in the frequency coding direction

[0057] G_y: field gradient in phase coding direction

[0058] G_z: field gradient in the slice selection direction

[0059] The coordinates x, y, z define the location of the representative of the corresponding volume element. The intensity at this location is then multiplied by the factor I_cin order to obtain a corrected intensity at the location of the reconstructed image element. As an approximation, the coordinates of a desired reference point can be employed instead of the coordinates of the nearest neighbor point.

[0060] Another MPR method of correction of the intensity values is known as the “trilinear algorithm”. The conventional method can only process original image elements which are situated at intersections of a cubic grid, however. To achieve this, a curvilinear coordinate system Q is defined, in which the original image elements p are situated on a cubic grid. The points p′ of the deformation-corrected image are also transformed according to Q by a transformation T_Q.

[0061] FIGS. 6 and 7 illustrate the transformation T_Q. FIG. 6 shows the original volume elements 14 and the reconstructed volume elements 16 in the Cartesian coordinate system. FIG. 7 likewise shows original volume elements 14′ and reconstructed volume elements 16′ in the curvilinear coordinate system Q. The local and generally non-linear transformation T_Qcan be defined as follows: ${\underset{\_}{p^{\prime }}}_i=T_Q\left({\underset{\_}{p^{\prime }}}_i\right)·{\underset{\_}{p}}_i=\left(\begin{array}{c}{p}_{\mathrm{ix}}+\left(\Delta B\left({\underset{\_}{p}}_i\right)+B_x\left({\underset{\_}{p}}_i\right)\right)/G_x\\ p_{\mathrm{iy}}+\left(B_x\left({\underset{\_}{p}}_i\right)\right)/G_c\\ p_{\mathrm{iz}}+\left(\Delta B\left({\underset{\_}{p}}_i\right)+B_x\left({\underset{\_}{p}}_i\right)\right)/G_x\end{array}\right)=\left(\begin{array}{c}{p}_{\mathrm{ix}}^{\prime }\\ p_{\mathrm{iy}}^{\prime }\\ p_{\mathrm{iz}}^{\prime }\end{array}\right),$

[0062] when there is phase coding in the y-direction. If the phase coding is in a different coordinate direction, the main field inhomogeneity ΔB(p_i) in this direction is not taken into consideration, because the main field inhomogeneity in the phase encoding direction basically does not deform the image elements. Gradient field non-linearities, however, do also affect the phase encoding direction.

[0063] The transformation T_Qdescribes the allocation of vectors p_iof the original vector field to the vectors p′_iof the target vector field. The locality and non-linearity of the transformation T_Qis contained in the inhomogeneous total field, which is additively composed of the inhomogeneous basic field ΔB(p_i) and the non-linearity of the gradient field B_x(p_i),B_y(p_i), B_z(p_i). The two inhomogeneous field portions are usually given in a series representation, e.g. using Legendre polynomials.

[0064] By the transformation T_Qof the image elements into the curvilinear coordinate system Q, the deformed measured slices 10 are transformed into planar slices 10′, as can be seen in FIG. 7. The planar slice 12 is transformed into a curved slice 12′ by the transformation T_Q.

[0065] Conventional MPR methods which require a cubical basic grid can thus be used in order to determine the image elements of arbitrary slices. FIGS. 8 and 9 depict the geometric relations in the conventional trilinear algorithm that can be employed to approximately define the image elements to be reconstructed, in a curvilinear coordinate system, provided the precision requirements are not very high.

[0066] As representatives of the transformed original volume elements 14′, points p₁to p₈represent the vertices of a grid element 20 of the cubic grid that exists after transformation, in which grid the image element p′ to be reconstructed is located. The points p₁to p₈are the nearest points to p′. In the trilinear algorithm, the intensities i(p₁) to i(p₈) define the intensity of the reconstructed image element p′ with corresponding weight factors. The equation

[0067] i′(p′)=a₁i(p₁)+ . . . +a₈i(p₈)

[0068] applies here, the weighting factors a_kbeing determined from the projections d_kx, d_ky, d_kzof the distance vector d_kon the grid lines g_x, g_y, g_z. FIG. 9 clarifies these relations. The projections in the direction of the edges of the transformed original volume elements are as follows:

[0069] d_kx=|(p′−p_k)e_x|

[0070] d_ky=|(p′−p_k)e_y|

[0071] d_kz=|(p′−p_k)e_z|

[0072] The terms e_x, e_yand e_zrepresent the corresponding mutually perpendicular unit vectors.

[0073] The weight factors can be determined from the projections according to the following relation:

[0074] a_k=(g_x−d_kx)(g_y−dky)(g_z−d_kz)

[0075] If the precision in the application of the conventional trilinear algorithm is not sufficient, a modified weight factor a_k^modis used, which accounts for the different volumes of the original image elements and of the reconstructed image elements. The weighting factor a_kis modified in this regard by multiplying it by the ratio of the volume contents of the reconstructed volume element and the volume contents of the original volume element, according to the relation:

[0076] a_k^mod=|V′|/|V_k|a_k.

[0077] Another variant which is based on the calculation of the planar slice in the curvilinear coordinate system Q is illustrated in FIG. 10. The original volume elements are still in the cubic grid; the reconstructed volume element is limited in a curved line. According to this variation, the weighting factors for the corresponding intensities are defined by the ratio of the volume contents of the overlap volume |V′109 V_k| relative to the original volume element |V_k|, according to:

[0078] a_k^mod=|V′∩V_k|/|V_k|.

[0079] The volumes |V_k| of the original volume elements 14′ are easier to determine. The reconstructed volume elements 16′ are generally bounded curvilinearly, so that approximating methods are applied here only given volume contents which cannot be determined analytically or which can be so determined only with difficulty. One possibility is to approximate the volume contents by small geometric base bodies, as explained above.

[0080] Finally, the intensity correction with the previously indicated factor I_c(x, y, z) can be performed. The intensities of the eight nearest original image elements are respectively multiplied by the corresponding intensity correction factor. As in the conventional trilinear algorithm, the intensity values of these points then enter into the intensity, which is to be calculated, of the reconstructed image weighted according to their distance.

[0081] Although modifications and changes may be suggested by those skilled in the art, it is the intention of the inventor to embody within the patent warranted hereon all changes and modifications as reasonably and properly come within the scope of his contribution to the art.