[attachment=3772]

---

Presented BY
M.SHALEM RAJ
III ECE (05R81A0435)
SRI SUNFLOWER COLLEGE OF ENGG. AND TECHNOLOGY

---

MAGNETIC RESONANCE IMAGING (MRI) SIMULATION
ON A GRID COMPUTING ARCHITECTURE

Abstract
In this paper, we present the implementation of a Magnetic Resonance Imaging (MRI) simulator on a GRID computing architecture. The simulation process is based on the resolution of Bloch equation in a 3D space. The computation kernel of the simulator is distributed to the grid nodes using MPICH-G2. The results presented show that simulation of 3D MRI data is achieved with a reasonable cost which gives new perspectives to MRI simulations usage.
INTRODUCTION
The simulation of Magnetic Resonance Imaging (MRI) is an important counterpart to MRI acquisitions. Simulation is naturally suited to acquire theoretical understanding of the complex MR technology. It is used as an educational tool in medical and technical environments. By offering an analysis independent of the multiple parameters involved in the MR technology, MRI simulation permits the investigation of artifact causes and effects. Likewise simulation may help in the development and optimization of MR sequences.
With the increased interest in computer-aided MRI image analysis methods (segmentation, data fusion, quantization...), there is a greater need for objective methods of algorithm evaluation. Validation of in vivo MRI studies is complicated by a lack of reference data (gold standard) and the difficulty of constructing anatomical realistic physical phantoms. In this context, an MRI simulator provides an interesting assessment tool since it generates 3D realistic images from medical virtual objects perfectly known.
In 1984, Bittoun et al. built one of the first 1D MRI simulators by numerically solving the Bloch equations at each point in an object. The same technique was later applied to 2D and 3D objects by Summers et al. and Olsson et al. Brenner et al. have proposed a distributed implementation of these computationally intensive simulations. In the domain of brain imaging, Kwan et al have limited the computation time by using templates of the main brain tissues. Nevertheless, MRI simulation is still under investigation in order to simulate high resolution MR images including MRI artifacts (Magnetic susceptibility, chemical shift, diffusion ¦) that may impact the quantitative measurements required for several MRI applications.
In this context, we develop, a 3D MRI simulator named SIMRI that is designed to simulate realistic high resolution 3D MR images. Since simulation of the MR physics is computationally very expensive, parallel implementation is mandatory to achieve performances compatible with the target applications. By offering a virtually unlimited computing power, grid technologies appear to be promising.
This paper focuses on the grid implementation of this simulator. This paper outlines the main components of the simulation process and details the "gridification" of the computation kernel and presents the main results.
MRI SIMULATION
Simulator overview
MR imaging involves the nuclear spin system of an object to be imaged, the magnetic fields generated by the imager, the resonance and relaxation phenomena arising from the interaction nuclei / magnetic fields, and the signal detection and reconstruction. All these components are taken into account by the SIMRI simulator (Figure 1). From a virtual object and a MRI sequence, the magnetization computation kernel computes the RF signals named k-space. The final MR image is reconstructed from the k-space with generally a 3D fast Fourier transform.

Virtual object description
The virtual 3D object is a discrete description of a real object spin system. Each voxel of the virtual object contains a set of physical values necessary to computes with the Bloch equations the local spin magnetization. These values are the proton density noted r and the two relaxation constants T1 and T2. Basically, a virtual object describes the nuclear spin system of one isochromat, the water hydrogen atom. Note that SIMRI is able to deal with many isochromats in order to simulate chemical shift artifact (Figure 2) or partial volume effect. In this case, each voxel contains the , T1, T2 values of each isochromat.
A local main field homogeneity is also attached to each virtual object voxel. This value is used to simulate inhomogeneous main field and/or susceptibility artifact.
A virtual 3D object of size S= X.Y.Z voxels is then described by (3p+1).S values where p is the number of isochromats modeling to each object point.
Figure 2: 256² simulation of the chemical shift artifact using a two isochromats (water and fat) virtual object. The white and black part at the left and the right of the central circle corresponds to the artifact.
MRI Sequence
During an MRI experience, an object is placed in a static magnetic field named B0 and is excited by electro-magnetic events of two types: Radio Frequency (RF) pulse referred as
B1 field and magnetic fields gradients. Excitation steps are followed by an RF acquisition step. It corresponds to the MRI system antennas acquisition of the magnetization state of the object which is stored as a complex signal in the k space [15]. The time chaining (Figure 3) of the magnetic events with RF signal acquisitions is named MRI sequence. It defines the way of filling the k-space and determines the final image characteristics. It must be noted that a MRI sequence in the simulator has the same components and the same time schedule than in a real system.

Magnetization kernel
The magnetization computation kernel called during the 3D Bloch equation [1]. It is a linear differential equation simulation of an MRI sequence is based on the solving of describing the time evolution of the spin magnetization and it is given by:
where is the spin magnetization vector, M0 is the spin magnetization equilibrium value which depends of the proton density are the relaxation constants and is the gyro magnetic constant of the considered isochromat (42.58 MHz/T for the Hydrogen proton). The magnetic field B is modeled as follows:
where B0 is the main static magnetic field, is the field in homogeneities, is the applied field gradient, is the RF pulse and is the spatial position.
The simulation kernel implements a discrete time solution [9] of the Bloch equation with the means of rotation matrices and exponential scaling which depend of the magnetic events of the MRI sequence. The magnetization vector evolution is computed from its previous time value according to the following equation:
where is a rotation matrix and it is defined by:
where is linked to the applied gradient by:
where is linked to the field in homogeneities by:
where R relax describes the relaxation effect by:
and where RRF represents the rotating effect of an RF pulse of phase angle leading to a flip angle defined by:
To design an MRI sequence, the magnetization computation kernel (based on equation 3) offers three functions which are application of RF pulse, application of gradient, and waiting which corresponds to the spin magnetization relaxation.
One can notice in equation 3 that the magnetization vector at a given localization depends only on the previous magnetization vector at the same position. This property will be used by the gridification of the simulation kernel.
The last functionality offered by the kernel is the RF signal acquisition (Figure 3) which corresponds to the acquisition by two antennas in the x, y plane of the magnetization state of the object after a given excitation. The RF signal is a one dimensional discrete complex signal that will fill in respect with the excitation sequence one line of the k-space volume.
One point k of the RF signal is obtained by summation of the local magnetization over the entire virtual object (Eq. 9). The next point is obtained after an evolution of the local magnetization respecting equation 3 with a time step DT equal to the sample period of the signal. This linearity property will be exploited by the grid architecture.
From equation 3 and equation 9, it appears that the simulation of one magnetization evolution is proportional to the virtual object dimension (X.Y.Z), while the simulation of the whole k-space acquisition is proportional to the object size and to the MRI image size (N.M.P). The whole simulation time of a 3D sequence (Ts) is then given by:
where c is a constant and Nex is the number of magnetic events simulated before an acquisition. In a MRI sequence, Nex is negligible comparing to the MRI image size. So, as an example, multiplying by two all the dimensions of the virtual object and the MR image leads to a simulation time multiplication by 16 in two dimensions and by 64 in three dimensions. If the simulation of small images (128x128 pixels) is done within a minute on a Pentium IV PC, equation 10 shows that high resolution MRI 3D simulation is not practicable on such a computer. Therefore, we turn toward Grid technologies that promise a virtually unlimited computing power.

GRID IMPLEMENTATION AND RESULTS
The whole code of the SIMRI simulator is written in ANSI C language. One module concerns virtual objects manipulation, one set of modules implements MRI sequences, one module deals with image reconstruction and the final set of modules is related to the magnetization computation kernel. Profiling the execution of a SIMRI simulation indicates
the percentage of time spent in each part of the simulation code. It appears clearly that most of the time (more than 90%) is spent in the computation kernel specifically on the rotational and exponential scaling needed by the Bloch equation resolution. Consequently, we focused our parallelization efforts on the kernel code in order to have a grid implementation of the simulator. This modification has to be almost transparent for MRI sequence developers .The gridification of the kernel must be done using the MPI [2] version of GLOBUS [16] since the middleware targeted for this application is the European Data Grid project [13] built on top of GLOBUS.
Gridification of the magnetization computation kernel
Because all the spin magnetization vectors are independents (Equation 3) and because the signal acquisition process is linear (Equation 9), many solutions for parallelization can be considered. One solution consists in distributing a minimal task which is to compute the new magnetization sate of a single spin vector. In this case, the granularity is maximized but the time of node communication will be predominant compared to the computation time.
The second approach consists in distributing the magnetization computation of a subset of spin vectors. This subset can be fixed to a given size or adapted to the number of active nodes. We have selected this last solution and since each computation node has a set of magnetization vectors, this node will compute the contribution of these vectors to the RF signal generated during an acquisition step. The data and process distribution to the grid nodes is given in Figure 4. All the computation nodes have the MRI sequence knowledge and they receive from the master node a part of the virtual object. They compute the magnetization evolution of the corresponding spin vector subset. At the end of each acquisition step, the master node collects and adds all the RF signal contributions and stores the RF signal in the k-space. At the end of the MRI sequence, the master node applies the reconstruction algorithm to generate the MRI simulated image. This solution is scalable and follows a simple parallelization scheme of type "divide & conquer" with a regular problem of space division. Two types of activities are defined: computation of subset for N-1 nodes and results collect for one node. The time dedicated to the node communication is minimized since each node received only one task. Nevertheless, this solution drawback is a low granularity introducing loss of efficiency on heterogeneous grid. Indeed, all the computation nodes have the same amount of data to process and then the lowest node will give the timing.
Results
The results presented here have been obtained on a 8 nodes cluster. Each node is a Pentium III-1Ghz with 1Go of RAM. We run the MPICH-G2 [2] version of MPI associated to Globus 2.2 [16] for the grid access control.
Figure 5: Nodes activity reported by XMPI during a 2562 Horizontal bars corresponds to the 8 nodes activity2D spin echo simulation using MPI instead of MPICH. between 0 and 16.94 seconds. Vertical white bars correspond to MPI message exchange at the end of RF signal acquisition.
First of all, we have compared the efficiency of using our grid compared to a single machine. When the virtual object size is over 210, the average efficiency of each grid node is over 87 % meaning that the simulation time is divided by 7 at least. Since one node out of eight is just collecting results and reconstructing the image, this means that the seven computation nodes have efficiency close to 100 %. This very good efficiency is due to the homogeneity of our grid and to the kernel gridification chosen which minimizes the number and the size of the exchanged messages as shown in figure 5. For a 2562 simulation, during an acquisition step, messages exchange and nodes synchronization takes 0.13 s compared to 2.5 s of computation time which corresponds to 5.2 % of the computation time. Table 1 shows that when the size of the virtual object increases this percentage decreases until a limit of 1.5 %. This limit is due to one node of our grid which is systematically late compared to the others.
The grid implementation of our simulator have been tested on 2D and 3D MRI Spin echo sequences applied on synthetic virtual objects of different sizes (Figure 6-7). Simulation time results are given in Tables 2 and 3. One can note that after a given size (1282 and 643) the time factor induced by the dimension change follows the rule of equation 10. Before these sizes, the communication and synchronization time is too important compared to the process time (Table 1).
These simulation results obtained are very encouraging since it shows that with a small cluster, MRI simulation of high resolution (1024²) 2D objects is possible within three days. Concerning 3D images, it is not realistic to simulate on such a small set of processors over 643 cubic MRI images. Nevertheless, it is possible to simulate within a week 3D multislice images (16 slices of 512x512 pixels). The simulation of high resolution 3D images should be tractable on full scale grids like the ones actually in development.

Figure 6 Figure 7
CONCLUSION
In this paper, we present a grid implementation of a 3D MRI simulator based on Block equations. The gridification is done in the magnetization computation kernel of the simulator and makes it transparent at the sequence level used by MRI sequence developers. The gridification strategy makes the simulator very efficient on a homogeneous grid. Results obtained on an 8 nodes cluster show that high resolution MRI simulation is possible at a reasonable cost. This is very encouraging for the different researches using MRI simulation like artifact analysis and correction or MRI post-processing [12]. Further works under development include analysis of MRI simulation performance on a larger grid and tests of different strategies for gridifying the computation kernel. Evolutions of the MRI simulation model are also under consideration in order to simulate more complex physical phenomena.
REFERENCES
http://globus.org
http://creatis.insa-lyon.fr/medigrid
http://edg.org
http://wwwsop. inria.fr/aci/grid/public/