When trying to employ the described iterative strategy for the reconstruction of real MRI data instead of simulated data, two additional difficulties arise. First, the observed object is usually complex-valued due to off-resonance effects and other phase perturba-tions. Although the desired result in the end is a real-valued image, it is impossible to remove the phase variations before resolving the spatial encoding. The removal of the phase is conventionally accomplished in a last step by calculating the magnitude
im-age, but this operation is non-linear and can not be integrated into the system matrix.
Therefore, it remains necessary to deal with the complex-valued nature of the object during the reconstruction procedure. Second, modern MRI systems use an array of re-ceive coils each having a varying intensity and phase profile, as described before. Again, a suitable combination of the signals from all individual coils can only be achieved after spatially resolving the object. As a consequence, any image reconstruction approach has to cope with the separate coil signals.
To this end, an iterative two-step reconstruction approach is proposed. The first step attempts to estimate sensitivity profiles for all coil elements, which are then utilized in the iterations of the second step to combine the individual coil channels and to remove the phase variations each time when mapping between the frequency domain and image space. Thus, the second step renders a combined and real-valued image. Noteworthy, because in this step a single image is found that complies with the observations from all individual coil elements at the same time, the method implicitly makes use of the locally varying coil sensitivities to gain additional object information in a similar manner to parallel imaging techniques like SMASH, GRAPPA, or SENSE [113,114,115,106]. How-ever, while most of these techniques try to extract the information in a concrete form (for example, as synthesized k-space samples), in the present approach the exploitation of the coil sensitivities manifests as an improved conditioning of the optimization prob-lem. Moreover, an attractive feature of the approach is that neither reference data nor shared data from prior measurements is needed to estimate the coil profiles, whereas in most parallel imaging methods additional acquisition steps are required for this task.
In the first step of the approach, the signals from all coil channels are handled separately and, thus, a set of individual images is calculated. Here, the real and imaginary parts are treated as independent parameters, leading to a complex-valued image estimate for every coil channel. It is well-known that MRI coil profiles are smooth functions that vary only slowly and do not have sharp edges. This knowledge is incorporated by using a quadratic regularization of the image derivatives, which leads to globally smooth images as discussed in the foregoing section
Rcoil(x) =X
i
(Dx(1)xi)2+ (D(1)y xi)2 , (5.14)
where Dx(1) and Dy(1) are the derivative operators of first order according to Eq. (5.11).
After finishing the iterations for all coils, a sum-of-squares image is computed. A division of the single channel images by the sum-of-squares image yields the respective coil profiles. Noteworthy, the estimated coil profiles also include other phase variations inside the object as the real-valued sum-of-squares image has been taken as a reference.
Because the penalty functionRcoil(x) depends quadratically onx, the line search, which
Figure 5.5: Schematic diagram of the procedural implementations of (left) the system matrix A and (right) the adjoint system matrix A] that are used to map from image domain to frequency domain and vice versa. For details see text.
is part of the conjugate gradient iteration, requires only one step and, thus, only a low number of function evaluations is needed to obtain a reasonable image. Therefore, the coil profile estimation step takes only moderate computational time.
For the second, combined reconstruction step, the raw data from all coil channels is stacked into one data vectory. Further, the system matrixAis extended by a multipli-cation with the corresponding coil profile before performing the Fourier transformation for every channel. Figure 5.5 shows a flow diagram of the operations that are executed by the system matrix A and the adjoint matrix A] to map between frequency and image space. By combining the data from all coil channels into the data vector y, the algorithm now renders one image estimate that matches to the observations from all coils. Moreover, removing the phase variations with the use of the estimated coil pro-files allows to discard the imaginary part of the image estimates and enables to apply constraints on negative values, which otherwise would not be possible. Noteworthy, a combined coil reconstruction also ensures that the TV constraint remains applicable in a multi-coil setup. Otherwise, the intensity modulation from the coil profiles would conflict with the idea of piecewise-constant images.
Outline of the Complete Algorithm
To estimate the coil sensitivity profiles, individual complex-valued imagesxj are recon-structed for each channelj by iteratively solving the optimization problem
xj = argmin whereyj denotes a vector containing all measured values from channel j. The penalty functionsRFOV andRcoil are used as defined in Eq. (5.6) and (5.14), where the weights of the penalty termsλFOV and λcoil have to be adjusted according to the value range of the measured data (see next section). The system matrixA comprises three operations
A= GKB· F ·M−1KB, (5.16)
where M−1KB denotes the pre-compensation for the roll-off effect according to Eq. (4.9), F is the two-dimensional fast Fourier transformation, and GKBdenotes the interpolation to the sampling locations in k-space using a convolution with the Kaiser-Bessel window and subsequent evaluation along the trajectory. A sum-of-squares image is calculated from these images with
which is then used to obtain the complex-valued coil sensitivities
Cj =C−1sos·xj . (5.18)
Here, · denotes a component-wise multiplication of the two vectors. In the final stage, the real-valued image x is calculated by solving the optimization problem
x= argmin where Aj denotes that the system matrix in the final stage includes a multiplication with the jth coil profile according to
Ajx=A(Cj· x). (5.20)
The penalty functions Rpos and RTV2 are given by Eq. (5.8) and (5.13). Because the total-variation term RTV2 depends on the modulus of the estimate components, a non-linear optimization technique is required for solving Eq. (5.19). When utilizing a technique based on the conjugate-gradient method, as done in the proof-of-principle implementation, it is necessary to evaluate the cost functions (i.e. the functions inside the brackets of Eq. (5.15) and (5.19)) as well as their gradients. The gradients of the
`2-norm terms can be deduced from Eq. (5.4), and the gradients of the penalty terms can be obtained by deriving the penalty functions with respect to all components of the estimate vector.