WO2009019722A1

WO2009019722A1 - Generation method of parametric functional maps from magnetic resonance images

Info

Publication number: WO2009019722A1
Application number: PCT/IT2007/000561
Authority: WO
Inventors: Marco Becchio; Giovanna Rizzo; Alberto Bert; Christian Bracco
Original assignee: Im3D S.P.A.
Priority date: 2007-08-03
Filing date: 2007-08-03
Publication date: 2009-02-12
Also published as: WO2009019722A8

Abstract

Method comprising the operations of acquiring a first plurality of magnetic resonance images, each comprising a plurality of voxels, and calculating for each voxel a first parameter (Tlpre) representing the time required for the longitudinal component of a magnetization vector of the tissue to recover a predetermined value; acquiring in a first instant of acquisition (to) a second plurality of magnetic resonance images, each comprising a plurality of voxels, and measuring a first signal intensity value (S0) for each voxel in the first instant of acquisition (to). The method also comprises the operations of acquiring, in second instants of acquisition (t1; t2, t3, t4, t5, t6), a third plurality of magnetic resonance images after the injection of a contrast medium, each image comprising a plurality of voxels, and calculating for each voxel second parameters (Tlpost) representing the time required for the longitudinal component of a magnetization vector present in the tissue to recover a predetermined value after the administration of the contrast medium. The method is characterized in that it also comprises the operations of calculating for each voxel of the said third plurality of images a plurality of second intensity values (ΔRn) as a function of the said first and second parameters, and of calculating, for each voxel, a concentration curve of the contrast medium as a function of the second instants of acquisition (t1, t2, t3, t4, t5, t6) on the basis of the second intensity values (ΔRn). The method also comprises the operations of selecting, as a function of the first intensity values (So), a first subset of voxels of the third plurality of images, and calculating, for each voxel, a third intensity value found by integrating the concentration curve; selecting, as a function of the third intensity value, a second subset of voxels of the first subset of voxels and estimating, for each voxel of the second subset of voxels, third parameters (Ktrans, Kep, vp) representing the kinetics of the contrast medium in the voxel, these parameters being capable of being represented in parameter maps relating to the tissue.

Description

Generation method of parametric functional maps from magnetic resonance images

The present invention relates to a generation method of parametric functional maps from magnetic resonance images, and more specifically to a method of automatically generating colour maps from a dynamic three-dimensional report of the breast.

Tumour masses with dimensions of more than 1-2 mm, such as those typically present in a breast, require a network of blood vessels to feed them. The vessels which perfuse them differ from those of healthy tissue in their size, shape and permeability. The capillary ball in a lesion is disproportionate to the size of the lesion itself and is excessively permeable, in as much as the endothelium of the vessels of the capillary ball has a different structure allowing the molecules of the contrast medium to flow out more rapidly than in normal vessels.

By injecting a contrast medium into the blood and acquiring image volume packets for the same anatomical region before, during, and after the injection of the contrast medium, it is possible to display the tumour masses which show altered perfusion due to their characteristics.

The method used to acquire images of the blood vessels of lesions is known as DCE-MRI (Dynamic Contrast Enhanced Magnetic Resonance Imaging).

In the case of DCE-MRI of the breast, the image volumes into which a magnetic resonance report is typically divided are acquired at a rate (temporal frequency) which is commonly kept low in the interests of an high spatial resolution. In particular, six or seven image volumes are acquired in a period of eight to ten minutes of the dynamic sequence.

The acquisition of an image volume can be two- or three-dimensional. In two-dimensional acquisition, successive anatomical sections are acquired and each section is represented as a two-dimensional image in which the single element of the image is called voxel. By successively acquiring a number of anatomically contiguous sections it is possible to generate an image volume formed by a plurality of stacked two-dimensional images; the x and y coordinates of the volume denote the position of the voxel considered in a two- dimensional image, the z coordinate corresponds to the position of the two-dimensional image in the image volume. In three-dimensional acquisition, the sections making up the image volume are acquired by an intrinsically three-dimensional sampling scheme.

In magnetic resonance images, the signal intensity of a single voxel depends on three physical parameters which characterize the tissue displayed by the voxel. The three parameters are:

- The longitudinal relaxation time (Tl): the time required for the longitudinal component of the magnetization vector, intrinsically present in the tissue represented by the voxel, to recover 63% of its value at thermodynamic equilibrium. The term "thermodynamic equilibrium" means the situation in which the magnetization vector is aligned and in agreement with the static magnetic field vector generated by the tomography machine used for the magnetic resonance. The direction of the magnetization vector of the tissue, which at thermodynamic equilibrium is in agreement with that of the static magnetic field, is modified by the interaction with the radio frequencies transmitted before the magnetic resonance acquisition, these frequencies generating a transverse component of this magnetization.

- The transverse relaxation time (T2): the time required for the transverse component of the magnetization vector to fall to 37% of its initial value. The transverse component is generated by the radio frequencies transmitted before the acquisition, these frequencies causing a forced phase alignment of all the magnetic spins present in the voxel.

- The proton density (D): the modulus of the tissue magnetization vector of the voxel at thermodynamic equilibrium.

The acquisition of magnetic resonance images allows to produce images in which the signal intensity of each voxel is weighted with respect to one or more of the preceding parameters, using specific acquisition sequences, particularly those of the "gradient echo" type, such as FAST SPGR GRE 3D.

An acquisition sequence consists in sending one or more radio frequency pulses spaced apart in time; each sequence is characterized by a repetition time of the radio frequency pulse (TR), which is the time that elaps between two consecutive radio frequency pulses, and an echo time (TE), which is the time elapsing between the sending of one radio frequency pulse and the maximum rephasing of the signal. The pulse repetition time and the echo time are selected by the operator, and the weight of the parameters Tl, T2 and D in the signal intensity of each voxel depends on them. In each sequence, the dependence between the parameters of the sequence (TR and TE) and the physical parameters (Tl, T2 and D) is expressed mathematically by an equation. For example, in the "gradient echo" acquisition sequence, the signal intensity (S) in the voxel (x,y,z) is expressed as: \ m (1)

In this case, the signal is mainly weighted in D and Tl; the weighting in T2 is negligible, because TE is typically set to be very short with respect to T2.

In the above formula, G is a constant which depends on the acquisition geometry and on signal amplification factors intrinsic to the scanner used, and α is the "flip angle", defined as the angle between the direction of the magnetization vector at the end of the sending of the radio frequency pulse and the static magnetic field vector of the scanner.

During an acquisition sequence, a plurality of image volumes is acquired, these volumes being sets of data, images or three-dimensional matrices, in which each element (voxel) is associated with three spatial coordinates: x and y, which represent the position in the image, and z, which indicates the number of the image.

A dynamic report consists of seven image volumes, each composed of 66 images. The first image volume is acquired at the moment of the injection of the contrast medium; the instant of time of the start of acquisition is called to.

The next six image volumes are acquired after the injection of the contrast medium. The instants of the start of acquisition of the individual volumes are called t^t^t^U,^ and ^.

The dynamic report, acquired by the FAST SPGR GRE 3D sequence, allows to do an acquisition with fat suppression, that is to say one in which the voxels corresponding to the fat have a low signal.

This feature is very important, since the breast is composed of a percentage of fat and a percentage of gland which is a dense tissue. The dense tissue tends to decrease and the fat component increases with age in a woman. It is therefore convenient to acquire images in which the signal from the gland is higher than that from the fat. Therefore, in the resulting images, the gland component is visible with a good contrast with respect to the surrounding tissues (the fat is suppressed). The dynamic report above described, in other words the set of seven image volumes, is associated with two static reports weighted in Tl, acquired before the injection of the contrast medium. Each of these two reports generates a single image volume, in other words a single three-dimensional matrix.

These static reports, produced in two-dimensional mode in order to make their acquisition faster, allow to derive, for each voxel, the Tl of the tissue contained therein before the injection of contrast. The pre-contrast Tl values (Ti_pre) for each voxel are found by using the following equation (E.M. Haake, R.W. Brown, M.R. Thompson, R. Venkatesan, "Magnetic resonance Imaging, Physical Principles and Sequence Design, Wiley Liss", 1999):

1 — Llπ since ■ COSCC₂

'o2 -1 S₀₂SVOa₁ coscc,

71 pre TR S_alsina₂ TR S_alsina₂ where S_αi(x,y,z) and S_α2(x,y,z) are signal intensity values of the single voxel (x,y,z) measured in the image volumes acquired at different "flip angles" αj and α₂, particularly with an angle Cc₁ of 50° and an angle α₂ of 9°.

When the pre-contrast Tl has been found in this way, the post-contrast Tl (l/T_lpOst) is acquired for all the voxels in the six image volumes of the dynamic view acquired after the injection of the contrast medium, using equation (3) (E.M. Haake, R.W. Brown, M.R. Thompson, R. Venkatesan, "Magnetic resonance Imaging, Physical Principles and Sequence Design,

where A is the following quantity:

where S₀ and S_n are signal intensity values of the same voxel measured during the dynamic report with respect to the time to and t_n.

During the acquisition of the six image volumes of the dynamic report, the parameters TR, TE and the "flip angle" are not modified. On the other hand, the signal intensity S varies, because Tl changes in the different acquisition instants t_ls ..., t₆, owing to the variations of concentration of the contrast medium in the same voxel. In equations (3) and (4), the "flip angle" of the dynamic report is used, particularly an angle α of 15°.

At this point, traditionally, the image processing takes place, the images of the dynamic report obtained in pre-contrast and post-contrast modes being compared to identify potential tumour masses.

For this procedure, however, a very large number of comparisons must be made, since there are numerous images to be analysed. Moreover, these comparisons are not completely objective, since they depend on the experience of the operator executing them.

One object of the present invention is therefore to provide a generation method of parametric functional maps which enables a small number of images with improved display characteristics to be obtained from the images obtained in pre-contrast and post- contrast modes, thus facilitating the reading of the images and making their evaluation more objective, regardless of the operator's ability.

Further characteristics and advantages of the invention will be made clear by the following detailed description, provided purely by way of non-limiting example, with reference to the appended drawings, in which: - Figure 1 is an example of a concentration curve of a voxel.

The first operation of the generation method of parametric functional maps according to the invention is the application of a two-dimensional median filter with a 3x3 mask to all the image volumes of the dynamic report, in order to improve the signal to noise ratio. The filter is applied to each of the images making up a given image volume. This operation, which is known, is carried out as described in R.C. Gonzalez, and R.E. Woods, "Digital image processing", 2nd edition, Prentice-Hall Inc., New Jersey, 2002, Chapt. 3.6.2, pp. 123-124.

The six image volumes of the dynamic report acquired after the injection of the contrast medium are then aligned with the first image volume of the same dynamic report acquired before the injection of the contrast medium, in order to recover and compensate possible movements of the patient during the acquisition. The two pre-contrast static volumes are also aligned with the said first image volume of the dynamic report. These alignment operations are carried out using a mutual information maximization method, for example of the type described in D. Mattes, D. R. Haynor, H. Vesselle, T. K. Lewellen, and W. Eubank, "Non-rigid multimodality image registration" in Proc. SPIE Medical Imaging 2001: Image Processing, vol. 4322, 2001, pp. 1609-1620.

For each voxel, in each of the six image volumes of the dynamic report acquired after the injection of the contrast medium, a modified intensity value ΔR_n

- R_Pre is calculated, where Rpostn is the longitudinal relaxation constant measured at t_n, and R_pre is the longitudinal relaxation constant measured at to. R_p0St is the inverse of the longitudinal relaxation time after the injection of contrast, i.e. R_poSt = 1/T_lpost. R_pre is the inverse of the longitudinal relaxation time before the injection of contrast, i.e. R_pre = 1/T)_pre.

The calculation of these two quantities is carried out as described above.

The values of modified intensity ΔR_n found in this way are stored in n three-dimensional images (ΔR_n matrices), which describe, for each voxel, the variation of the concentration of the contrast medium with time. Since six post-contrast image volumes are acquired, there will be a total of six ΔR_n matrices of longitudinal relaxation difference.

When a voxel having the position (x,y,z) is selected in each of the six ΔR_n matrices, a concentration curve 1 (see Figure 1) for the contrast medium over time is obtained for that voxel, to which various pharmacokinetic models can be applied for the quantitative characterization of the tumour tissue in the selected voxel. This concentration curve 1 is obtained by plotting the modified intensity values ΔR_n as a function of the instants of acquisition ti, t^ for the voxel concerned.

A first predetermined threshold, with a value of 300 for example, is applied to the first image volume of the dynamic report, in other words that which has the acquisition time to. All the voxels having an intensity value So below said first threshold are identified as "not relevant" and are disregarded in all the subsequent stages of processing.

The value of said first threshold is slightly greater than the signal intensity value which is found in the voxels containing fat and can be used to eliminate all the voxels corresponding to background, lungs, air and fat (which are all clearly irrelevant regions), thus enabling the subsequent calculations to be limited to a smaller number of voxels and consequently reducing the processing time.

The ΔR_n matrices are used to generate a new AUGC matrix in which the semi-quantitative parameter "Area Under Gadolinium Curve" (AUGC) is calculated for the voxels marked as "relevant", as described in DJ. Collins, A.R. Padhani, "Dynamic Magnetic Resonance Imaging of Tumor Perfusion" in /EEE English Medical Biological Magazine, vol. 23, 2004, pp 65-83. The irrelevant voxels are set to 0.

In this new AUGC matrix, each individual voxel has an AUGC intensity value obtained found by processing the modified intensity values present for this single voxel in the n ΔR_n matrices. In the ΔR_n matrices, in each voxel, the modified intensity values describe the variation of the concentration of the contrast medium with time in that voxel. The AUGC intensity value is then obtained, for each voxel, as the area under the concentration curve.

The AUGC matrix is used to generate a binary mask, by setting all the voxels with an AUGC intensity value above 0 to 1 and setting all the other voxels to 0. The mask therefore contains only the voxels belonging to relevant regions and having a contrast capture of more than zero. The said mask is therefore a matrix of the same dimensions as those of the AUGC matrix, but with all the voxels equal to 1 or 0. This mask can be used for carrying out the subsequent operations, which are more expensive in computational terms, only on the relevant voxels (i.e. those marked 1 in the mask).

After a position voxel (x, y, z) has been selected in the ΔR_n, the value of the position voxel (x, y, z) in the mask is analysed. If this voxel has a value of 0, the application of the model described below is terminated; otherwise, if the value of the voxel is 1, the application continues. The procedure is then iterated on all the voxels of the ΔR_n matrices.

The aforesaid operation of applying the model is carried out using, for example, the two- compartment model described in P. S. Tofts, G. Brix, D.L. Buckley, J. L. Εvelhoch, Ε. Henderson, M. V. Knopp, H. B.W. Larsson, T. Lee, N. A. Mayr, G. J.M. Parker, R.Ε. Port, J. Taylor, R. M. Weisskoff, "Estimating Kinetic Parameters From Dynamic Contrast- Enhanced Tl -Weighted MRI of a Diffusable Tracer: Standardized Quantities and Symbols", in Journal of Magnetic Resonance Imaging, vol. 10, 1999, pp.223-232.

This model estimates, from the corresponding concentration curve, for each voxel marked 1, the value of three parameters, characterizing the kinetics of the contrast medium in this voxel, namely K_trans (constant of transfer from plasma to extravascular-extracellular space (EES)), k_ep (constant of transfer from EES to plasma), and v_p (plasma volume per unit of tissue volume). The parameters are estimated by using an optimizer which modifies the parameters of the model in such a way as to minimize a cost function, particularly the quadratic sum function of the residues. An example of an optimizer is the Variable Metric Minimizer. The optimizer iterates the operations of minimizing the cost function until a tolerance threshold on the variations of the minimum achieved, for example a threshold equal to le-7, is satisfied, in other words until the minimum of the cost function is found not to vary beyond said tolerance, or until a maximum number of iterations, for example 1000, is reached.

At this point, the optimizer restores the parameters of the model. Each parameter is stored in a corresponding parameter matrix or map, in other words in a map which describes the spatial distribution of the parameter.

During the estimation, a determination coefficient for each voxel is also calculated, as described in E. Furman-Haran, D. Grobgeld, H. Degani, "Dynamic Contrast-Enhanced Imaging and Analysis at High Spatial Resolution of MCF7 Human Breast Tumors", in Journal of Magnetic Resonance 1997, vol. 128, pp. 161-171. This determination coefficient is used as an index of the "closeness of fit".

A second threshold, for example 0.75, is applied to the parameter maps of spatial distribution of each parameter, and all the voxels having a determination coefficient below the said second threshold are set to 0.

At this point, a further median filter is applied to the parameter maps to eliminate isolated noise voxels. This operation is optional.

A third threshold is then applied to each of the three parameter maps to increase the contrast: the threshold relating to the K_trans parameter is fixed, for example, between 0 and 1 min^"1, and the voxels with value of less than 0 are set to 0 and those with value of more than 1 are set to 1 ; the threshold relating to the K^, parameter is also fixed, for example, between 0 and 1 min^"1, and the threshold relating to the v_p parameter is set for example, between 0 and 0.4.

Finally, each of the three parameter maps is encoded in a different colour channel, using a known RGB display system. The three colours are then combined to form a unique colour map in which the colour of a single voxel depends on the relative weight of the three parameters K_trans, K^₉ and v_p in the voxel concerned. Clearly, provided that the principle of the invention is retained, the forms of application and the details of construction can be varied widely from what has been described and illustrated purely by way of example and without restrictive intent, without thereby departing from the scope of protection of the invention as defined in the attached claims.

Claims

1. Generation method of parametric maps from magnetic resonance images of a tissue, comprising the operations of:

- acquiring in two-dimensional mode a first plurality of data representing magnetic resonance images, each image comprising a plurality of voxels;

- calculating for each voxel a first parameter (Tl_pre) representing the time required for the longitudinal component of a magnetization vector of the tissue to recover a predetermined value;

- acquiring in three-dimensional mode, in a first instant of acquisition (to), a second plurality of data representing magnetic resonance images, each image comprising a plurality of voxels;

- measuring a first signal intensity value (So) for each voxel in the said first instant of acquisition (to);

- acquiring in three-dimensional mode, in second instants of acquisition (U, t₂, t₃, t₄, t₅, U,) following the said first instant of acquisition (to), a third plurality of data representing magnetic resonance images after the administration of a contrast medium, each image comprising a plurality of voxels;

- calculating for each voxel second parameters (TI_p0St) representing the time required for the longitudinal component of a magnetization vector present in the tissue to recover a predetermined value after the administration of the contrast medium; the method being characterized in that it also comprises the operations of:

- calculating for each voxel of the said third plurality of data a plurality of second intensity values (ΔR_n) as a function of the said first and second parameters;

- calculating, for each voxel, a concentration curve of the contrast medium, as a function of the second instants of acquisition (ti, t₂, t₃, U, h, U) on the basis of the said second intensity values (ΔR_n);

- selecting, as a function of the said first intensity values (So), a first subset of voxels of the said third plurality of data;

- calculating, for each voxel of the said first subset of voxels, a third intensity value found by integrating the said concentration curve;

- selecting, as a function of the said third intensity value, a second subset of voxels of the said first subset of voxels;

- estimating, for each voxel of the said second subset of voxels, third parameters (Kt_rans, K_ep, V_p) representing the kinetics of the contrast medium in the voxel, these parameters being capable of being represented in parameter maps relating to the tissue.

2. Generation method of parametric maps according to Claim 1, additionally comprising the operation of applying a two-dimensional median filter with a 3x3 mask to the said second and third pluralities of data.

3. Generation method of parametric maps according to Claim 1 or 2, additionally comprising the operation of aligning the said third plurality of data with the said second plurality of data.

4. Generation method of parametric maps according to any one of the preceding claims, additionally comprising the operation of aligning the said first plurality of data with the said second plurality of data.

5. Generation method of parametric maps according to any one of the preceding claims, in which the said plurality of second intensity values (ΔR_n) are calculated by subtraction between the said second parameters and the said first parameter.

6. Generation method of parametric maps according to any one of the preceding claims, in which the operation of selecting a first subset of voxels comprises the operation of:

- checking for each voxel that the said first intensity value (S₀) is greater than a first predetermined threshold.

7. Generation method of parametric maps according to any one of the preceding claims, in which the operation of selecting a second subset of voxels comprises the step of:

- selecting the voxels having a third intensity value which is greater than a second predetermined threshold.

8. Generation method of parametric maps according to any one of the preceding claims, in which the operation of estimating third parameters comprises the steps of: a) minimizing a cost function; b) repeating the said minimization step until a minimum value of the said cost function is found not to vary beyond a predetermined tolerance; c) calculating a determination coefficient; d) selecting a third subset of voxels having a determination coefficient above a third predetermined threshold; e) comparing the third parameters of the voxels of the said third subset of voxels with corresponding final thresholds and setting each voxel to 0 or 1, relative to each third parameter, according to the outcome of these comparisons.

9. Generation method of parametric maps according to Claim 8, in which step b) is replaced by the following step f): f) repeating the said minimization step until a maximum number of iterations is reached.

10. Generation method of parametric maps according to Claim 8 or 9, in which the cost function is a quadratic sum function of the residues.

11. Generation method of parametric maps according to any one of the preceding claims, additionally comprising the operation of encoding the said parameter maps in a different colour channel in an RGB display system.

12. Generation method of parametric maps according to any one of the preceding claims, additionally comprising the operation of encoding the said parameter maps in different colour channels according to the RGB display system and generating a single RGB multiple-parameter map.