The estimation of a source of activity in the brain is much more useful if it is combined with structural information of the brain. Furthermore, the structural information can help in creating better models of the head. This information can be obtained from Magnetic Resonance Images (MRI). However, for these images to be useful, they first need to be processed. This chapter will describe the data processing of MRI-scans for one purpose only: The processing is dedicated for electrical and magnetic source imaging and is designed such that the processing is as simple as possible and can be carried out fully automatically, without any supervision.
The processing steps as developed consist of the automatic segmentation of the data set, the generation of models of the head as a volume conductor, and of the generation of 3D-views of the head. In the segmentation step, the information in the MRI is classified according to its type of tissue. This step is necessary to enable the generation of both the surface of the head, interfaces between tissues, and the 3D-views. This step allows the hydrogen-density image to be converted to a map of the areas in the head with different conductivities, as used in the inverse calculations. Usually the head is described by a compartment model where each compartment has a homogeneous, isotropic conductivity. The compartments considered in the processing of MRI-data are the brain, the cerebro-spinal-fluid, the skull and the scalp because the conductivities of these tissues differ significantly. The ratio of these conductivities is of the order of 1 : 3 : 1180 : 1. The most often used compartment model for the head consists of spheres, which are concentric with a sphere which is locally fitted to the surface of the scalp or to the inside of the skull. This sphere can be obtained from the MRI and so can the thicknesses of the various shells. The segmentation allows also the use of more realistically shaped models of the head, which are expected to improve the estimation of the source of the electric and magnetic fields. A triangulation of the surface of the head and of the various interfaces between the compartments can be used in the Boundary Element Method. Finally, the 3D-views can be used to display and evaluate the source in various ways. In order to decide whether the orientation and location of an equivalent dipole is acceptable, it is necessary to know how the cortex is folded at that location. This is only possible if in addition to the transversal slice one has to have other views like a sagital or coronal slice to ones disposal. 3D-views can further assist to identify the structure which contains the source of activity. This chapter will describe all these fully automated processing steps in detail.
When taking an MRI-scan, the subject is placed in a magnetic field with a strength usually between 0.5 - 1.5 T. This field aligns the atoms with a magnetic moment, such as hydrogen atoms. The hydrogen atoms are then perturbed by an rf signal. As the atoms return to their equilibrium position, they emit their own rf signal, with a frequency that depends on the local magnetic field strength. By introducing variable gradients in the applied fields, and taking phase shifts in the signal into account, the hydrogen density within a particular volume can be reconstructed from the rf-signals. The information obtained in this way is a measure of the amount of magnetization of the hydrogen in each volume element. This value differs for different types of tissues due to differences in the hydrogen density and its chemical bonds.
The two most important types of scanning techniques are the T11-weighted and
the T2-weighted. Both techniques, or combinations of them, give varying
contrasts between different tissues. For instance, in a T2-scan the
cerebrospinal fluid (CSF) is highlighted, but in a TI-scan it is hardly
distinguishable. On the other hand a TI-scan shows a clear difference
between the white matter and the gray matter, whereas a T2-scan does not.
For our purposes the complete head is scanned using a 3D TI-contrast enhanced
scan. The data of this scan can be processed without much difficulty
leading to a naturally looking display of the brain. Our scans are made on
a Philips T15 Gyroscan at the local hospital, Medisch Spectrum, in Enschede.
A complete data set consists of 256 transversal slices, each containing 256x256
voxels of lxlxl mm. The slices are contiguous.
The scanner uses a VAX/VMS system to reconstruct the data set. The data is
transformed to the ARC-NEMA standard (National Electronics Manufacturers
Association, 1989). In this standard each slice is described in one file,
so that the complete set consists of 256 files, plus one additional file with
personal information. This takes up approximately 22 Mbytes of disk
space. The data is then written an tape, and read into our computer, the
VAXstation 3520. The system is used for processing both MRI-data as well
as MEG/EEG data.
Slices can be displayed by translating the value of a voxel, which is usually scaled between 0 and 4095, into a gray value. This is not necessarily a linear sr-ale. Contrast can be enhanced in certain regions by using a non-linear scaling. When looking at the voxel values of a data set, it is obvious that a certain type of tissue does not have a unique range of values. Moreover, this range of values can shift throughout the data set, which makes it difficult to classify the voxels to the type of tissue.
Whereas a human being does not have any difficulty discriminating the brain from other parts in an MRI slice, this is not easy for a computer. To analyse the slices by computer, a histogram of each slice has to be made. In a histogram of a slice the number of voxels having a certain value is displayed (see figure IV.1)
The histograms show peaks for voxel values which occur relatively more often than others in the slice. A (large) peak which represents the background is always present. This consists of voxel values outside the head which are black in an image of the slice. These are low-valued voxels. The head consists of all the pixels except the background ones. This would suggest that all voxels with a value above a certain threshold can be considered to belong to the head. This works reasonably well, except that due to noise in the image, some pixels outside the head are also identified as part of the head.
The histograms of most slices show, apart from the background peak, two more peaks. In slices with poorer contrast both peaks are merged together. All the voxels of the gray matter in the slice have their values in the area of the first peak, all the voxels of the white matter in the area of the higher valued, second peak. Because also a number of voxels which are part of the skin and other tissues have a value within the range of these peaks, thresholding cannot be applied. To take into account the differing voxel values for certain tissues in different slices the histogram is analysed by choosing certain "keyvalues", such as the high end of the background peak, the lower end of the first peak, the high end of the second peak and, if possible, a value between the peaks. To obtain these values the histograms are smoothed and the derivative is determined, with peaks at the beginning and end of each of the original peaks, and zero at the extreme (see figure IV.2). The keyvalues are found, by scanning this derivative for appropriate values. These values are used in the segmentation of each slice.
The histogram is calculated for each slice of a data set and the keyvalues
are determined in the method described above. The first step in the actual
segmentation of a slice is to identify which voxels are part of the head and
which are not. This is carried out by selecting all voxels with a value
higher than the lower end of the first peak (one of our keyvalues). These
voxels will be the seeds of a regiongrowing operation, where neighbours of a
voxel which is already part of a set, are examined. If such a neighbour
satisfies particular criteria it is added to the set. This step is
repeated until no more voxels can be added to the set. At the start of the
procedure the set only consists of the seeds.
Note that a voxel's neighbour has to be defined. Each voxel can have
either four or eight neighbours. The speed of the region-growing operation
depends on the order in which the voxels are scanned. A better performance
in speed is usually reached if the order is reversed after each pass.
IV.3.1 The head
The criterion used for this region-growing operation is that the voxel value
has to be higher than the value at the upper end of the background peak in the
histogram. In this way the complete head is found, without noise elements
around it. Because such noise elements have no connection with the head
elements the region-growing operation cannot reach them. In case of severe
noise in the picture it is possible that voxels outside the head are selected as
seed points. Possible solutions are to set the seed-selection criterion to
a higher value, to perform an erosion after seed selection or to check
afterwards for isolated structures in the image. Performing an erosion
means that, from a set of voxels, all those that have a neighbour which is not
part of the set are removed. This ran be repeated more than once. Of
the above mentioned solutions the last option is used for the removal of
isolated structures. The risk of using the first two options is that the
complete head is not obtained if certain structures in the head are not
connected in a slice.
After the region-growing, the complete image is checked for 'holes' i.e., for
voxels which are not part of the set which makes up the head, but that are
completely surrounded by voxels that do. These voxels are also considered
to be part of the head. In practice, this applies mainly to the dark
points which represent the skull.
IV.3.2 The brain
The next step is the determination of the voxels which represent the
brain. First the voxels which are part of the head and have values between
the two following keyvalues are selected. The lowest keyvalue is found in
the histogram at a quarter of the distance between the upper end of the
background peak and the lower end of the first peak. This keyvalue is
called the lowerbound (see figure IV.2). The second keyvalue is at the upper end
of the second peak. In this way almost the entire brain and large parts of
the skin are selected. At least two erosion operations are performed on
this set of voxels. The purpose of this is to disconnect the brain voxels
from the skin voxels. The set of voxels obtained in this way, will act as
a 'mask'. A small region in the middle of the slice is chosen as the seed
under the condition that it is part of the mask. These voxels are usually
part of the brain. If the slice does not include any brain tissue, the set
of voxels making up the mask should have been empty, and no seed is selected.
A region-growing operation is again started, under the condition that the
neighbouring voxels are part of the mask. In this way region-growing
should not have been able to reach points outside the brain. If this is
not the case, because the erosion did not break all the connections between the
brain and the skin, this can easily be detected by the computer because all the
brain voxels should be surrounded by 'head-only' voxels. If this is not
the case, and there are voxels classified as brain on the outside of the head,
the operation has to be completely repeated, but now with an additional
erosion. This can be repeated as long as necessary or for a limited number
of times. If the number of attempts is too high, it has to be assumed that
there is no brain tissue in the slice in question, and that some other tissue
was initially assumed to be brain tissue. This happens mainly with
skin/fat tissue in a slice at the top of the head.
Next the number of erosions used to separate the brain from other tissues has to be compensated with the same number of dilations. A dilation is the addition of each neighbour of voxels that are already part of the set to that set, without any further restriction. The definition of "neighbour" is relevant in this respect, and the same definition has to be used in the dilation as was used in the erosion. Next, a region-growing operation is started, adding only voxels with values between the lowerbound and the lower end of the first peak to the region. This has to add the CSF pixels to the brain pixels if they were not already included. Any pixel inside the brain which is not yet considered as being part of it, should now be classified as if it were. Finally, a dilation followed by an erosion, also known as a closing operation, is performed to close small holes on the surface of the brain. In this way, the brain has a more or less smooth surface.
IV.3.3 CSF
All voxels classified as brain, but with a value lower than the start of the first peak are considered to be fluid.
IV.3.4 The skull
The last step is the detection of the skull. An upperbound value is established, which has a value halfway between the end of the background peak and the beginning of the first peak. Of all voxels which are inside of the head, but which are not part of the brain, those having a value below this upperbound are selected. An erosion removes single odd pixels and then a region-growing operation is again started, where the added voxels have a value lower than the upperbound. A closing operation completes this stage. It is well known that the detection of the skull from MRI images is difficult. Our segmentation does give usable results for modeling purposes, but turned out not to be good enough for 3D-images.
Figure IV.3: Segmented slice. Each tissue category is displayed in a
different color
Figure IV.3 is an example of the segmentation of a slice. Slices are always viewed from underneath. The segmentation information is coded in 4 bits. The original voxel values are contained in 12 bits and the 4 bits are added to them. The segmentation is usually carried out at night, needs no supervision and gives satisfactory to good results. These results can be improved afterwards by minimal user-editing, but in our research this never turned out to be necessary.
A graphical version of the program has been constructed to show the various segmentation steps for a single slice. This can be very helpful when improving the algorithms.
Points on the surface of the brain and the head can be generated from the segmented data set. These points are useful for fitting spheres in the head or for the generation of more realistic surface descriptions, which can be applied to solve the biomagnetic forward and inverse problems. When fitting a sphere locally, only the region of interest of the head is sampled, while for a surface description, the complete head has to be sampled.
A central point in the compartment of interest is selected in each slice. An imaginary line is drawn from this point towards the back of the head, voxel by voxel, until one is reached that is no longer part of the compartment of interest. This voxel is the first point taken on that slice. If the number of desired points on this slice is known, the angle between each line from the central point to the edge can be calculated. A new line is then initiated from the centre towards the edge. If it reaches the edge a new point is found. This repeats itself until points are found through 360 degrees (see figure IV.4).
Figure IV.4: Generation of surface points on a slice
Because the circumference of the head is smaller in some slices than in others the number of points per slice will have to vary. Therefore the slice that is approximately in the middle of the data set is sampled first. The length of the line from its central point towards the back of the head is an indication of the circumference of the head. By choosing the desired number of points on the slice, the length of the arc between points can be approximated. Looking at the length of the same line in other slices the number of points needed in that slice can be estimated.
All the points are written in a file. The order in which the points are found on the slice has to be the same for all slices. In our case, looking from the top, this order is anti-clockwise. The first points on each slice are aligned on the back of the head. A triangulated surface description can be generated from the points found in the described way. This is not very difficult to achieve, because both the head and the smoothed surface of the brain have convex shapes.
In order to triangulate the surface, two points are added, namely a top and a bottom point. Both can be considered to be the only point on a slice, and are not treated as exceptions. Their location is usualy in the middle of a slice, and the vertical distance to the next slice ran be taken as equal to the other inter-slice distances. The correct slice for the top point can be determined from the data set by tracing from the centre of the head to the top.
The triangulation is performed by looking at two adjacent slices, starting from the bottom or the top. The points of each slice form contours. Figure IV.5 demonstrates the principle. The first point of each contour is always at the back of the head. The first points of two adjacent contours (Ai and Ai+1) are connected by a line, the baseline. This baseline forms the start of a triangle. On each contour the next point is taken, and for each of these points (Bi and Bi+1) the distance to the point of the baseline on the other contour is calculated. The point with the shortest distance (Bi) is taken as the third point of the triangle. The other point is not used, but it becomes the 'forward point' (Bi+1), and is the most forward one in the direction of triangulation.
Now a new base-line has to be established. If the distance between the third point added to the triangle (Bi) and the old baseline point on the other contour (Ai+1) is smaller than the distance between this third point and the forward point, then this line will be the new base-line. If not, the line between the third point of the last triangle and the forward point will be taken as the baseline. The space skipped also forms a triangle consisting of one side of the previous triangle, the new baseline and a connecting line along the contour of the forward point (Ai+1, Bi and Bi+1). A new triangle is generated on the new base-line. The procedure stops when the first point on a contour is met. Then the next contour can be used to make a new band of triangles until all contours are used. The size of the triangles varies with the number of points selected per slice.
Figure IV.6: Triangulated surface of the head with a 3D-view of the
brain inside
An example of the triangulation of the head can be seen in figure IV.6.
To
present the results of any localisation, it would seem advantageous to use 3D
images of the head and brain. This allowed the visualization of a source
in its context, and should aid in finding the correct structure in which the
activity takes place.
Using ray-tracing or volume-rendering techniques (Robb et al, 1989, Tiede et al, 1990) it is possible to generate 3D-images from a segmented data set. For this purpose, the MRI data set is considered to be part of a larger virtual 3D-space. The object of interest, in our case the head, is considered to be illuminated by one or more imaginary light sources. The object is viewed on the computer screen. Light from an imaginary light source is reflected by the object to the eye of the observer.
Each pixel has to be illuminated in accordance with the light reflected from a 3D object. To be able to calculate the intensity of a pixel, a ray is directed from it into the virtual space until it hits the object of interest (figure IV.7). This can easily be established because of the segmentation described above. The computer ascertains how this point is illuminated by the light source. A method called Phong shading is used for this purpose. In this shading method, reflection of light is considered to consist of three elements. This is described by the following formula:
The first term, the ambient reflection, represents the light shining indirectly on the object. This term comprises the intensity of the ambient light Ia, and the reflection constant for ambient light Ka.
The
third term represents light that hits the surface of the object from the light
source directly, and is reflected in all directions (see figure IV.8). The angle
at which the light hits the surface is important for the amount of light that is
reflected. This is represented by taking the inner product of the surface
normal N and the unit vector in the direction from which the light comes L,
multiplied by a diffuse reflection constant Kd and the intensity of the light
source Is. If the inner product is negative, the diffuse term is assumed to be
zero, because no light can hit the surface at that point.
The second term represents the effect that most of the light is reflected at
the same angle as that at which it hits the surface. This is called
specular reflection.
The main term is the inner product between the unit vector in the viewing
direction V and that in the direction of the main reflection H. By raising this
term to the power of n, the reflection becomes more specular. If the
surface has to be shiny, n is chosen high and otherwise n is chosen low.
The main term is multiplied by a specular reflection contant Ks and the
intensity of the light source Is.
Hence, the directions of the light source L, the viewing direction V and the
normal vector N to the surface have to be known. N can best be calculated
from the original MRI data. The gradient of the values of the voxels at
the surface point of reflection is usually a good measure for the normal vector
on the surface. All the constants can be chosen more or less freely, and
together they determine the appearance of the image, and give an impression of
what the 'material' is made of. The values we use are as follows; n=10,
Is..Ks.=100,Is.Kd=150, Ia.Ka.=50.
Figure IV.9: 3D-View of the brain using integrative shading
In this way any object can be visualized. In figure IV.10 an example is
shown. Improvement can be made by taking shadows into consideration.
This involves tracing back from the surface point toward the light source to see
that nothing is blocking the light. However, this significantly enhances
the already high computation times, and does not seem to be of great importance
for displaying the head. Another improvement on the 3D effect is to
suggest motion by rapidly displaying different views of the head in
sequence. The best results are obtained by keeping the light source at the
same position, relative to the screen, although it is computationally less
expensive to keep the light source at a constant position with respect to the
head. In the latter case most of the calculations only have to be carried
out once for all views.
Although this rendering method works well for the head, it is less suitable for
generating images from the brain because it is difficult to obtain its exact
shape, including all the gyri and sulci. Looking back to our segmentation
method of the brain, it is so designed that it generates a more or less smooth
surfaced brain.
Using a different method called transparant rendering, it is still possible to
obtain a good 3D-image of the brain (Bomans et al., 1990). In this method,
a ray is traced from the screen towards the brain. When the ray hits the
surface, the original gray value is stored, and the ray proceeds into the
brain. In each step the gray value is added to the precious value, and
afterwards the average value is calculated. The values of voxels at a
gyrus are low due to the fluid, but on the edges of a gyrus the values are
high. Hence, averaging gives good 3D images. Bomans used an
integration of 6 voxels, but in our case this rendered low contrast images of
the brain with a lot of surface details such as veins. Since we are not
interested in the veins, we use the average of the squared values of nine voxels;
this works very well. This is demonstrated in figure IV.9.
Because most of the neuromagnetic sources are inside the brain, a 3D-display
of the outside of the brain may not be the most useful one. One of the
possibilities of obtaining a better display is to intersect the head by an
imaginary plane, and to show the original MRI data on the head-plane
intersection.
The easiest way is to start from a full 3D-display of the head which has been
generated earlier. The intersecting plans can be chosen in any direction
or location. The points where rays cast from every pixel on the screen
intersect the chosen plane are calculated. If a point lies inside the
head, the original voxel value has to be retrieved and displayed. If a
point does not lie within the head, it has to be determined whether it lies in
front of the head or behind the back, as seen from the viewing direction.
In the first case, the original information on the screen, a pixel of the
3D-display of the head, has to be retained. Otherwise the pixel has to be
made equal to the background, effectively erasing any part of the image that may
have been there. See figure IV.10 for an example.
Figure IV.10: 3D-View of the head, with a part cut away
(c) MEG, EEG and the integration with Magnetic Resonance Images, H.J. Wieringa, 1993
[<<Chapter III][Contents][Home][Chapter V: Integration of the MEG and EEG with the MRI>>]