WO2008151380A1

WO2008151380A1 - Method of parallel magnetic resonance image processing

Info

Publication number: WO2008151380A1
Application number: PCT/AU2008/000852
Authority: WO
Inventors: Jingxin Zhang; Zhaolin Chen
Original assignee: Monash University
Current assignee: Monash University
Priority date: 2007-06-14
Filing date: 2008-06-13
Publication date: 2008-12-18
Anticipated expiration: 2009-12-14

Abstract

A method, an apparatus, and a computer-implemented system are provided for performing reconstruction of an image of a subject within a magnetic resonance imaging (MRI) machine having a plurality of receiver coils. In particular, parallel magnetic resonance imaging (pMRI) reconstruction is performed by receiving pMRI signal data comprising digitised samples acquired from the receiver coils of the MRI machine. A dynamic input-output system of a pMRI detection process of the MRI machine is constructed, wherein the model embodies a plurality of estimated spatial sensitivity functions corresponding with each of the plurality of receiver coils. An inverse model is determined, which corresponds with the system model, using an optimisation method which is adapted to minimise a measure of reconstruction error. More particularly, the measure of reconstruction error accounts for accuracy of image reconstruction by the inverse model in the absence of additive noise and estimation error, in combination with amplification by the inverse model of additive noise and estimation error. An image of the subject is reconstructed by applying the inverse model to the pMRI signal data. A reconstructed image output may then be provided to a user.

Description

METHOD OF PARALLEL MAGNETIC RESONANCE IMAGE PROCESSING FIELD OF THE INVENTION

The present invention relates generally to magnetic resonance imaging (MRI) and, more particularly, to improved methods of image reconstruction in parallel magnetic resonance imaging (pMRI) systems. BACKGROUND OF THE INVENTION

Magnetic resonance imaging (MRI) is a non-invasive method which can be used to acquire and record images of the inside of a subject, which finds application in medical imaging where it is commonly used to identify pathological or other physiological alterations in living tissues, for example of a human subject. In conventional MRI systems, an extremely strong, uniform, static magnetic field is applied in order to polarise nuclear spins within molecules of the imaging subject. Static magnetic fields typically in the range of 0.5 to 2.0 Tesla may be generated using permanent magnets or resistive electromagnets, however in modern MRI apparatus the static fields are most commonly generated using superconducting electromagnets.

Typical MRI apparatus also includes three orthogonally-oriented gradient magnets. The gradient magnets are conventionally resistive electromagnets which are of relatively low strength compared to the main static magnetic field. As the name implies, the purpose of the gradient magnets is to generate a spatial variation in the magnetic field, ideally being a linear variation, which is controllable along three orthogonal axes. The gradient magnets enable spatial selection of the imaging function within the subject volume.

The generation of nuclear magnetic resonance (NMR) signals for MRI data acquisition is accomplished by exciting the magnetic moments of the polarised nuclear spins with a uniform radio frequency (RF) magnetic field, typically in the form of an RF pulse. This RF field is centred on a selected frequency that corresponds with a resonance frequency of nuclei within the subject under the influence of a magnetic field with a corresponding magnitude (ie for which the Larmor frequency is resonant with the frequency of the RF field).

During excitation, the nuclear spin system absorbs magnetic energy, resulting in a precession of magnetic moments of the nuclei around the direction of the main magnetic field. When the excitation is removed, the precessing magnetic moments undergo free induction decay (FID), releasing their absorbed energy and returning to a static state. During FID, NMR signals are detected using one or more receiver RF coils.

While various techniques for MRI image acquisition may be employed, for present purposes it is sufficient briefly to describe just one particularly common and effective representative method, namely Fourier Transform (FT) tomographic imaging. The simplest FT imaging sequence commences by applying a slice selection gradient field, simultaneously with an RF excitation pulse. The slice selective gradient results in one particular plane, or slice, within the subject volume experiencing a magnetic field whereby the nuclei within the selected plane have a Larmor frequency corresponding with the centre frequency of the RF excitation pulse. Accordingly, only the nuclei within the subject volume located within the selected slice undergo excitation under the influence of the RF pulse. An orthogonal phase-encoding gradient field is then applied, along one axis of the selected slice. The magnetic moment of the excited nuclei within the slice are thus caused to precess at a frequency corresponding with the local magnetic field under the influence of the phase-encoding gradient, such that the frequency of oscillation ideally varies linearly along the corresponding axis of the selected slice. Due to the variation in frequency, when the phase-encoding gradient field is removed there has been established a corresponding phase variation across the selected slice.

A frequency-encoding gradient field is then applied along the other, orthogonal, axis of the selected slice. This again results in a variation in frequency of precession, in the perpendicular direction across the selected slice. During application of the frequency-encoding gradient pulse the decaying response field is detected using the receiver RF coil, or coils. The resulting acquired signal comprises the sum of the response of all excited nuclei, each of which emits a field having a phase and frequency corresponding to its particular location within the two-dimensional plane of the selected slice.

While an enormous amount of information in relation to the MRI response of the selected plane within the subject volume is contained within the signal acquired via the receiver coil, there is insufficient information contained within a single measurement to reconstruct the complete image of the selected slice. In accordance with the FT imaging technique, the additional information required may be obtained by repeating the measurement using different values of phase-encoding. The applied phase shift may be changed, for example, by changing the magnitude of the phase-encoding gradient field, or by changing the length of time for which the phase-encoding gradient field is applied. The end result is a plurality of acquired waveforms corresponding with the plurality of applied levels of phase-encoding. Generally, each waveform is sampled and digitised in order to produce a complete data set that may be represented as a two-dimensional array, or matrix, comprising a plurality of rows and columns wherein each row corresponds with a particular acquisition and corresponding phase-encoding, and each column corresponds with a particular sample time during the acquisition of the waveform received by the receiver coil.

The two-dimensional array, or matrix, of phase/frequency encoded data is commonly known as a "k-space" representation of the MRI image. The corresponding "image-space" matrix, ie the excited spin density function which is effectively the MRI image itself, may be obtained by computing the Discrete Fourier Transform (DFT) of the k-space data. That is, k-space and image-space form a (two-dimensional) Fourier transform pair. A particular problem with MRI medical imaging is that it is necessary to make a large number of individual measurements (eg by varying the phase-encoding gradient) in order to acquire a sufficiently high-resolution image of the selected slice of the imaging subject. During the measurements, the subject, ie potentially a living human patient, must remain extremely still. It is therefore recognised as being highly desirable to accelerate the overall imaging process.

Accordingly, sophisticated multi-element coils have recently been developed, which are capable of acquiring multiple channels of data in parallel. This "parallel imaging", or parallel MRI (pMRI), technique enables accelerated imaging to be performed. In general, a "pMRI signal" comprises a plurality of signal components, one for each of a corresponding plurality of receiver coils, wherein each component may be represented as an individual matrix of digitised samples (either in k-space or image-space) as described above. In order to accelerate acquisition, the pMRI signal may be under-sampled in the phase-encoding direction, as a result of which every pixel (Ze matrix element) in the image-space components represents data from multiple points in space. The image arrays are accordingly "compressed", and subject to "wrap-around" artefacts (or aliasing) in a corresponding dimension.

In principle at least, reconstruction of the complete image from the individual signal components is straightforward. The general approach in image-space may be understood with the assistance of a specific example, in which a subject plane S is imaged with four coils A, B, C and D, and pMRI data is acquired, under-sampling by a factor of four. Wrap-around artefacts (aliasing) will occur, such that each point in the image arrays PA, PB, PC and PD corresponding with the four receiver coils represents data from four points in S, one in each of four distinct sections of the plane S-i, S₂, S₃ and S₄. Furthermore, each receiver coil has a distinct spatial sensitivity function CA, CB, C_C and CD. Ideally, the relationship between the measured image arrays PA, PB, PC and P₀ and the actual content of the sections Si, S₂, S₃ and S₄ of the overall plane S, is given by the following simultaneous (matrix) equations:

P, ^=C ₄(S₁ +S₂ +S₃ +S₄) P₅ =C₅(S_{1 +}S_{2 +}S_{3 +}S₄) P₀ =C^S₁ +S₂ +S₃ +S₄) P₀ =C₅(S_{1 +}S_{2 +}S_{3 +}S₄)

Solving the above equations for S-i, S₂, S₃ and S₄, and piecing these sections together, thus (in principle) enables the complete image of S to be obtained.

A similar analysis to the foregoing is also possible utilising k-space representations, resulting (again in principle) in a matrix solution which may be transformed using the DFT in order to obtain the reconstructed image. In practice, however, reconstructing the complete image in accelerated pMRI is not so straightforward. In a real MRI machine, the spatial sensitivity functions are not exactly known, /e there will always be some discrepancy, or estimation error, ΔC, between the actual spatial sensitivity functions and an estimate, denoted as C, which is used in image reconstruction. Furthermore, in real MRI apparatus there will always be additive noise in the acquired data signals. Due to both of the foregoing factors, there will be propagation of errors in any attempt to solve the simultaneous equations required to reconstruct the acquired image.

In many cases, due to the abovementioned problems, it may be difficult to obtain a satisfactory solution, ie to reconstruct an image of acceptable quality.

Nonetheless, there exist a number of prior art methods for pMRI image reconstruction. These include SENSE, PILS, SMASH, AUTO SMASH, and GRAPPA, each of which is described briefly below.

The SENSE method operates in the image domain. It requires an initial estimation by pre-imaging, and subsequent inversion of a spatial sensitivity matrix C, via the pseudo inverse C⁺. In many cases, this problem is ill- conditioned, particularly when significant acceleration is employed, resulting in amplification of noise and estimation errors, such that image reconstruction may be unacceptably poor. The PILS technique is also an image domain method, and is a special case of SENSE. In PILS, the sensitivity functions C are assumed to be ideally localised. This avoids the need to compute a corresponding inverse matrix, but is, in practice, of limited application. For example, the assumption of ideal, localised sensitivity functions may not be valid for real MRI machines, and additionally it is not possible to take advantage of techniques, such as over-sampling, in order to improve the overall signal-to-noise ratio (SNR) of the final image.

The SMASH method is a k-space technique with similar limitations to the (image domain) SENSE method. The AUTO SMASH method utilises specific assumptions in order to estimate the receiver coil spatial sensitivity functions C during real imaging, and accordingly generally requires a period of unaccelerated operation during initial measurements for this purpose. The GRAPPA method is a generalisation of AUTO SMASH which requires similar assumptions to be made. None of the prior art techniques allow any control over the accuracy with which the receiver coil responses (ie sensitivity functions) are estimated and inverted, in relation to the extent to which additive noise and errors are propagated and amplified within the reconstructed image. It is also noteworthy that a trade-off may exist between degree of acceleration and quality of reconstructed image, but that no existing technique enables quantification or optimisation of this trade-off.

Accordingly, there is an ongoing need for the development of new and improved reconstruction methods for accelerated pMRI imaging. SUMMARY OF INVENTION

In view of the shortcomings of the prior art, the present invention provides, in one aspect, a parallel magnetic resonance imaging (pMRI) reconstruction method for reconstructing an image of a subject within a magnetic resonance imaging (MRI) machine having a plurality of receiver coils, wherein accelerated imaging is performed by reducing a number of distinct MRI measurements performed while simultaneously receiving measurement data from said plurality of receiver coils, the method comprising the steps of: receiving pMRI signal data comprising digitised samples acquired from the plurality of receiver coils of the MRI machine; constructing a dynamic input-output system model of a pMRI detection process of the MRI machine, said model at least embodying a plurality of estimated spatial sensitivity functions corresponding with each of said plurality of receiver coils; determining an inverse model corresponding with said system model in accordance with an optimisation method which is adapted to minimise a measure of reconstruction error, wherein said measure of reconstruction error accounts for accuracy of image reconstruction by the inverse model in the absence of additive noise and estimation error, in combination with amplification by the inverse model of additive noise and estimation error; reconstructing an image of the subject by applying the inverse model to the pMRI signal data; and providing a reconstructed image output.

By comparison with prior art methods, the present invention advantageously enables an overall improvement in the image reconstruction by optimising over both the accuracy with which an inverse function of the spatial sensitivity of the receiver coils is established, and the corruption of the final reconstructed image due to amplification of additive noise and estimation errors. In particular,, there are many cases of practical interest, particularly where substantial acceleration is employed, in which improvements in the estimation of an inverse function for the MRI detection system are correlated, in pMRI image reconstruction, with the existence of ill-conditioned matrices (or functions), the inversion of which results in propagation and amplification of estimation errors and additive noise components. That is, in many practical cases a method that seeks to perform improved image reconstruction by more accurately estimating, and inverting, the spatial sensitivity functions of the receiver coils, may actually produce inferior results due to the corresponding amplification of noise and estimation errors. The present invention, by contrast, seeks for the first time to account for both of these contributions to overall reconstruction error, and to minimise the error resulting from the combination of both effects, using suitable optimisation techniques.

In various embodiments of the invention, the spatial sensitivity functions may be embodied either via image-space representations, or via k-space representations. In one preferred embodiment, the spatial sensitivity functions are embodied in an image-space representation as an aliasing component (AC) matrix. In alternative embodiments, however, the spatial sensitivity functions may be embodied in a k-space representation as a polyphase matrix. In particularly preferred embodiments, the optimisation method includes minimising, according to specified criteria, a measure of performance that comprises: (a) a term corresponding with a perfect reconstruction (PR) condition; and (b) a term corresponding with a gain of the inverse model. Advantageously, this approach seeks to avoid producing a reconstructed image based solely upon accuracy of reconstruction, but without consideration of the effects on image quality of amplified noise and propagated errors.

In accordance with an exemplary implementation, the optimisation method includes using the HL norms of the reconstruction error system and the inverse system to measure, respectively, the deviation from PR and the gain of the inverse model, and performing an HL norm optimisation. Preferably, the HL norm optimisation utilises an adaptively estimated weighting factor. However, it is also possible to utilise a predetermined or prescribed weighting factor in HL norm optimisation. In another aspect, the present invention provides a computer-implemented system for performing a parallel magnetic resonance imaging (pMRI) reconstruction method for reconstructing an image of a subject within a magnetic resonance imaging (MRI) machine having a plurality of receiver coils, wherein accelerated imaging is performed by reducing a number of distinct MRI measurements performed while simultaneously receiving measurement data from said plurality of receiver coils, the system comprising: means for receiving pMRI signal data comprising digitised samples acquired from the plurality of receiver coils of the MRI machine; means for constructing a dynamic input-output system model of a pMRI detection process of the MRI machine, said model at least embodying a plurality of estimated spatial sensitivity functions corresponding with each of said plurality of receiver coils; means for determining an inverse model corresponding with said system model in accordance with an optimisation method which is adapted to minimise a measure of reconstruction error, wherein said measure of reconstruction error accounts for accuracy of image reconstruction by the inverse model in the absence of additive noise and estimation error, in combination with amplification by the inverse model of additive noise and estimation error; means for reconstructing an image of the subject by applying the inverse model to the pMRI signal data; and means for providing a reconstructed image output. Preferably, the means for receiving pMRI signal data, the means for constructing a system model, the means for determining an inverse model, the means for reconstructing an image, and the means for providing a reconstructed image output, all comprise computer software and/or hardware components. For example, the means for receiving pMRI signal data may include a suitable input peripheral interface of a computer system, along with associated software. Suitable input interfaces include network interfaces, whereby the pMRI signal data is transferred via a data communications network, direct serial or parallel interfaces to an MRI machine, storage device interfaces for retrieving the signal data from corresponding storage media, and/or a user interface device such as a keyboard. Similarly, the means for providing a reconstructed image output may comprise an output peripheral interface of a computer system, such as a display, printer, or the like. Alternatively, or additionally, a reconstructed image, in digital form, may be output via a network interface or a storage device interface for external storage, display and/or processing. Various suitable computer arrangements, including conventional personal computer (PC) hardware, will be readily apparent to those skilled in the art.

The means for constructing a system model, determining an inverse model, and reconstructing an image will typically comprise software components adapted for performing suitable signal processing and/or numerical computation functions. Suitable means for implementing such software components, for any given embodiment of the invention, will be readily available to those skilled in the relevant art.

In a further aspect, the present invention provides an apparatus for reconstructing an image of a subject within a magnetic resonance imaging (MRI) machine having a plurality of receiver coils, wherein accelerated imaging is performed by reducing a number of distinct MRI measurements performed while simultaneously receiving measurement data from said plurality of receiver coils, the apparatus comprising: at least one processor; at least one input peripheral interface operatively coupled to the processor; at least one output peripheral interface operatively coupled to the processor; and at least one storage medium operatively coupled to the processor, the storage medium containing program instructions for execution by the processor, said program instructions causing the processor to execute the steps of: receiving via the input peripheral interface parallel magnetic resonance imaging (pMRI) signal data comprising digitised samples acquired from the plurality of receiver coils of the MRI machine; constructing a dynamic input-output system model of a pMRI detection process of the MRI machine, said model at least embodying a plurality of estimated spatial sensitivity functions corresponding with each of said plurality of receiver coils; determining an inverse model corresponding with said system model in accordance with an optimisation method which is adapted to minimise a measure of reconstruction error, wherein said measure of reconstruction error accounts for accuracy of image reconstruction by the inverse model in the absence of additive noise and estimation error, in combination with amplification by the inverse model of additive noise and estimation error; reconstructing an image of the subject by applying the inverse model to the pMRI signal data; and providing via the output peripheral interface a reconstructed image output.

"Comprises/comprising" when used in this specification is taken to specify the presence of stated features, integers, steps or components but does not preclude the presence or addition of one or more other features, integers, steps, components or groups thereof.

Further preferred features and advantages of the invention will be apparent to those skilled in the art from the following description of preferred embodiments of the invention, which should not be considered to be limiting of the scope of the invention as defined in the preceding statement, or in the claims appended hereto.

BRIEF DESCRIPTION OF THE DRAWINGS

Preferred embodiments of the invention are described with reference to the accompanying drawings, wherein:

Figure 1 is a block diagram representing pMRI image reconstruction as a model-matching problem in accordance with preferred embodiments of the present invention;

Figure 2 is a schematic diagram illustrating a processing system suitable for implementing embodiments of the invention;

Figures 3A and 3B illustrate results of simulated pMRI image reconstruction in accordance with the prior art and with a preferred embodiment of the present invention, respectively; Figures 4A, 4B, and 4C illustrate results of image reconstruction resulting from an MRI scan in accordance with the prior art and with an embodiment of the present invention; and

Figure 5 illustrates further results of MRI image reconstruction in accordance with the prior art and with an embodiment of the present invention. DETAILED DESCRIPTION OF PREFERRED EMBODIMENTS

In accordance with embodiments of the present invention, the pMRI image reconstruction is cast as a model-matching problem. A block diagram 100 of one such embodiment, representing the reconstruction problem in image space, is depicted in Figure 1.

As illustrated in the block diagram 100, an excited spin density function P is the (two-dimensional) response function of the excited nuclei within a selected slice of a subject volume under MRI imaging. The excited spin density function P (102) is detected by a plurality of receiver coils, collectively having spatial sensitivity functions represented by the blocks C (104) and ΔC (106). The spatial sensitivity functions are not precisely known, and must be estimated. Accordingly, the system model includes an estimated part C (104) and an unknown part ΔC (106) which is the estimation error. The detected signal is also accompanied by additive noise μ (108). A resulting pMRI signal S (110) incorporates the desired image function

P (102) modified by the estimated spatial sensitivity functions C (104) and "corrupted" by random noise μ (108) and noise υ(109) which is generated by the estimation error ΔC (106) driven by the desired image P (102).

The model-matching problem is thus to find the transfer function (matrix) F (112) which generates an "optimal" reconstruction P (114) of the "true" image P. The model represented by the block diagram 100 includes an error signal

^ = P-P (116). In a real pMRI image reconstruction system, the error signal ε cannot, of course, be measured directly, since the input P (102) is unknown, but it can be related to C, ΔC and F to derive optimal reconstruction of P. In the general case, in which the estimation error ΔC (106) and noise μ (108) are unknown, perfect reconstruction (PR) is not possible, and thus the error signal ε is non-zero Consider an embodiment of an MRI system with L receiver coils, and pMRI acceleration M, such that L > M. It may be shown, in accordance with a preferred approach, that this system is equivalent to an over-sampled parallel filter bank with uniform decimation. Image reconstruction in pMRI may be considered equivalent to signal reconstruction in a cyclic filter bank.

In accordance with preferred approaches, signal reconstruction may therefore be performed either in the image domain using aliasing component (AC) matrices, or in k-space using polyphase matrices.

For image domain reconstruction, the analysis AC matrix has the form:

while the synthesis AC matrix F has the form:

0, 1,-, L-I (3)

In k-space, an analysis polyphase matrix E may be denoted as :

0, l,---, M-I (4) while a corresponding synthesis polyphase matrix H may be denoted as:

0,\,-, L-I (5)

The relationship between the analysis AC and polyphase matrices is given by:

Ε(«) = C(n)WΛ (6) where Λ = diag{\, W^, - , wf^~'⁾ⁿ }, and W is the M x M DFT matrix.

In accordance with the AC matrix (image domain) approach, the image P is constructed as:

where i^>(n):= [p(n), P{n + K), -, P{n + {M -I)K)J (8) is a blocked version of P , and

is the received pMRI signal assembled from the signals received by the L receiver coils. According to the polyphase matrix (k-space) approach, a Fourier transformed image may be reconstructed using: p(M:) = H(n)s(M;), O≤k≤K-l (10) wherein

(11 ) and s(Mc) := [S₁ [Mk), S₂ [Mk), -, S_L [Mk)J (12)

The perfect reconstruction (PR) conditions for the polyphase and the AC matrix representations are respectively:

F(n)C(n) = I_M (14)

Considering the image-space (AC matrix) approach, which may be understood by direct reference to the block diagram 100 illustrated in Figure 1 , in the absence of any estimation error ΔC (106) and noise μ (108) the reconstructed image P is given by

P = FCP (15)

Accordingly, the reconstructed image P (114) will be a perfect reconstruction of the input image P (102) if F = C⁺ , where C⁺ is a pseudo-inverse of C. However, in the presence of both estimation error ΔC (106) and noise μ (108) the error signal ε is given by ε = (FCP + FΔCP + Fμ) - P

= (FC - I) P + F(ΔCP + /z)

In accordance with the prior art SENSE method, for example, a left pseudo

(least squares) inverse C⁺ is computed based upon the estimated C. In accordance with the model-matching representation illustrated in Figure 1 , even in the best case the residual error ε is non-zero, and depends upon the estimation error ΔC, the true image P and noise μ: ε = FΔCP + Fμ

= F(ΔCP + /z) (17)

= F(υ + μ) with the synthesis AC matrix F = C⁺. Since ΔC is the unknown part of the sensitivity functions, it generally has the same system structure as C(n) given in Equation (2), which is a dynamic input-output system. The input to this system is the true image P and the output is y = ΔCP(109), which is an additive noise input to the synthesis AC matrix C⁺.

As can be seen, from equation (17), the pseudo-inverse C⁺ operates on both the noise υ generated by estimation error ΔC and the random noise μ. Accordingly, the SENSE method is highly sensitive to the additive noises when the condition number of C is poor, resulting in the pseudo-inverse C⁺ having high gain. Unfortunately, this is quite likely to be the case when high acceleration is used for fast imaging.

In order to address this problem, preferred embodiments of the present invention provide a computerised pMRI reconstruction method for reconstructing an image of a subject within an MRI machine. In terms of the foregoing notation, the MRI machine has a plurality L of receiver coils, and imaging is accelerated by a factor M, /e the number of distinct MRI measurements is reduced by this factor while simultaneously receiving measurement data from the plurality L of receiver coils. In practice, therefore, L > M.

In general, the method commences by a computer receiving pMRI signal data including digitised samples acquired from the L receiver coils of the MRI machine. A system model of a pMRI detection process of the MRI machine is constructed. In particular, the system model embodies at least a plurality L of estimated spatial sensitivity functions corresponding with each of the L receiver coils. For example, in accordance with the image domain model represented by the block diagram 100 of Figure 1 , the estimated special sensitivity functions are embodied within the analysis AC matrix C (104). As will be appreciated, the estimated spatial sensitivity functions may equivalents be embodied via other means, for example, a synthesis polyphase matrix H.

It will be understood that the spatial sensitivity functions of the receiver coils of the MRI machine may be estimated using various techniques known in the art. For example, the spatial sensitivity functions of the receiver coils may be estimated by pre-imaging, such as in the prior art SENSE and PILS methods. In the AUTO SMASH and GRAPPA methods, the synthesis polyphase matrix H is estimated directly during the actual imaging process, and techniques of this type may also be applicable to the present invention. As will be appreciated, there may be various other existing techniques, as well as methods yet to be developed, that could equally be applied within embodiments of the present invention.

Embodiments of the inventive method then proceed by determining, for example via the model-matching approach previously described, an inverse model corresponding with the system model. In terms of the image-space model represented by the block diagram 100 of Figure 1 , the inverse model is the transfer function (ie matrix) F (112), the purpose of which, as previously discussed, is to generate an optimal reconstruction P (114) of the subject image P (102). In accordance with the invention, an optimisation method is employed which is adapted to minimise a measure of reconstruction error, such as the energy or the variance of error signal ε (116) of the model depicted in Figure 1. Importantly, the measure of reconstruction error accounts for the accuracy of image reconstruction by the inverse model in the absence of additive noise and estimation error, in combination with amplification by the inverse model of additive noises υ and μ . That is, embodiments of the present invention seek to improve over the prior art by optimising not only over the accuracy with which an inverse function of the spatial sensitivity of the receiver coils is established, but also taking into account the extent to which the final reconstructed image may be corrupted due to amplification of additive noise and estimation errors, which is a drawback of existing methods, such as SENSE as described above. Once the inverse model has been determined, subject to the relevant optimisation constraints, the subject image is reconstructed by applying the inverse model to the pMRI signal data S (110). A reconstructed image output may then be provided to an end-user.

More particularly, in accordance with preferred embodiments of the invention, the step of determining the inverse model may be formulated as an HL norm optimisation problem. For example, the SENSE method may be improved utilising the teaching of the present invention by seeking a solution to the HL norm optimisation problem defined by: min{ ||l-F(n)C(n)||² _∞ -i-

(18)

wherein β is a weighting factor, which may be predetermined, or determined automatically by adaptive estimation procedures. As will be appreciated, the optimisation problem defined by equation (18) includes two components, the first of which corresponds with the PR condition addressed by the prior art SENSE method, while the second term corresponds with the gain of the transfer matrix F (112) which may be responsible for undesirable amplification of additive noises.

The I-L norm optimisation problem defined by equation (18) can be solved by Linear Matrix Inequality (LMI) optimisation, as now outlined.

Defining {Ac, Bc, Cc, D_c} and (AF, BF, C_F, DF} be the state space realisations of C and F respectively, then the state space realisation of the system [I-FC, F] has the following form:

It can then be shown that the solution to the problem defined by equation (18) can be found by the following LMI optimisation:

subject to

^~AA^r7PPAA--PP AA^T7PPBB_NN

P=P⁷ > 0 where γ_a and γ_b are the upper bounds on the first and second terms of the expression in equation (18), respectively. The foregoing LMI optimisation procedure may be performed either using a predetermined value of β, or using a value of β determined using any suitable adaptive estimation method.

It will be appreciated that various computational techniques are available for performing the optimisation defined by the foregoing equations. There is accordingly shown in Figure 2 a schematic diagram illustrating a processing system 200 which is suitable for the implementation of embodiments of the present invention. In particular, the system 200 is a computing system having at least one processor 202 interfaced to at least one storage medium 204, typically being a suitable type of memory, such as Random Access Memory, for containing program instructions and transient data related to the operation of the processing system 200, as well as the implementation of operations and methods embodying the invention. The storage medium 204 may also include persistent storage, such as magnetic or optical drive storage, for longer-term storage of programs and data. The processor 202 is also typically interfaced to a peripheral bus 206, or equivalent, which provides access to one or more peripheral interfaces 208. The peripheral interfaces 208 may include one or more input/output devices, including user interface devices such as a keyboard, pointing device, display and so forth, as well as communications interface devices, such as network interface devices, serial/parallel data interface devices, and so forth. As will be appreciated, peripheral interface devices 208 may be provided enabling the MRI signal data to be entered into the processing system 200, and to enable the output of a reconstructed image, for example on a display device, hard-copy device and/or via a communications data interface. The memory 204 further contains a body of program instructions 210 particularly implementing the various method steps, and computational/numerical procedures required in accordance with an embodiment of the invention, such as those developed in the foregoing description. As will be appreciated, various development and execution environments are available that would be suitable for implementation of such methods. These include native code for direct execution by the processor 202, but also encompass dedicated numerical and computational processing environments, such as that provided by the product MATLAB®. Various numerical and signal processing libraries are also available, compatible with a variety of programming languages and environments, which could be used in the implementation of embodiments of the invention. It will accordingly be appreciated that the invention is not limited to any particular means of implementation, in this respect. A processing system generally in accordance with that represented in

Figure 2 has been implemented and utilised to illustrate the effectiveness of the present invention, particularly in the embodiment represented by the block diagram 100 of Figure 1 , and equations (18) to (21) above. The H∞ optimised reconstruction obtained from this implementation has been compared with the prior art SENSE method, in the simulation described by the following.

In particular, an array of eight receiver coils has been simulated by numerical calculation (utilising the Biot-Savart law) in order to generate each (simulated) coil sensitivity function. The simulated sensitivity functions are then multiplied with a magnetic resonance image, and Fourier transformed to produce corresponding sampled k-space representations. The k-space representations have been down-sampled in order to simulate accelerated acquisition with M = 8. Additionally, additive noise μ has been simulated by adding white Gaussian noise samples to the down-sampled k-space data to produce an overall simulated pMRI signal S having a signal-to-noise ratio of 37dB. The above-described H∞ optimisation has been performed, to obtain the H∞ optimised inverse model F, which is compared with the least squares PR solution of the prior art SENSE algorithm. In each case, a relative reconstruction error ε_r is defined according to:

wherein Δlj,j and lj,j are respectively the absolute reconstruction error and the simulated input at a corresponding pixel of the MRI image.

Figure 3A illustrates the results of simulate pMRI image reconstruction 302, in accordance with the prior art SENSE method. Figure 3B illustrates the simulated pMRI image reconstruction 304 in accordance with the aforedescribed embodiment of the present invention. The improvement provided in accordance with this embodiment is clearly visually apparent in the images. This improvement is quantified by the following table of results.

As can be seen, the relative reconstruction error in accordance with the embodiment of the present invention is 15.82%, as compared with 81.48% for the prior art SENSE method. It is notable that while the reconstruction in accordance with the present invention does not approximate as closely the perfect reconstruction condition (9.9867x10^"\ as compared to 1.0247x10^'9 in accordance with the prior art method), overall a significant improvement is achieved due to the reduction, by approximately a factor of 100, in the gain (as measured by the Hoo norm) of the inverse function F (3.9529x10⁴ in accordance with the present method, as compared with 2.1665x10⁶ for the prior art).

While the images in Figure 3A and 3B are based upon simulated data, Figures 4A, 4B and 4C show images derived from a real-world MRI example. Figure 4A shows and MRI image 402 of a slice through the head of a human subject. The image 402 has been generated using four receiver coils without acceleration (M = 1 ).

The image 404 shown in Figure 4B corresponds with pMRI imaging having an acceleration factor of M = 4. The image 404 is reconstructed using the prior art (PR) SENSE method, and the quality and clarity of the image 404 are clearly inferior to the unaccelerated image 402. The effects of amplified noise in the reconstructed image 404 are also apparent.

Finally, the image 406 shown in Figure 4C corresponds with the above- described embodiment of the present invention, again utilising an acceleration factor of M = 4. The quality and clarity of the image 406 are more closely comparable to the unaccelerated image 402, and are significantly superior to the prior art pMRI reconstruction 404. Accordingly, taking into consideration the acceleration factor of 4, reducing the imaging time for the human subject by a factor of four, pMRI imaging in accordance with the present invention provides the potential for clear benefits over the prior art. In particular, embodiments of the present invention will enable relevant health professionals to select between a wide variety of tradeoffs between imaging time and image quality via pMRI techniques, while generally having improved outcomes as compared with prior art pMRI methods. Embodiments of the present invention may accordingly find application in a wide variety of situations, including those in which a dynamic picture of internal subject activity is required, including reasonably rapid real-time updates of MRI imagery. Figure 5 shows further images derived from another real-world MRI example, namely a cardiac image generated using an eight channel cardiac coil array. The image 502 has been generated using the SENSE method without acceleration (M = 1 ). Images 504, 508 correspond with pMRI imaging having an acceleration factor of M = 6, while images 506, 510 correspond with pMRI imaging having an acceleration factor of M = 8. The images 504, 506 have been reconstructed using the prior art SENSE method, whereas the images 508, 510 are H∞ optimal reconstructions generated in accordance with the above-described embodiment of the present invention. Once again, it is apparent that the quality and clarity of the prior art SENSE reconstructed images 504, 506 are markedly inferior to the unaccelerated image 502. The H_∞ optimal reconstructions 508, 510 are visibly superior, and many features of the images 508, 510 remain visible, as compared with the unaccelerated image 502, despite the significant acceleration factors of M = 6 and M = 8 respectively.

The observation that image reconstruction in pMRI may be considered equivalent to signal reconstruction in a cyclic filter bank (FB) enables additional improvements and benefits to be achieved in the case of auto-calibrating image reconstruction methods, such as GRAPPA. As has previously been noted, such methods do not use pre-imaging in order to estimate the spatial sensitivity functions of the receiver coils, but rather utilise auto-calibrating signals (ACS) during an initial, unaccelerated, measurement period. In particular, the ACS are acquired by full excitation of a sub-area of k-space, without downsampling, for a calibration period. As will be appreciated, the number of ACS lines used in the calibration process affects both image quality, and calibration time. The greater the number of ACS lines utilised, the more accurately the coil sensitivity functions may be estimated, but also the greater time that will be required for calibration.

For a given number of ACS lines, the particular lines chosen for calibration may be selected in order to optimise estimation, and in particular it can be shown that the use of central k-space harmonics provides an optimal choice, since these include low-frequency harmonics containing the most energy of the associated image signal.

However, it is desirable to provide improved methods, enabling superior image reconstruction for any given choice of number and position of ACS lines. As noted above, the observation that pMRI image reconstruction is equivalent to signal reconstruction in a cyclic FB enables such improvements to be achieved. For example, the GRAPPA method may be improved in accordance with preferred embodiments of the invention by reformulating its underlying reconstruction process within the FB framework. More particularly, the synthesis polyphase matrix, or FB, may be decomposed into L subsystems as follows:

It may then be shown that the synthesis FB subsystems may be estimated using

A key insight of the FB interpretation of pMRI reconstruction is that the conventional GRAPPA estimation effectively assumes that the FB is causal, and has finite impulse response (FIR). However, it is known that the synthesis FB may, in general, be noncausal infinite impulse response (MR), even if the analysis FB (Ze polyphase matrix E) has causal FIR characteristics. Accordingly, image artefacts will inherently be generated in the conventional GRAPPA reconstruction due to the artificial structural constraint that the FB is causal FIR. Preferred embodiments of the invention utilising ACS techniques for sensitivity function estimation accordingly avoid this artificial constraint, by utilising the following methods. Assuming that the synthesis FB takes a general noncausal MR structure, it may be decomposed into causal and anti-causal subsystems as follows:

H(n) = H_c(rc)+ H» (25)

In the foregoing equation, and in the subsequent analysis, the symbol "c" indicates a causal component, whereas "ac" indicates an anti-causal component. The following relationships then apply:

% [Mk

In equation (28), measured quantities are indicated by the tilde ("~"), the input s(Mk) is defined analogously to equation (12), and the output is: s, [Mk)

(29)

Additionally, the causal and anti-causal subsystems may be written in left matrix fraction form:

It is then necessary to determine the output of the causal and anti-causal synthesis subsystems. An approach similar to the conventional GRAPPA reconstruction is used, whereby it is assumed that all of the first row elements of both H_c(n) and H_ac(n) are equal to Vz. This is equivalent to the key assumption of GRAPPA, and indeed all ACS methods, namely that H_0/(n) = 1. Physically, this corresponds with an assumption that received coils are disposed in a spatially symmetrical arrangement, which is often the case in practical MRl machines. The abovementioned assumption gives:

- s: (Mc) = H^n) s (Mc) (32)

i s,(Mc) = Hl_c(n) s(M/c) (33) Accordingly, the outputs for both subsystems are the same, but with opposite indices. Again, this may be given a physical interpretation, namely that if one models the spatial response of a certain coil in the receiver array, then the amount of "prediction" contributed by its left coils (causal) and right coils (anti-causal) should be the same.

Equation (28) may be written in the form of the following linear parameter model

in which

and

wherein Ri, Qi and R₂, Q₂ are the orders of the causal and anti-causal subsystems respectively, and the A' and B' are the parameter coefficient matrices of subsystem i, as defined by equations (30) and (31).

The least squares solution for an optimal set of parameters is then given by:

Once all of the subsystem FBs defined by equations (30) and (31) have been estimated, the final output can be calculated using equations (32) and (33) and:

L-I p(Mc) = ∑_S/(M/c) (36)

/^•=1

The foregoing method has been tested using a human head dataset measured on a 1.5 tesla Siemens Avanto system with a four channel head coil. The matrix size for each channel is 224 x 256. Each channel's k-space data was downsampled to simulate accelerated imaging, with a number of the central k-space lines fully retained as the ACS data. The sum of squares of all channel images was utilised as a reference. An experimental analysis was then performed comparing conventional GRAPPA with a preferred embodiment of the present invention, with an acceleration factor of M = 2, and the number of ACS lines chosen as 112, 56 and 28 respectively. The order of the conventional GRAPPA reconstruction was chosen between 4 and 8, based on the reconstruction quality. The modelling orders utilised in the experimental embodiment (Ze the Rs and Qs) were experimentally determined according to a normalised reconstruction error calculated in a similar manner to equation (22). The resulting orders of the IIR model in all experiments were relatively low, with values of R in the range of 1 to 2, and values of Q in the range of 2 to 5.

A comparison of the reconstruction errors for the three values of ACS, between conventional GRAPPA and the experimental noncausal HR FB embodiment, is summarised in the following table.

It is apparent that the experimental embodiment outperforms conventional GRAPPA, and this improvement was also visible in the corresponding reconstructed images (not shown). As can be seen from the foregoing table, the experimental embodiment performed better than conventional GRAPPA in all cases, and is able to produce a relatively consistent, and low, level of reconstruction error even when the number of ACS lines is significantly reduced, from 112 to 28. Conventional GRAPPA, on the other hand, exhibits higher reconstruction error generally, and significantly increased reconstruction error when the number of ACS lines decreases. In summary, the experimental embodiment demonstrates the superiority of noncausal HR FB methods over conventional GRAPPA, resulting from removal of the structural constraints in the prior art methods. The use of a more general noncausal HR structure for identification of the synthesis FB gives more accurate reconstruction. Furthermore, it is anticipated that by extending embodiments of the present invention to multi-slice imaging, further improvements may be achieved, since additional data will be available for FB identification.

While preferred embodiments of the invention have been described herein, it will be understood that these are not intended to limit the scope of the invention, which is defined by the claims appended hereto.

Claims

CLAIMS:

1. A parallel magnetic resonance imaging (pMRI) reconstruction method for reconstructing an image of a subject within a magnetic resonance imaging (MRI) machine having a plurality of receiver coils, wherein accelerated imaging is performed by reducing a number of distinct MRI measurements performed while simultaneously receiving measurement data from said plurality of receiver coils, the method comprising the steps of: receiving pMRI signal data comprising digitised samples acquired from the plurality of receiver coils of the MRI machine; constructing a dynamic input-output system model of a pMRI detection process of the MRI machine, said model at least embodying a plurality of estimated spatial sensitivity functions corresponding with each of said plurality of receiver coils; determining an inverse model corresponding with said system model in accordance with an optimisation method which is adapted to minimise a measure of reconstruction error, wherein said measure of reconstruction error accounts for accuracy of image reconstruction by the inverse model in the absence of additive noise and estimation error, in combination with amplification by the inverse model of additive noise and estimation error; reconstructing an image of the subject by applying the inverse model to the pMRI signal data; and providing a reconstructed image output.

2. The pMRI reconstruction method of claim 1 wherein the step of constructing a dynamic input-output system model comprises embodying the spatial sensitivity functions in an image-space representation as an aliasing component (AC) matrix.

3. The pMRI reconstruction method of claim 1 wherein the step of constructing a system model comprises embodying the spatial sensitivity functions in a k-space representation as a polyphase matrix. 52

27

4. The pMRI reconstruction method of any one of claims 1 to 3 wherein the step of determining an inverse model comprises utilising an optimisation method which includes minimising, according to specified criteria, a measure of performance that comprises: (a) a term corresponding with a perfect reconstruction (PR) condition; and (b) a term corresponding with a gain of the inverse model.

5. The pMRI reconstruction method of claim 4 which further comprises (c) an Hoc norm metric of the said terms (a) and (b).

6. The pMRI reconstruction method of any one of claims 1 to 5 wherein the step of determining an inverse model comprises utilising an Hoc norm optimisation method.

7. The pMRI reconstruction method of claim 6 wherein the H∞ norm optimisation method utilises an adaptively estimated weighting factor (β).

8. A computer-implemented system for performing a parallel magnetic resonance imaging (pMRI) reconstruction method for reconstructing an image of a subject within a magnetic resonance imaging (MRI) machine having a plurality of receiver coils, wherein accelerated imaging is performed by reducing a number of distinct MRI measurements performed while simultaneously receiving measurement data from said plurality of receiver coils, the system comprising: means for receiving pMRI signal data comprising digitised samples acquired from the plurality of receiver coils of the MRI machine; means for constructing a dynamic input-output system model of a pMRI detection process of the MRI machine, said model at least embodying a plurality of estimated spatial sensitivity functions corresponding with each of said plurality of receiver coils; means for determining an inverse model corresponding with said system model in accordance with an optimisation method which is adapted to minimise a measure of reconstruction error, wherein said measure of reconstruction error accounts for accuracy of image reconstruction by the inverse model in the absence of additive noise and estimation error, in combination with amplification by the inverse model of additive noise and estimation error; means for reconstructing an image of the subject by applying the inverse model to the pMRI signal data; and means for providing a reconstructed image output.

9. An apparatus for reconstructing an image of a subject within a magnetic resonance imaging (MRI) machine having a plurality of receiver coils, wherein accelerated imaging is performed by reducing a number of distinct MRI measurements performed while simultaneously receiving measurement data from said plurality of receiver coils, the apparatus comprising: at least one processor; at least one input peripheral interface operatively coupled to the processor; at least one output peripheral interface operatively coupled to the processor; and at least one storage medium operatively coupled to the processor, the storage medium containing program instructions for execution by the processor, said program instructions causing the processor to execute the steps of: receiving via the input peripheral interface parallel magnetic resonance imaging (pMRI) signal data comprising digitised samples acquired from the plurality of receiver coils of the MRI machine; constructing a dynamic input-output system model of a pMRI detection process of the MRI machine, said model at least embodying a plurality of estimated spatial sensitivity functions corresponding with each of said plurality of receiver coils; determining an inverse model corresponding with said system model in accordance with an optimisation method which is adapted to minimise a measure of reconstruction error, wherein said measure of reconstruction error accounts for accuracy of image reconstruction by the inverse model in the absence of additive noise and estimation error, in combination with amplification by the inverse model of additive noise and estimation error; reconstructing an image of the subject by applying the inverse model to the pMRI signal data; and providing via the output peripheral interface a reconstructed image output.