- Email: [email protected]
Automated registration of multimodal brain image sets using computer vision methods
Computers in Biology and Medicine 29 (1999) 333±359
www.elsevier.com/locate/compbiomed
Automated registration of multimodal brain image sets using computer vision methods Gleb Secretta, Peter H. Gregson* Computer Vision and Image Processing Laboratory, Department of Electrical and Computer Engineering, DalTech, Dalhousie University, P.O. Box 1000, Halifax, NS, Canada B3J 2X4 Received 7 August 1998; accepted 13 May 1999
Abstract We present a new method of registering three dimensional image volumes of the brain of a given patient acquired at dierent times or with dierent imaging modalities, or both. Registration is an essential requirement for fusing the data from the two image sets so as to either increase the available information by exploiting complementary imaging modalities, or to measure small changes over time for prognostication, disease assessment, etc. The new technique exploits an external, removable, remountable reference frame which is attached to the head. Computer vision techniques are used to determine the positions of ®ducial marks in every image. The transformation required to map each image of one image volume onto the other image volume is developed using the theory of quaternions. The results indicate that the new technique is robust and practical in a clinical setting. # 1999 Elsevier Science Ltd. All rights reserved. Keywords: Image registration; Fiducial reference; MRI; SPECT; 3D imaging
1. Introduction Clinicians use comparative two- and three-dimensional imaging to quantify changes in tissue structure or function and to accurately localize functional information with respect to anatomical information. For tracking changes in structure such as the growth of some lesions, the clinician compares the apparent volume of the lesion as shown in MRI (magnetic * Corresponding author. Tel.: +1-902-494-6050; fax: +1-902-422-7535. E-mail address: [email protected] (P.H. Gregson) 0010-4825/99/$ - see front matter # 1999 Elsevier Science Ltd. All rights reserved. PII: S 0 0 1 0 - 4 8 2 5 ( 9 9 ) 0 0 0 1 3 - X
334
G. Secretta, P.H. Gregson / Computers in Biology and Medicine 29 (1999) 333±359
resonance imaging) image sets acquired at dierent times. In order for this comparison to yield clinically useful information early in the course of the disease, this comparison must result in accurate estimates of changes in lesion volume. Comparison of image sets is exceedingly dicult unless the two image sets are acquired with identical orientation and position. This is dicult to do in practice. SPECT (single photon emission computed tomography) imaging is used to detect areas of high metabolic rate such as tumors. The low spatial resolution of SPECT prevents accurate localization of the tumor site, however. An attractive solution is to use SPECT to acquire functional information and to use MRI to determine anatomical location. In this instance, the task is to present the functional information on the anatomical image, requiring accurate registration of dissimilar image sets. This task is especially demanding because dierent modalities usually have dierent resolutions as well as dierent coordinate system orientations. We present a new head-mounted reference frame and an algorithm for registration of two sets of brain images acquired with either dierent imaging modalities or at widely dierent times, or both. The new approach uses a new, light-weight, adjustable reference frame designed to ®t inside various imaging systems. Once adjusted for the particular patient, it can be repeatedly re-mounted on that patient with little positional error. Techniques from computer vision are used to detect the ®ducial marks in the images. A new algorithm which exploits the redundancy of the data is used for determining the vector-quaternion transformation which maps the coordinate system of one image set into the coordinate system of the other. Multilinear interpolation is applied to reduce quantization eects caused by dierent spatial resolutions of the data sets and by non-cubic voxels. The balance of this paper is organized as follows: Section 2 discusses the background to multi-modal image registration and describes SPECT and MRI imaging methods along with the restrictions imposed by each method. Section 3 presents the proposed approach to image registration. Experimental results and error analyses are given in Section 4. Conclusions and suggestions for future work are presented in Section 5.
2. Background In clinical practice and research, it is often necessary to determine both function and anatomy of brain structures of interest. Approaches for accomplishing this include functional MRI (fMRI) and `fusion' of SPECT and MRI image volumes. fMRI has the potential for providing both anatomical and functional information simultaneously. With fMRI, it is possible to determine the concentration of certain substances such as lactate that are the products of neurophysiological activity [1]. Detection of changes of lactate within speci®c brain regions in response to controlled stimulation should make it possible to use spectroscopic imaging to investigate stimulated brain function. There is much on-going research into fMRI but it is not currently a clinical tool. Further, the requirement for controlled stimulation makes fMRI an excellent research tool for establishing brain function but it is not suitable for localization of tumors and lesions. The most promising method of gaining both functional and anatomical information is based on `fusing' MRI and SPECT image sets [2]. However, the two modalities have dierent spatial
G. Secretta, P.H. Gregson / Computers in Biology and Medicine 29 (1999) 333±359
335
resolutions in all three directions and since the image sets are taken at dierent times with dierent imaging systems it is highly unlikely that the image planes of the two modalities would even be parallel. While some of the desired information can be determined from the two image sets without further processing, improved sensitivity to small lesions requires registering the two image sets by transforming the images of one set to match the three-dimensional spatial coordinates of the other set. In principle, three-dimensional image registration involves interpolating full three-dimensional information from the set to be transformed and then resampling it with the sampling parameters of the other set. We propose to incorporate these two steps into a single step so as to reduce computational burden and memory requirements. Applications for three-dimensional image registration include determining the anatomical locations of tumors found in SPECT imaging and comparing MRI image sets taken at signi®cant intervals for assessing disease progression. We present a short description of the two imaging modalities so as to establish the context for the proposed approach. 2.1. MR imaging MRI produces images based on the density of isolated and loosely-bound protons in the tissues of the body by measuring the intensity of a radio-frequency (RF) signal emitted by the tissue after an exciting magnetic ®eld and RF pulse are removed. This technique is based on the fundamental property that protons possess a magnetic moment and spin. When placed in a magnetic ®eld, the proton precesses about this ®eld. Initially all the protons are aligned either parallel or antiparallel to the magnetic ®eld. When a RF signal having an appropriate strength and frequency is applied to the object, the protons absorb energy, and more of them precess to the antiparallel state. When the applied RF signal is removed, the absorbed energy is reemitted and detected by an RF receiver. The proton density in the environment can be determined from the characteristics of the received signal. By suitably controlling the applied RF signal and the surrounding magnetic ®eld, and by appropriately processing the recovered RF signal, a three-dimensional image of the body can be formed as a set of slices. MRI is now widely used in clinical practice to image internal anatomy. Typical resolution of MRI machines is about 0.5 mm/pixel in the slice plane and about 5 to 7 mm between slices [3]. An example of an image acquired using MRI imaging is shown in Fig. 1. 2.2. SPECT imaging The goal of SPECT is to reconstruct, from projections, the three-dimensional distribution of radioactivity inside the patient's body. Dierent projections are acquired by counting the number of photons collected on a planar, collimated detector placed at various positions around the body. These detectors have point-spread functions (PSF) due to non-perfect collimation and to the attenuation of the photon's energy during its ¯ight from source to collimator. In addition, Compton scattering, gamma-ray attenuation and other physical phenomena make the data noisy thereby lowering SNR (signal to noise ratio). Various pre-reconstruction and post-reconstruction techniques have been used to enhance SPECT data such as constrained deconvolution, Wiener ®ltering and count-dependent Metz
336
G. Secretta, P.H. Gregson / Computers in Biology and Medicine 29 (1999) 333±359
®ltering [4]. Maximum entropy and homomorphic ®lters have also produced good results [5]. However, anatomical localization with SPECT is dicult because the method is designed to study physiological activity of the imaged organ and not its structure. In addition, SPECT images are often noisy and their spatial resolution ranges from 4 to 10 mm/pixel. SPECT data are characterized as being noisy, having poor resolution and showing function rather than anatomy [6]. An example of an image acquired using a SPECT scanner is shown in Fig. 1.
Fig. 1. Imaging with MRI and SPECT. Images a, b and c show MRI images. (a) Placement of the imaging planes. (b and c) Two adjacent slices. Images d, e and f show SPECT images. (d) Placement of the imaging planes. (e and f) Two adjacent slices. The dierence in spatial resolution is readily apparent.
G. Secretta, P.H. Gregson / Computers in Biology and Medicine 29 (1999) 333±359
Fig. 1 (continued)
337
338
G. Secretta, P.H. Gregson / Computers in Biology and Medicine 29 (1999) 333±359
2.3. The registration problem It is obvious that to solve the registration problem for sets of brain images, a transformation that maps one of the sets into the other must be found [7]. Since the brain doesn't appreciably deform or move inside the skull, we model it as rigid body. Thus, the transformation required is a simple translation±rotation combination and doesn't require anisotropic scaling. The desired transformation may be found using either anatomical landmarks in the two image sets or ®ducial marks based on an external reference [8]. Use of anatomical landmarks is not appropriate because the poor spatial resolution of SPECT imaging as seen in Fig. 1 can lead to a huge registration error. Even detection of the landmarks is highly problematic. Use of external ®ducial marks requires axing a reference frame external to the skull. The reference frame must be removable and re-axable in the same position so as to permit registration of images acquired at signi®cantly dierent times (perhaps months or years apart), as well as to avoid imaging-system scheduling problems (diculties with scheduling patients for both modalities on the same day). Ethical and practical considerations (availability and cost of operating suite time) prohibit the use of surgically-implanted markers. While this external frame is inconvenient for both the patient and the imaging process, it provides reliable information as to the orientation of the patient's head with respect to the imaging system.
3. The proposed approach The proposed registration procedure is: 1. 2. 3. 4.
Capture the two image sets with the external reference frame ®tted to the patient. Determine orientation of the head in each image set. ompute the transformation from one image set to the other. Use the transformation and interpolate voxels to determine values for the lower-resolution image set at the spatial locations of the higher-resolution image set.
The proposed reference system will appear on each slice in both sets as four bright spots outside the region of the patient's head (Fig. 9a). With the coordinates of the reference points in dierent slices it is possible to compute the orientation of the reference system with respect to the imaging system. Determining the transformation which maps one image set into another is straightforward with the aid of quaternion mathematics. Finally, bilinear interpolation is used to generate an image set derived from the lower-resolution set with the voxel coordinates (relative to the head) of the higher-resolution set. 3.1. Design of the external reference frame The external reference system consists of an adjustable frame and four detachable triangles. The head frame ®ts on the head and provides one keyed mounting location for a large triangle on each side of the head (Fig. 2). The frame is mounted on the head at the bridge of the nose and the auditory meati (ear canals). This provides a stable, reliable, repeatable mounting of the reference system on every
G. Secretta, P.H. Gregson / Computers in Biology and Medicine 29 (1999) 333±359
339
Fig. 2. (a) The head frame with triangles. Triangles are mounted in the sagittal plane. The head frame is adjustable and may be removed from the patient without disturbing adjustments. See text for explanation. (b) Detail of one triangle showing the mounting of the tube.
340
G. Secretta, P.H. Gregson / Computers in Biology and Medicine 29 (1999) 333±359
skull when properly adjusted. MRI imaging requires that no ferromagnetic parts be employed in the reference system and so the head frame and triangles are made of Plexiglas. Positioning of the head frame is provided by ®ve points of contact with the patient, two of which are adjustable. Non-adjustable points provide support of the reference system at the bridge of the nose (A) and two points on the upper part of the skull (B and C). Adjustable points D and E are `earplugs', each having three degrees of freedom to ®t into the auditory meati of the patient. After points D and E are properly adjusted for the patient, the reference system can be removed by undoing bolt F without disturbing the adjustment of points D and E. The frame may be remounted by ®tting it onto patient's head and retightening bolt F. Around the perimeter of each triangle is cut a narrow groove into which a ¯exible tube is mounted. This tube is ®lled with either a radio-emitting source to make its lumen easily detected in SPECT imaging or cod liver oil which renders its lumen easily detected in MRI. The triangles are otherwise identical and are made easily removable and replaceable so as to permit the appropriate pair of triangles to be used for each imaging modality. Accurate triangle positioning is assured because each triangle is keyed to the head frame. By making triangles removable, radiation exposure to patients and sta is minimized. The triangles are mounted in the sagittal plane as shown in Fig. 2. This is a matter of convenience for use with head coils and for the comfort of the patient. It also permits the use of two large triangles which maximizes the accuracy of head orientation measurements due to both triangle size and triangle redundancy. For both coronal and trans-axial imaging, sagittal placement of the triangles will result in ®ducial marks in most images. This placement of the triangles is not optimal for sagittal imaging, however, because very few images will contain ®ducial marks. A simple modi®cation to the current design of the head frame would permit sagittal imaging. The head frame is adjustable to ®t each patient. While our prototype does not have scales on each adjustment, scales could easily be provided on a production version of the device. This would make it possible to accurately re®t the head frame to patients on repeat visits for acquisition of sequential images.
3.2. Determining head orientation Consider the situation depicted in Fig. 3. In Fig. 3, one imaging plane is shown intersecting with both triangles. The points of intersection are shown as dots. Other imaging planes will intersect the triangles at dierent points, with the result that each side of each triangle will be described by a series of dots. By ®tting lines to the dots, the locations of the sides of each triangle can be estimated. Since the triangles are fastened to the patient's head by the head frame, head position and orientation can be determined. Clearly the method requires only one triangle; the second triangle provides redundancy to improve accuracy although in the current work it was used to check the results. The method depends on accurate determination of the intersection points of the reference triangles with the image slice planes; speci®cally, we must accurately detect the centers of the dots appearing in each image slice so as to accurately estimate the locations of the triangle sides.
G. Secretta, P.H. Gregson / Computers in Biology and Medicine 29 (1999) 333±359
341
Fig. 3. The intersections of the imaging plane with the perimeter tubes of the two reference triangles are indicated by the four large dots.
3.2.1. Finding dot centers Intersection of one image slice plane with the sides of the two triangles results in four dots appearing outside the head (Fig. 9a). To detect these dots, the operator selects either the left or the right side of the image set to provide the reference. This is done once per image set and is the only operator intervention required. Data from the other side are currently used for comparison, but may be incorporated into orientation determination in future. The dots are then isolated from the remainder of the image data. Since the reference system is designed so that the dots will appear about 30 mm from the head, a simple masking operation is used. The original image is low-pass-®ltered to achieve two goals. It must remove high-frequency noise and blur the image to `attach' those parts of the image that appear disconnected from the head (the ears or the nose can produce disjoint segments). The blurred image is thresholded at 10% of maximum intensity to produce a binary image. In this image, the head should be a single segment and the four dots are readily visible, disjoint from the head as required. Finally, area thresholding is employed to discard segments with areas either too small or too large, resulting in a binary `mask' image. This mask image is ANDed with the original image to produce a grey-level image in which only the dots appear. The lower and upper area thresholds (determined empirically) for MRI are 2 and 100 pixels and for SPECT are 5 and 50 pixels. Tubing that is not perpendicular to the plane of an image slice appears not as a disc but rather as a ®gure consisting of two elliptic and two straight segments due to the thickness of the slice (Fig. 4). Due to symmetry considerations, we assume that the center of a reference dot
342
G. Secretta, P.H. Gregson / Computers in Biology and Medicine 29 (1999) 333±359
represents the point that lies closest to the mid-line of the tubing, and thus can be used for determining the orientation of the corresponding triangle side. We have found that this procedure yields acceptable results, however we note that morphological processing and/or various `dot detectors' might yield better performance. 3.2.2. Validating dot images Since there can still be noise in the grey-level dot image, the algorithm must discriminate between valid and invalid candidate reference dots. We threshold the dot image at 30% of its maximum intensity. This results in suppression of some noise and makes the dots smaller. Resultant binary dot images are then analyzed further for the shape, number and relative position of dots. The program detects all the dots in the image and assigns each dot a unique label using a `connected components' algorithm. Coordinates of dot centroids are found for each dot as coordinates of its center of gravity. Binary dot images are used for this purpose, since the grey-level variation in the original image throughout the dot interior is very small. An acceptable dot image must meet the following criteria: 1. There must be exactly two dots in the reference half of the image. 2. The distance between the centroids of the putative reference dots must exceed 5 mm for MRI and 15 mm for SPECT images. 3. Areas of valid reference dots must be at least 3 pixels and not greater than 15 pixels for both image sets. 4. The standard deviations, measured in pixels in both the X and Y directions about each dot centroid, must not exceed 5 pixels. The ®rst three criteria were devised to detect various conditions that can appear in marginal slices, those slices which either do not contain dots because they do not intersect the reference triangles, or which contain images of a reference triangle near a vertex. Criterion 4 handles the
Fig. 4. View of the tubing crossing a slice. The tubing will appear as an ellipse due to the angle made by the tubing centerline with the slice and the thickness of the slice.
G. Secretta, P.H. Gregson / Computers in Biology and Medicine 29 (1999) 333±359
343
situation in which the imaging slice plane intersects the tubing wrapped around a triangle vertex, resulting in elongated dots which are unreliable for dot centroid calculations. 3.2.3. Estimating triangle sides The set of points due to the detection of the sides of one triangle in adjacent image slices constitutes a cloud of points. To estimate the location and orientation of this side, we ®t a straight line to this cloud. We assume a coordinate frame in which the X and Y axes lie in the image slice plane and Z is orthogonal to it. Physical constraints imposed by the possible orientations of the patient's head within the imaging system, the placement of the reference triangles on the head frame and the positioning of the head frame on the patient cause two sides of the reference triangle to meet the imaging plane at a large angle. This allows us to determine the orientation of the sides of each reference triangle by ®tting lines to the projections of its cloud of points on the Z±X and Z±Y planes. The parameters of these lines provide the information that we require. We use orthogonal regression to ®t lines to the projections of the cloud. With the general equation of the line in the 2D Z±X space in which z is the independent variable: az bx c 0
1
Line ®tting to a set of points by orthogonal regression in 2D assumes that each point is aected by independent identically distributed noise in both x and z directions. In that case the analysis tells us that: 1. The solution line must intersect the centroid of the points xi, zi, i = 1, . . ., n. The centroid coordinates xc, zc are x c , zc
N 1 X x i , zi N i1
2
P 2. The normal vector is the `smallest' eigenvector of the covariance matrix xz of the points, i.e. the eigenvector corresponding to the smallest eigenvalue. The covariance matrix is de®ned as X xz
N X x i ÿ x c 2 x i ÿ x c zi ÿ zc x i ÿ x c zi ÿ zc zi ÿ zc 2 i0
The covariance can be rewritten as X p q , q r xz P Then the smallest eigenvalue of xz is q p r ÿ p ÿ r2 4q2 l 2
3
4
5
344
G. Secretta, P.H. Gregson / Computers in Biology and Medicine 29 (1999) 333±359
The corresponding eigenvector is q lÿp Then the line ax+bz+c = 0 that has normal vector (q l ÿ p )T and intersects the point (xc, zc) is given by the parameters aq blÿp c ÿ ax c bzc :
6
3.2.4. Finding the reference vertex Only two sides of each triangle are used for determination of head orientation. Accuracy and precision of the orientation measurement requires that each triangle side is as close as possible to being parallel to the surface normal of the imaging plane, and further, that the sides meet at an angle that is not too acute. These con¯icting requirements and operational requirements on triangle size are satis®ed by only one vertex of each triangle, hereafter referred to as the reference vertex. In order to compute the translational component of the desired transformation, it is necessary to determine the coordinates of the reference vertex. This vertex should be at the point of intersection of the two lines which are estimates of the triangle sides. However, due to imaging errors, noise, and errors in locating dot centers, the two lines will rarely intersect in three-dimensional space. Thus, we take as the reference vertex the mid-point of the line orthogonal to both line estimates. 3.3. Computing the transformation The transformation desired consists of a translational component and a rotational component. The translational component is easily computed from the coordinates of one triangle vertex in each of the image sets. We convert pixels and voxels to millimeters in both image sets so as to use real units and also to avoid diculties with conversion of measurements. In our case we are dealing with a three-dimensional problem. Position of a rigid body in 3D space is determined by the position of two non-parallel lines belonging to the body, hence it is sucient to use only one triangle for the reference. Since the brain is assumed to be a rigid body for this research, we need not be concerned with anisotropic scaling and so we use unit vectors with the orientations of the lines. To determine the rotational part of the transformation, we need to ®nd two corresponding pairs of normalized, non-parallel vectors and the combination of rotations that transforms one pair into the other. The Euler theorem states that any two states of a rigid body with a single ®xed point are
G. Secretta, P.H. Gregson / Computers in Biology and Medicine 29 (1999) 333±359
345
separated by only one rotation about some axis. To ®nd this axis and the angle of the rotation we use the algebra of quaternions because it ®ts this class of problems extremely well [10]. Any quaternion q can be written in the following form: q q0 q1 i q2 j q3 k qs q~v ,
7
where i, j and k are the triple de®ning the Cartesian coordinate system, qs=q0 is the scalar part of the quaternion q and q~v =q1i+q2j+q3k is the vector part of the quaternion q. The norm of a quaternion q, written Nq is de®ned as q 8 Nq q20 q21 q22 q23 Thus, any quaternion q can also be written as: q Nq cos Q e~ sin Q,
9
where cos Q q0 =Nq 1 q q2 q22 q23 sin Q Nq 1 q1 i q2 j q3 k e~ 2 q q21 q11 q23 The normalized quaternion (q0=1) actually represents conical rotation of any vector a~ about an axis q~v passing through the origin of a~ for the angle of 6 q~v6 (Fig. 5). The fundamental theorem is given in Ref. [9]: Theorem 1. If q and r are any non-scalar quaternions, then
Fig. 5. Rotation in terms of quaternions.
G. Secretta, P.H. Gregson / Computers in Biology and Medicine 29 (1999) 333±359
346
r 0 qrqÿ1
10
is a quaternion whose norm and scalar are the same as for r. The vector V~r 0 is obtained by revolving V~r conically about V~q through twice the angle of q with respect to r (Fig. 5). Thus, if q Nq cos Q e~ sin Q,
11
then V~r 0 is obtained by revolving V~r conically about e~ through an angle 2Q. The normalized quaternion has Nq=1. Eq. (11) for a normalized quaternion is simpli®ed to q cos Q e~ sin Q:
12
Consider now the special case when vectors V~r and V~r 0 lie in the plane of rotation. Knowing vectors V~r and V~r 0 allows the angle between them and the orientation of vector e~ to be
Fig. 6. Finding the rotational part of the registration transformation as a combination of two rotations. (a) First rotation. (b) Second rotation.
G. Secretta, P.H. Gregson / Computers in Biology and Medicine 29 (1999) 333±359
347
computed. Vector e~ is simply their vector product. Thus, a quaternion can be de®ned using a pair of vectors. Now, a succession of rotations must be de®ned in terms of quaternions. Since the product of quaternions is also a quaternion then any number of rotations is equivalent to a single rotation. The succession of two rotations de®ned by quaternions p and q, respectively, corresponds to the operator r 0 qp rpÿ1 qÿ1 qp r qpÿ1
13
Thus, we can ®nd the quaternion which corresponds to the rotational part of the transformation that transforms one pair of vectors in one image set into the corresponding pair in the other image set. We de®ne our vectors in each set as vectors with origin at the reference vertex and with the orientation of each triangle side adjacent to that vertex. The estimates of the sides of the triangle are used to provide the pair of vectors in each image set. The problem is depicted in Fig. 6. There are two pairs of vectors V~1 , V~1 0 , V~2 and V~2 0 . Vector V~1 0 corresponds to V~1 , and vector V~2 0 corresponds to V~2 . The task at hand is to ®nd a normalized quaternion given two pairs of corresponding vectors. The solution is sought in the form of two successive rotations. Quaternions representing rotations can only be computed if corresponding vectors have the same length. To satisfy this requirement the longest vector is sought and the three other vectors are scaled up to have the same length. The ®rst rotation maps V~1 into V~1 0 . The axis of this rotation coincides with vector product V~1 V~1 0 . According to Eq. (12) the normalized quaternion that represents this rotation can be written as: q1 cos b
V~1 V~1 0 sin b jV~1 V~1 0 j
v !2 u V~1 V~1 V~1 V~1 u V~1 V~1 0 t 1ÿ q1 jV~1 jjV~1 0 j jV~1 V~1 0 j jV~1 jjV~1 0 j 0
0
14
15
The ®rst rotation results in V~2 going to V~2int . Now the second unit quaternion q2 that conically rotates V~2int into V~2 0 about axis V~1 0 can be calculated. As for the ®rst quaternion, two corresponding vectors that lie in the plane of rotation are used. The vectors P~2int and P~2 0 that are projections of V~2int and V~2 0 on the plane of rotation are used. Quaternion q2 that conically rotates V~2int into V~2 0 also rotates P~2int into P~2 0 and can be determined from them. Projection vectors P~2int and P~2 0 are found as V~2int V~1 0 ~ 0 P~2int V~2int ÿ V1 jV~1 0 j2
16
348
G. Secretta, P.H. Gregson / Computers in Biology and Medicine 29 (1999) 333±359
V~2 0 V~1 0 ~ 0 P~2 0 V~2 0 ÿ V1 jV~1 0 j2 Unit quaternion q2 is found according to Eq. (12) as v !2 u u 0 0 0 0 V~1 P~2 P~2int V~1 t P~2 P~2int sin g 1ÿ q2 cos g jV~1 0 j jP~2 0 jjP~2int j jV~1 0 j jPV~2 0 jjP~2int 0 j
17
18
Finally, according to Eq. (13), multiplication of quaternion q2 by quaternion q1 results in a quaternion q that represents rotation of vectors V~1 and V~2 into V~1 0 and V~2 0 respectively. Quaternion q is also a unit quaternion because it is the product of unit quaternions q1 and q2. Once this resultant quaternion q is found the rotation is de®ned. 3.4. Applying the registration transformation Since the transformation is carried out on a voxel-by-voxel basis, the low-resolution image is transformed to match the high-resolution image so as to retain as much image information as possible [11]. Thus, we transform the SPECT image set to match the MRI image set. The transformation is carried out by forming a vector from the reference vertex to the voxel under consideration in the MRI image set. If this vector is V~MRI , then by applying V~SPECT qV~MRI qÿ1 t~
19
where q is the quaternion representing the rotational part, t~ is the translational component of the transformation and vector V~SPECT in the SPECT image refers to the same anatomical location as V~MRI . The value found at this voxel in the SPECT image set is plotted at the MRI coordinates of the voxel under consideration in the transformed SPECT volume. Thus, the transformed SPECT image resulting from this procedure will have the SPECT data plotted with respect to the MRI coordinate axes, with the same anatomical points appearing at the same coordinates in each image set. The image sets are now in register. Two problems remain. 1. SPECT voxels are much larger than MRI voxels and so adjacent voxels in the transformed SPECT image have identical values, resulting in a strongly `pixelated' image. 2. Since the imaging planes of the SPECT data are not parallel to those of the MRI data, the transformed SPECT data exhibits visible edges due to data being taken from dierent planes. Both of these problems are evident in Fig. 8. To solve these problems, we use a multi-linear interpolation. Collections of new data points which all arise due to the same voxel in the original SPECT data are termed `supervoxels'. The value of the transformed voxel is computed as follows (Fig. 7): 1. Assume that the x and y coordinate axes lie in the slice plane and z is orthogonal to the slice plane. 2. Label the supervoxel containing the new voxel as K and the two nearest supervoxels on the
G. Secretta, P.H. Gregson / Computers in Biology and Medicine 29 (1999) 333±359
349
Fig. 7. Computing the grey scale value of a new voxel. (a) Labeling of the supervoxels. (b) Interpolation in the plane of the new voxel.
same plane as L and M. 3. Label supervoxels on the nearest plane to the new voxel and adjacent to K, L and M as K', L' and M ', respectively. 4. Using the x and y coordinates of the new transformed voxel, 4.1. Compute G as the linear interpolation of K and the average of L and M. 4.2. Compute G' as the linear interpolation of K'and the average of L' and M '. 5. Compute the new value of the transformed voxel as the linear interpolant of G and G'. This interpolation solves both pixelation and interplane transitions eects (Fig. 8).
350
G. Secretta, P.H. Gregson / Computers in Biology and Medicine 29 (1999) 333±359
Fig. 8. Results of interpolation. (a) Transformed SPECT image without interpolation illustrating `pixelation' and plane transitions. (b) Transformed SPECT image after multi-linear interpolation.
G. Secretta, P.H. Gregson / Computers in Biology and Medicine 29 (1999) 333±359
351
4. Experimental results 4.1. Head frame mounting accuracy Experiments were conducted to determine the registration accuracy of the new technique. Accurate registration requires that both the rotational and translational components of the transformation be correct and that the voxel values of the transformed image set are correctly computed. The ®rst step was to determine the accuracy with which the head frame could be repositioned. To perform this experiment, a volunteer's head was placed in a known position with the aid of a custom-®tted bite plate molded to ®t the volunteer's teeth. The bite plate was axed to a Plexiglas bar (5 cm wide, 5 mm thick, 15 cm long) which was clamped in a vise that was bolted to a table. The head frame was equipped with a forward-projecting laser. The volunteer positioned his head accurately by biting on the bite plate, moving his head until his teeth indexed positively on the bite plate. An assistant mounted the head frame on the volunteer, and the position of a small spot created by the laser on a tangent screen placed 1.53 m in front of the laser was recorded. The head frame was removed and remounted 30 times, with the location of the laser spot being recorded each time. After acquiring the data, the mean and standard deviation of the dot coordinates were computed. These values were computed only for the horizontal plane (motions due to trans-axial rotations) for two reasons. Firstly, the geometry of the head frame suggests, and the volunteer con®rms, that re-positioning accuracy should be lower with respect to transaxial misalignment than with respect to misalignments in the other two planes. Users report that the head frame mounting seems less ®rmly anchored to the head for coronal and transaxial rotations than for sagittal rotations. Secondly, the volunteer found it extremely dicult and very time consuming to bite down on the bite plate without causing its Plexiglas mounting to bend. For lateral (transaxial) rotations, the Plexiglas bar resisted bending very well as it was 5 cm wide. Further, lateral mis-alignments of the head were very obvious to the volunteer as they lead to discomfort in the jaws due to the teeth not ®tting the bite plate properly. In contrast, vertical mis-alignments resulted in signi®cant bending of the Plexiglas bar with relatively low forces applied because the bar was only 5 mm thick. Additionally, vertical motions did not result in loss of the ®t between teeth and the bite plate and so the volunteer had very little indication of vertical misalignment. Coronal rotations also lead to discomfort as under these conditions the Plexiglas bar is being twisted and the jaws feel as if they are being pried open. The results of the experiment are found in Table 1. The right-most column refers to the error distance measured on a cylinder of 200 mm radius, suciently large to enclose a head. Table 1 Results of the head frame remounting experiment
Standard deviation Maximum mag.
Distance at 1.53 m (mm)
Angular error (min, s)
Distance at 200 mm
2.41 5.97
5 '240 13 '230
0.31 0.78
352
G. Secretta, P.H. Gregson / Computers in Biology and Medicine 29 (1999) 333±359
We note that the maximum error magnitude is in the order of one voxel at this radius. Misalignment of the head frame will result in, at most, a one-voxel shift in the placement of the ®ducial reference marks in the images and, most of the time, a shift considerably less than this. Thus, the portion of the error due to head frame mis-alignment is generally greatly below one voxel. (Note that brain structures are more central and so the errors in alignment will be proportionally smaller.) We conducted MRI-to-MRI, SPECT-to-MRI and SPECT-to-SPECT studies to determine the accuracy of transformation. For assessing MRI-to-MRI registration, the reference triangles mounted on the head frame were ®lled with cod liver oil to generate bright spots in MRI images. To permit accurate re-positioning, we used a phantom consisting of a U-shaped piece of ABS pipe of 2-inch diameter, ®lled with water and rigidly attached to the head frame. This phantom shape allowed us to position the head frame precisely with the aid of a laser, while simultaneously providing and image-able object with sucient spatial complexity that we could verify our test data. For SPECT-to-SPECT and SPECT-to-MRI comparisons, we used a more physiologicallyrelevant phantom consisting of a plastic canister reinforced with plastic rods and containing ®llable compartments to emulate structures in the brain. The main part of the canister was ®lled with radioactive technetium (Tc). A phantom was used because ethical considerations precluded the use of a volunteer in these studies. For the SPECT-related studies, the tubes of the two reference triangles were ®lled with a radioactive solution with a higher concentration than the solution in the canister. Three SPECT data sets at diering orientations were acquired. Finally, two MRI image sets of a volunteer were acquired to illustrate algorithm performance with real images. 4.2. Estimating accuracy and precision To assess accuracy, we need two sets of images of one object in dierent poses. We call these set 1 and set 2. We wish to evaluate the dierences between set 1 and set 2 after set 2 has been transformed to match the pose of set 1. Ideally, there would be no dierence. The ®rst study assessed the accuracy of MRI-to-MRI registration. Our MRI scanner produced in-plane resolution of 0.87 mm/pixel and interslice spacing of 5.2 mm. For this study, the phantom was mounted on a test jig that permitted only out-of-image-plane rotation (the worst case) of the phantom. V-shaped notches were made in the test jig to denote angles to which the phantom with the attached headframe was to be rotated. After imaging in the MRI scanner, the actual angles were determined by ®tting a laser to the head frame to provide an optical lever, reproducing the angular positions of the head frame, noting the positions of the resulting laser spots on a tangent screen and computing the resulting angles. Repeated repositioning trials produced a positioning error of 21 '300. Four MRI scans of the phantom in four dierent angular positions were performed. Rotations were computed using each of the two triangles for each of the six pair-wise combinations of the four images (Table 2). From the data, it is clear that the rotations estimated by the two triangles are in close agreement with actual rotation angles and with each other. In Table 2, MR1 refers to MRI image 1. The columns headed `left triangle' and `right
G. Secretta, P.H. Gregson / Computers in Biology and Medicine 29 (1999) 333±359
353
Table 2 Angles of rotation in degrees for MR-to-MR registration Registration
Left triangle
Right triangle
Actual angle
Maximum error
MR1±MR2 MR1±MR3 MR1±MR4 MR2±MR3 MR2±MR4 MR3±MR4
11.81 8.59 11.25 7.42 19.46 12.47
11.78 8.38 11.39 7.37 19.10 11.98
11.75 8.22 11.30 6.98 19.72 12.93
ÿ 0.06 ÿ 0.37 ÿ 0.09 ÿ 0.44 0.62 0.95
triangle' present the angular dierence between the poses of the object (the phantom) in the two images as computed using the left triangle and the right triangle, respectively. The actual rotation angle and the maximum error are also given. The small rotation angles used represent physiologically-relevant positioning dierences between imaging sessions and so are realistic. The worst-case error on the surface of a cylinder with a radius of 200 mm is thus 3.3 mm or 3.8 pixels. The absolute average of the worst-case errors was 0.428 which is 1.47 mm or 1.69 pixels on the cylinder. The angular standard deviation of all estimates was 0.448 which is 1.52 mm or 1.7 pixels on the 200 mm cylinder. While this is not strictly the precision of the system (it does not look at the `spread' of the errors for all pixels within one pair of image sets but rather the spread for six imageset pairs), it does provide some indication of repeatability. These results are clinically acceptable as changes that might be overlooked due to errors of this size are not considered clinically signi®cant. The second study looked at SPECT-to-MRI registration using the canister phantom. One MRI image set and three SPECT image sets were acquired. The SPECT system captures 64 64 pixel images, with pixels about 4 mm square. The consistency with which the SPECT image sets were brought into register with the MRI set are shown in Table 3. We were not able to determine accuracy because to do so we need to accurately determine the inherent coordinate systems of both scanners. The coordinate system of the MRI scanner may be modi®ed by changes to pulse sequences and it is also aected by minor disturbances in its magnetic ®eld and so the reliability with which we could determine its coordinate frame was considered dubious.
Table 3 Angles of rotation for the SPECT-to-MR registration Registration
Left triangle
Right triangle
Angle error
SP1±MR SP2±MR SP3±MR
5.2 5.3 11.3
5.4 5.6 11.9
0.2 0.3 0.6
354
G. Secretta, P.H. Gregson / Computers in Biology and Medicine 29 (1999) 333±359
Table 4 Angles of rotation for the SPECT-to-SPECT registration Registration
Left triangle
Right triangle
Angle error
SP1±SP2 SP1±SP3 SP2±SP3
7.8 11.6 7.2
7.5 11.1 7.9
0.3 0.5 0.7
Comparison of the transformation angles estimated with right and left triangles shows that they agree well. This shows that the estimates are consistent. For this study, the worst case discrepancy is 0.68 which is 2.1 mm or about 0.5 pixels on the surface of the 200 mm cylinder. The consistency data for the SPECT-to-SPECT registration study appear in Table 4. The worst-case dierence in angular position is 0.78. This corresponds to 2.4 mm on the 200 mm cylinder or 0.6 pixels. This is clinically acceptable precision. Finally, real data were analyzed to illustrate the eectiveness of registration. Two MRI image sets of a volunteer were acquired and registered. Real data are useful because they have higher spatial frequencies than our phantom and so any mis-registration is clearly visible. Fig. 9(a) shows one image from the ®rst image set of the pair. Fig. 9(b) shows the corresponding plane, after the registration transformation, taken from the second image set of the pair. The algorithm's ability to retain relatively high spatial frequencies is clear. The sulci between the gyri are retained to a large degree and the distinction between grey matter and white matter is also clear. There is some loss of resolution; this is unavoidable if we are to prevent `pixelation' and inter-plane discontinuities in the image slices. The small black wedge at the bottom of the image represents the amount in the x±y plane by which the registration algorithm was required to transform the second image set. Not shown is the out-of-plane rotation of the head for which compensation was also accomplished. Fig. 9(c) shows the dierence or error image. The grey levels in this image are not meaningful as our display stretches contrast. The degree of mis-registration is visible in Fig. 9(c), appearing to be a sort of `double exposure' of the head. Note that should registration be perfect, the image would be a uniform grey level. While not perfect, the registration error is very small, certainly compared to the error without registration. Typical maximum deviation from the mean grey level inside the region of interest (the brain) in the error image is 230 grey levels, or 212% of the whole range of 256 grey levels, while the standard deviation for the same region equals 5 grey levels or 2%.
5. Conclusions and future work A new reference frame and a new algorithm have been developed to enable rapid, clinicallyuseful registration of image volumes both between imaging modalities and within a given modality for sequential scanning. Applications of the these developments include integrating
G. Secretta, P.H. Gregson / Computers in Biology and Medicine 29 (1999) 333±359
Fig. 9. MR-to-MR registration. (a) Benchmark image. (b) Registered image. (c) Result of subtraction.
355
356
G. Secretta, P.H. Gregson / Computers in Biology and Medicine 29 (1999) 333±359
Fig. 9 (continued)
functional and anatomical data, monitoring disease progression and planning invasive therapies. The new approach has proved to be an accurate tool for registration of unimodal and multimodal three-dimensional data sets. The head frame provides sucient adjustment range while preserving an acceptable degree of positional repeatability when remounted. Performance of the system, consisting of the head frame and the algorithm, is sucient for many applications where 3D registration is required. MRI-to-MRI registration error was under 2 mm for the region of interest (the head, as exempli®ed by the 200 mm cylinder). These errors are not considered clinically meaningful. The maximum discrepancy between estimates of angular transformation for the SPECT-toSPECT registration was 2.4 mm for our SPECT scanner. SPECT-to-MRI registration exhibited similar performance. Among the approaches to the registration problem that are reported in current literature there are two that reach similar accuracy. One of them proposes use of a combination of ®ducial markers and structured light to achieve registration [11] and has RMS error 1.1720.45 mm. This method implies use of certain auxilliary equipment such as special light projector and videocamera to capture the pattern of light on the patient's body. The other method does not require any additional equipment [12]. The algorithm minimizes a mutual information metric in an iterative way and provides suboxel accuracy. One drawback of this approach is the considerable amount of processing time needed for the algorithm to converge, which typically
G. Secretta, P.H. Gregson / Computers in Biology and Medicine 29 (1999) 333±359
357
is between 30 and 60 min. The other problem with this method is that it requires that images of both modalities have high resolution. Our algorithm uses minimal auxilliary equipment and is independent of the image resolutions in the data sets. However, further work is required to more accurately determine the angle of rotation, as we suspect that much of our error derives from noise and mislocalization of reference dots in cases when most of them fall close to the vertex of the triangle. One way to overcome this drawback is to increase the dimensions of the triangles and mount them closer to the patient's head. The system should also be modi®ed to permit image acquisition in all three planes. Finally, the transformation should be estimated using both triangles instead of only one. This would provide more measurements. It would also permit making measurements over greater displacements (between triangles) thereby increasing accuracy.
6. Summary We present a new method of registering three dimensional image volumes of the brain of a given patient acquired at dierent times or with dierent imaging modalities (MRI and SPECT), or both. Registration is an essential requirement for fusing the data from the two image sets so as to either increase the available information by exploiting complementary imaging modalities, or to measure small changes over time for prognostication, disease assessment, etc. The new technique exploits an external, removable, re-mountable reference frame which is attached to the head. The head frame is suciently adjustable to ®t most patients, the positional accuracy with which it may be remounted indicates that it contributes negligibly to the registration errors. The head frame is equipped with two removable triangles which have perimeter tubes containing contrast material (MRI) or radioactive material (SPECT). These triangles generate four dots in each image which intersects them. The position and spacing of these dots are dependent on the orientation of the head frame, and therefore the head, with respect to the imaging plane. These dots form the ®ducial references. Using computer vision techniques, the dots are automatically located in each image. The locations of the triangle sides are estimated from the set of dots displayed in all the images of the image set. The location of a reference vertex of the triangle is determined. This procedure is followed for both image sets. Quaternion mathematics is used to determine the rotational transformation which carries one head pose into the other. This transformation is then applied to each pixel in the image set. Bilinear interpolation is used to remove `pixelation' eects due to non-parallel imaging planes in the two image sets and to dierent spatial resolutions of imaging modalities. The new algorithm is tested with MRI-to-MRI, SPECT-to-MRI, and SPECT-to-SPECT registration tasks. Performance is deemed to be clinically acceptable, with registration errors below 2 mm throughout the head volume. SPECT-based registration is consistent to within 6 mm.
358
G. Secretta, P.H. Gregson / Computers in Biology and Medicine 29 (1999) 333±359
Finally, MRI-to-MRI registration was performed on the images of a volunteer. These results show that the method preserves much of the anatomical information present in these images. References [1] M. Singh, Towards proton MR spectroscopic imaging of stimulated brain function, IEEE Trans. Nucl. Sci. 39 (4) (1992) 1161. [2] G.K. Matsopoulos, S. Marshall, J.N.H. Brunt, Multiresolution morphological fusion of MR and CT images of the human brain, IEEE Proc.-Vis. Image Signal Process. 141 (3) (1994) 137±143. [3] R. Sauter, M. Shneider, K. Wicklow, H. Kolem, Current status of clinically relevant techniques in magnetic resonance spectroscopy, Siemens Electro-Medica 60 (29-56) (1992) 32. [4] D. Boulfelfel, R.M. Rangayan, L.J. Hann, R. Kloiber, Prereconstruction restoration of myocardial single photon emission computed tomography images, IEEE Trans. Med. Imaging 11 (3) (1992) 336±399. [5] W.A. Kalender, Quo Vadis CT? CT in the year 2000, Siemens Electro-Medica 61 (2) (1993) 30. [6] O. Kuler, G. Gerig, Computer-analysis of 3D medical images, Siemens Rev. (1991) 4. [7] P.J. Besl, N.D. McKay, A method for registration of 3D shapes, IEEE Trans. Pattern Anal. Mach. Intell. 14 (2) (1994) 239±256. [8] W.R. Fright, A.D. Linney, Registration of 3D head surfaces using multiple landmarks, IEEE Trans. Med. Imaging 12 (3) (1993) 515. [9] L. Brand, Vector and Tensor Analysis, in: John Wiley and Sons, NY, 1948, pp. 403±421. [10] R.J. Moorhead II, Z. Zhu, Signal processing aspects of scienti®c visualization, IEEE Signal Process. Mag. 12 (5) (1995) 20. [11] R.B. Backman, A temporal 3D-registration framework for computer-integrated surgery, Ph.D. thesis, Department of Computer Science, Robotics and Vision Group, University of Western Australia, 1998. [12] C.R. Meyer, J.L. BoeCR, B. Kim1, P.H. Bland, K.R. Zasadny, P.V. Kison, K. Koral, K.A. Frey, R.L. Wahl, Demonstration of accuracy and clinical versatility of mutual information for automatic multimodality image fusion using ane and thin-plate spline warped geometric deformations, Medical Image Analysis 1(3) (1996± 1997) 195±206. Dr. Peter Gregson (B.Eng. 74, M.Eng. 77, Ph.D. 88) is a Professor in the Department of Electrical Engineering at DalTech and is Director of the Computer Vision and Image Processing Laboratory. He has founded three companies designing and manufacturing data acquisition and control systems and was employed by the Defense Research Establishment. Dr. Gregson's research is focused on the development of new theory, architectures and algorithms for `early' computer vision and image processing and on applications of this technology to medical imaging, robot navigation and industrial control and measurement. He is developing theory, algorithms and analog VLSI circuitry for real-time, vision based short-range robot navigation, for automated surveillance monitoring and for motion estimation. Current medical imaging projects include automated assessment of diabetic retinopathy; automated determination of refractive error in pre-verbal children; automated correction of motion artifact in CT images; detection and quanti®cation of prostrate tumors and benign prostate hyperplasia in 3D ultrasound; multi-modal image registration and detection, localization and quanti®cation of multiple sclerosis lesions in 3D MRI. Dr. Gregson's laboratory has several ongoing industrial collaborations. Dr. Gregson has a major interest in eective, relevant teaching. He was the recipient of the Wighton Fellowship in 1996. This annual award is made to one professor of Engineering and Applied Science in Canada for innovative undergraduate laboratory presentation.
Mr. Gleb Secretta is a Ph.D. student in the department of Electrical Engineering at DalTech. He is a member of the Computer Vision and Image Processing Laborator. Mr. Secretta was born in 1967 in Moscow, Russia. He graduated with Bachelor of Engineering and Master of Engineering degrees in Systems and Control from the Moscow State Technical University in 1991. His Russian degree thesis proposed a new type of proportional electro-magnetic actuator for the fuel injection pump of a diesel engine. Before coming to the laboratory, he was an engineer with a Moscow-based company designing and producing embedded systems for inductrial applications. Mr. Secretta completed his M.A.Sc. degree at DalTech. His thesis proposed a new method for registration of multi-modal and time-sequential three-dimensional images. Using this method, two or more medical images acquired with dierent modalities are merged to construct a single image with higher information content. The method also allows
G. Secretta, P.H. Gregson / Computers in Biology and Medicine 29 (1999) 333±359
359
uni-modal three-dimensional registration of time-sequential images for detecting changes in pathology. These changes are used in both research and clinically to assess disease progression for prognostication and treatmen planning. Two other projects with which Gleb is heavily involved are the development of theory and algorithms for near-real-time counting and classi®cation of zooplankton and for counting and classifying droplets and bubbles in fast-¯owing, dynamic oil-water-air mixtures.
www.elsevier.com/locate/compbiomed
Automated registration of multimodal brain image sets using computer vision methods Gleb Secretta, Peter H. Gregson* Computer Vision and Image Processing Laboratory, Department of Electrical and Computer Engineering, DalTech, Dalhousie University, P.O. Box 1000, Halifax, NS, Canada B3J 2X4 Received 7 August 1998; accepted 13 May 1999
Abstract We present a new method of registering three dimensional image volumes of the brain of a given patient acquired at dierent times or with dierent imaging modalities, or both. Registration is an essential requirement for fusing the data from the two image sets so as to either increase the available information by exploiting complementary imaging modalities, or to measure small changes over time for prognostication, disease assessment, etc. The new technique exploits an external, removable, remountable reference frame which is attached to the head. Computer vision techniques are used to determine the positions of ®ducial marks in every image. The transformation required to map each image of one image volume onto the other image volume is developed using the theory of quaternions. The results indicate that the new technique is robust and practical in a clinical setting. # 1999 Elsevier Science Ltd. All rights reserved. Keywords: Image registration; Fiducial reference; MRI; SPECT; 3D imaging
1. Introduction Clinicians use comparative two- and three-dimensional imaging to quantify changes in tissue structure or function and to accurately localize functional information with respect to anatomical information. For tracking changes in structure such as the growth of some lesions, the clinician compares the apparent volume of the lesion as shown in MRI (magnetic * Corresponding author. Tel.: +1-902-494-6050; fax: +1-902-422-7535. E-mail address: [email protected] (P.H. Gregson) 0010-4825/99/$ - see front matter # 1999 Elsevier Science Ltd. All rights reserved. PII: S 0 0 1 0 - 4 8 2 5 ( 9 9 ) 0 0 0 1 3 - X
334
G. Secretta, P.H. Gregson / Computers in Biology and Medicine 29 (1999) 333±359
resonance imaging) image sets acquired at dierent times. In order for this comparison to yield clinically useful information early in the course of the disease, this comparison must result in accurate estimates of changes in lesion volume. Comparison of image sets is exceedingly dicult unless the two image sets are acquired with identical orientation and position. This is dicult to do in practice. SPECT (single photon emission computed tomography) imaging is used to detect areas of high metabolic rate such as tumors. The low spatial resolution of SPECT prevents accurate localization of the tumor site, however. An attractive solution is to use SPECT to acquire functional information and to use MRI to determine anatomical location. In this instance, the task is to present the functional information on the anatomical image, requiring accurate registration of dissimilar image sets. This task is especially demanding because dierent modalities usually have dierent resolutions as well as dierent coordinate system orientations. We present a new head-mounted reference frame and an algorithm for registration of two sets of brain images acquired with either dierent imaging modalities or at widely dierent times, or both. The new approach uses a new, light-weight, adjustable reference frame designed to ®t inside various imaging systems. Once adjusted for the particular patient, it can be repeatedly re-mounted on that patient with little positional error. Techniques from computer vision are used to detect the ®ducial marks in the images. A new algorithm which exploits the redundancy of the data is used for determining the vector-quaternion transformation which maps the coordinate system of one image set into the coordinate system of the other. Multilinear interpolation is applied to reduce quantization eects caused by dierent spatial resolutions of the data sets and by non-cubic voxels. The balance of this paper is organized as follows: Section 2 discusses the background to multi-modal image registration and describes SPECT and MRI imaging methods along with the restrictions imposed by each method. Section 3 presents the proposed approach to image registration. Experimental results and error analyses are given in Section 4. Conclusions and suggestions for future work are presented in Section 5.
2. Background In clinical practice and research, it is often necessary to determine both function and anatomy of brain structures of interest. Approaches for accomplishing this include functional MRI (fMRI) and `fusion' of SPECT and MRI image volumes. fMRI has the potential for providing both anatomical and functional information simultaneously. With fMRI, it is possible to determine the concentration of certain substances such as lactate that are the products of neurophysiological activity [1]. Detection of changes of lactate within speci®c brain regions in response to controlled stimulation should make it possible to use spectroscopic imaging to investigate stimulated brain function. There is much on-going research into fMRI but it is not currently a clinical tool. Further, the requirement for controlled stimulation makes fMRI an excellent research tool for establishing brain function but it is not suitable for localization of tumors and lesions. The most promising method of gaining both functional and anatomical information is based on `fusing' MRI and SPECT image sets [2]. However, the two modalities have dierent spatial
G. Secretta, P.H. Gregson / Computers in Biology and Medicine 29 (1999) 333±359
335
resolutions in all three directions and since the image sets are taken at dierent times with dierent imaging systems it is highly unlikely that the image planes of the two modalities would even be parallel. While some of the desired information can be determined from the two image sets without further processing, improved sensitivity to small lesions requires registering the two image sets by transforming the images of one set to match the three-dimensional spatial coordinates of the other set. In principle, three-dimensional image registration involves interpolating full three-dimensional information from the set to be transformed and then resampling it with the sampling parameters of the other set. We propose to incorporate these two steps into a single step so as to reduce computational burden and memory requirements. Applications for three-dimensional image registration include determining the anatomical locations of tumors found in SPECT imaging and comparing MRI image sets taken at signi®cant intervals for assessing disease progression. We present a short description of the two imaging modalities so as to establish the context for the proposed approach. 2.1. MR imaging MRI produces images based on the density of isolated and loosely-bound protons in the tissues of the body by measuring the intensity of a radio-frequency (RF) signal emitted by the tissue after an exciting magnetic ®eld and RF pulse are removed. This technique is based on the fundamental property that protons possess a magnetic moment and spin. When placed in a magnetic ®eld, the proton precesses about this ®eld. Initially all the protons are aligned either parallel or antiparallel to the magnetic ®eld. When a RF signal having an appropriate strength and frequency is applied to the object, the protons absorb energy, and more of them precess to the antiparallel state. When the applied RF signal is removed, the absorbed energy is reemitted and detected by an RF receiver. The proton density in the environment can be determined from the characteristics of the received signal. By suitably controlling the applied RF signal and the surrounding magnetic ®eld, and by appropriately processing the recovered RF signal, a three-dimensional image of the body can be formed as a set of slices. MRI is now widely used in clinical practice to image internal anatomy. Typical resolution of MRI machines is about 0.5 mm/pixel in the slice plane and about 5 to 7 mm between slices [3]. An example of an image acquired using MRI imaging is shown in Fig. 1. 2.2. SPECT imaging The goal of SPECT is to reconstruct, from projections, the three-dimensional distribution of radioactivity inside the patient's body. Dierent projections are acquired by counting the number of photons collected on a planar, collimated detector placed at various positions around the body. These detectors have point-spread functions (PSF) due to non-perfect collimation and to the attenuation of the photon's energy during its ¯ight from source to collimator. In addition, Compton scattering, gamma-ray attenuation and other physical phenomena make the data noisy thereby lowering SNR (signal to noise ratio). Various pre-reconstruction and post-reconstruction techniques have been used to enhance SPECT data such as constrained deconvolution, Wiener ®ltering and count-dependent Metz
336
G. Secretta, P.H. Gregson / Computers in Biology and Medicine 29 (1999) 333±359
®ltering [4]. Maximum entropy and homomorphic ®lters have also produced good results [5]. However, anatomical localization with SPECT is dicult because the method is designed to study physiological activity of the imaged organ and not its structure. In addition, SPECT images are often noisy and their spatial resolution ranges from 4 to 10 mm/pixel. SPECT data are characterized as being noisy, having poor resolution and showing function rather than anatomy [6]. An example of an image acquired using a SPECT scanner is shown in Fig. 1.
Fig. 1. Imaging with MRI and SPECT. Images a, b and c show MRI images. (a) Placement of the imaging planes. (b and c) Two adjacent slices. Images d, e and f show SPECT images. (d) Placement of the imaging planes. (e and f) Two adjacent slices. The dierence in spatial resolution is readily apparent.
G. Secretta, P.H. Gregson / Computers in Biology and Medicine 29 (1999) 333±359
Fig. 1 (continued)
337
338
G. Secretta, P.H. Gregson / Computers in Biology and Medicine 29 (1999) 333±359
2.3. The registration problem It is obvious that to solve the registration problem for sets of brain images, a transformation that maps one of the sets into the other must be found [7]. Since the brain doesn't appreciably deform or move inside the skull, we model it as rigid body. Thus, the transformation required is a simple translation±rotation combination and doesn't require anisotropic scaling. The desired transformation may be found using either anatomical landmarks in the two image sets or ®ducial marks based on an external reference [8]. Use of anatomical landmarks is not appropriate because the poor spatial resolution of SPECT imaging as seen in Fig. 1 can lead to a huge registration error. Even detection of the landmarks is highly problematic. Use of external ®ducial marks requires axing a reference frame external to the skull. The reference frame must be removable and re-axable in the same position so as to permit registration of images acquired at signi®cantly dierent times (perhaps months or years apart), as well as to avoid imaging-system scheduling problems (diculties with scheduling patients for both modalities on the same day). Ethical and practical considerations (availability and cost of operating suite time) prohibit the use of surgically-implanted markers. While this external frame is inconvenient for both the patient and the imaging process, it provides reliable information as to the orientation of the patient's head with respect to the imaging system.
3. The proposed approach The proposed registration procedure is: 1. 2. 3. 4.
Capture the two image sets with the external reference frame ®tted to the patient. Determine orientation of the head in each image set. ompute the transformation from one image set to the other. Use the transformation and interpolate voxels to determine values for the lower-resolution image set at the spatial locations of the higher-resolution image set.
The proposed reference system will appear on each slice in both sets as four bright spots outside the region of the patient's head (Fig. 9a). With the coordinates of the reference points in dierent slices it is possible to compute the orientation of the reference system with respect to the imaging system. Determining the transformation which maps one image set into another is straightforward with the aid of quaternion mathematics. Finally, bilinear interpolation is used to generate an image set derived from the lower-resolution set with the voxel coordinates (relative to the head) of the higher-resolution set. 3.1. Design of the external reference frame The external reference system consists of an adjustable frame and four detachable triangles. The head frame ®ts on the head and provides one keyed mounting location for a large triangle on each side of the head (Fig. 2). The frame is mounted on the head at the bridge of the nose and the auditory meati (ear canals). This provides a stable, reliable, repeatable mounting of the reference system on every
G. Secretta, P.H. Gregson / Computers in Biology and Medicine 29 (1999) 333±359
339
Fig. 2. (a) The head frame with triangles. Triangles are mounted in the sagittal plane. The head frame is adjustable and may be removed from the patient without disturbing adjustments. See text for explanation. (b) Detail of one triangle showing the mounting of the tube.
340
G. Secretta, P.H. Gregson / Computers in Biology and Medicine 29 (1999) 333±359
skull when properly adjusted. MRI imaging requires that no ferromagnetic parts be employed in the reference system and so the head frame and triangles are made of Plexiglas. Positioning of the head frame is provided by ®ve points of contact with the patient, two of which are adjustable. Non-adjustable points provide support of the reference system at the bridge of the nose (A) and two points on the upper part of the skull (B and C). Adjustable points D and E are `earplugs', each having three degrees of freedom to ®t into the auditory meati of the patient. After points D and E are properly adjusted for the patient, the reference system can be removed by undoing bolt F without disturbing the adjustment of points D and E. The frame may be remounted by ®tting it onto patient's head and retightening bolt F. Around the perimeter of each triangle is cut a narrow groove into which a ¯exible tube is mounted. This tube is ®lled with either a radio-emitting source to make its lumen easily detected in SPECT imaging or cod liver oil which renders its lumen easily detected in MRI. The triangles are otherwise identical and are made easily removable and replaceable so as to permit the appropriate pair of triangles to be used for each imaging modality. Accurate triangle positioning is assured because each triangle is keyed to the head frame. By making triangles removable, radiation exposure to patients and sta is minimized. The triangles are mounted in the sagittal plane as shown in Fig. 2. This is a matter of convenience for use with head coils and for the comfort of the patient. It also permits the use of two large triangles which maximizes the accuracy of head orientation measurements due to both triangle size and triangle redundancy. For both coronal and trans-axial imaging, sagittal placement of the triangles will result in ®ducial marks in most images. This placement of the triangles is not optimal for sagittal imaging, however, because very few images will contain ®ducial marks. A simple modi®cation to the current design of the head frame would permit sagittal imaging. The head frame is adjustable to ®t each patient. While our prototype does not have scales on each adjustment, scales could easily be provided on a production version of the device. This would make it possible to accurately re®t the head frame to patients on repeat visits for acquisition of sequential images.
3.2. Determining head orientation Consider the situation depicted in Fig. 3. In Fig. 3, one imaging plane is shown intersecting with both triangles. The points of intersection are shown as dots. Other imaging planes will intersect the triangles at dierent points, with the result that each side of each triangle will be described by a series of dots. By ®tting lines to the dots, the locations of the sides of each triangle can be estimated. Since the triangles are fastened to the patient's head by the head frame, head position and orientation can be determined. Clearly the method requires only one triangle; the second triangle provides redundancy to improve accuracy although in the current work it was used to check the results. The method depends on accurate determination of the intersection points of the reference triangles with the image slice planes; speci®cally, we must accurately detect the centers of the dots appearing in each image slice so as to accurately estimate the locations of the triangle sides.
G. Secretta, P.H. Gregson / Computers in Biology and Medicine 29 (1999) 333±359
341
Fig. 3. The intersections of the imaging plane with the perimeter tubes of the two reference triangles are indicated by the four large dots.
3.2.1. Finding dot centers Intersection of one image slice plane with the sides of the two triangles results in four dots appearing outside the head (Fig. 9a). To detect these dots, the operator selects either the left or the right side of the image set to provide the reference. This is done once per image set and is the only operator intervention required. Data from the other side are currently used for comparison, but may be incorporated into orientation determination in future. The dots are then isolated from the remainder of the image data. Since the reference system is designed so that the dots will appear about 30 mm from the head, a simple masking operation is used. The original image is low-pass-®ltered to achieve two goals. It must remove high-frequency noise and blur the image to `attach' those parts of the image that appear disconnected from the head (the ears or the nose can produce disjoint segments). The blurred image is thresholded at 10% of maximum intensity to produce a binary image. In this image, the head should be a single segment and the four dots are readily visible, disjoint from the head as required. Finally, area thresholding is employed to discard segments with areas either too small or too large, resulting in a binary `mask' image. This mask image is ANDed with the original image to produce a grey-level image in which only the dots appear. The lower and upper area thresholds (determined empirically) for MRI are 2 and 100 pixels and for SPECT are 5 and 50 pixels. Tubing that is not perpendicular to the plane of an image slice appears not as a disc but rather as a ®gure consisting of two elliptic and two straight segments due to the thickness of the slice (Fig. 4). Due to symmetry considerations, we assume that the center of a reference dot
342
G. Secretta, P.H. Gregson / Computers in Biology and Medicine 29 (1999) 333±359
represents the point that lies closest to the mid-line of the tubing, and thus can be used for determining the orientation of the corresponding triangle side. We have found that this procedure yields acceptable results, however we note that morphological processing and/or various `dot detectors' might yield better performance. 3.2.2. Validating dot images Since there can still be noise in the grey-level dot image, the algorithm must discriminate between valid and invalid candidate reference dots. We threshold the dot image at 30% of its maximum intensity. This results in suppression of some noise and makes the dots smaller. Resultant binary dot images are then analyzed further for the shape, number and relative position of dots. The program detects all the dots in the image and assigns each dot a unique label using a `connected components' algorithm. Coordinates of dot centroids are found for each dot as coordinates of its center of gravity. Binary dot images are used for this purpose, since the grey-level variation in the original image throughout the dot interior is very small. An acceptable dot image must meet the following criteria: 1. There must be exactly two dots in the reference half of the image. 2. The distance between the centroids of the putative reference dots must exceed 5 mm for MRI and 15 mm for SPECT images. 3. Areas of valid reference dots must be at least 3 pixels and not greater than 15 pixels for both image sets. 4. The standard deviations, measured in pixels in both the X and Y directions about each dot centroid, must not exceed 5 pixels. The ®rst three criteria were devised to detect various conditions that can appear in marginal slices, those slices which either do not contain dots because they do not intersect the reference triangles, or which contain images of a reference triangle near a vertex. Criterion 4 handles the
Fig. 4. View of the tubing crossing a slice. The tubing will appear as an ellipse due to the angle made by the tubing centerline with the slice and the thickness of the slice.
G. Secretta, P.H. Gregson / Computers in Biology and Medicine 29 (1999) 333±359
343
situation in which the imaging slice plane intersects the tubing wrapped around a triangle vertex, resulting in elongated dots which are unreliable for dot centroid calculations. 3.2.3. Estimating triangle sides The set of points due to the detection of the sides of one triangle in adjacent image slices constitutes a cloud of points. To estimate the location and orientation of this side, we ®t a straight line to this cloud. We assume a coordinate frame in which the X and Y axes lie in the image slice plane and Z is orthogonal to it. Physical constraints imposed by the possible orientations of the patient's head within the imaging system, the placement of the reference triangles on the head frame and the positioning of the head frame on the patient cause two sides of the reference triangle to meet the imaging plane at a large angle. This allows us to determine the orientation of the sides of each reference triangle by ®tting lines to the projections of its cloud of points on the Z±X and Z±Y planes. The parameters of these lines provide the information that we require. We use orthogonal regression to ®t lines to the projections of the cloud. With the general equation of the line in the 2D Z±X space in which z is the independent variable: az bx c 0
1
Line ®tting to a set of points by orthogonal regression in 2D assumes that each point is aected by independent identically distributed noise in both x and z directions. In that case the analysis tells us that: 1. The solution line must intersect the centroid of the points xi, zi, i = 1, . . ., n. The centroid coordinates xc, zc are x c , zc
N 1 X x i , zi N i1
2
P 2. The normal vector is the `smallest' eigenvector of the covariance matrix xz of the points, i.e. the eigenvector corresponding to the smallest eigenvalue. The covariance matrix is de®ned as X xz
N X x i ÿ x c 2 x i ÿ x c zi ÿ zc x i ÿ x c zi ÿ zc zi ÿ zc 2 i0
The covariance can be rewritten as X p q , q r xz P Then the smallest eigenvalue of xz is q p r ÿ p ÿ r2 4q2 l 2
3
4
5
344
G. Secretta, P.H. Gregson / Computers in Biology and Medicine 29 (1999) 333±359
The corresponding eigenvector is q lÿp Then the line ax+bz+c = 0 that has normal vector (q l ÿ p )T and intersects the point (xc, zc) is given by the parameters aq blÿp c ÿ ax c bzc :
6
3.2.4. Finding the reference vertex Only two sides of each triangle are used for determination of head orientation. Accuracy and precision of the orientation measurement requires that each triangle side is as close as possible to being parallel to the surface normal of the imaging plane, and further, that the sides meet at an angle that is not too acute. These con¯icting requirements and operational requirements on triangle size are satis®ed by only one vertex of each triangle, hereafter referred to as the reference vertex. In order to compute the translational component of the desired transformation, it is necessary to determine the coordinates of the reference vertex. This vertex should be at the point of intersection of the two lines which are estimates of the triangle sides. However, due to imaging errors, noise, and errors in locating dot centers, the two lines will rarely intersect in three-dimensional space. Thus, we take as the reference vertex the mid-point of the line orthogonal to both line estimates. 3.3. Computing the transformation The transformation desired consists of a translational component and a rotational component. The translational component is easily computed from the coordinates of one triangle vertex in each of the image sets. We convert pixels and voxels to millimeters in both image sets so as to use real units and also to avoid diculties with conversion of measurements. In our case we are dealing with a three-dimensional problem. Position of a rigid body in 3D space is determined by the position of two non-parallel lines belonging to the body, hence it is sucient to use only one triangle for the reference. Since the brain is assumed to be a rigid body for this research, we need not be concerned with anisotropic scaling and so we use unit vectors with the orientations of the lines. To determine the rotational part of the transformation, we need to ®nd two corresponding pairs of normalized, non-parallel vectors and the combination of rotations that transforms one pair into the other. The Euler theorem states that any two states of a rigid body with a single ®xed point are
G. Secretta, P.H. Gregson / Computers in Biology and Medicine 29 (1999) 333±359
345
separated by only one rotation about some axis. To ®nd this axis and the angle of the rotation we use the algebra of quaternions because it ®ts this class of problems extremely well [10]. Any quaternion q can be written in the following form: q q0 q1 i q2 j q3 k qs q~v ,
7
where i, j and k are the triple de®ning the Cartesian coordinate system, qs=q0 is the scalar part of the quaternion q and q~v =q1i+q2j+q3k is the vector part of the quaternion q. The norm of a quaternion q, written Nq is de®ned as q 8 Nq q20 q21 q22 q23 Thus, any quaternion q can also be written as: q Nq cos Q e~ sin Q,
9
where cos Q q0 =Nq 1 q q2 q22 q23 sin Q Nq 1 q1 i q2 j q3 k e~ 2 q q21 q11 q23 The normalized quaternion (q0=1) actually represents conical rotation of any vector a~ about an axis q~v passing through the origin of a~ for the angle of 6 q~v6 (Fig. 5). The fundamental theorem is given in Ref. [9]: Theorem 1. If q and r are any non-scalar quaternions, then
Fig. 5. Rotation in terms of quaternions.
G. Secretta, P.H. Gregson / Computers in Biology and Medicine 29 (1999) 333±359
346
r 0 qrqÿ1
10
is a quaternion whose norm and scalar are the same as for r. The vector V~r 0 is obtained by revolving V~r conically about V~q through twice the angle of q with respect to r (Fig. 5). Thus, if q Nq cos Q e~ sin Q,
11
then V~r 0 is obtained by revolving V~r conically about e~ through an angle 2Q. The normalized quaternion has Nq=1. Eq. (11) for a normalized quaternion is simpli®ed to q cos Q e~ sin Q:
12
Consider now the special case when vectors V~r and V~r 0 lie in the plane of rotation. Knowing vectors V~r and V~r 0 allows the angle between them and the orientation of vector e~ to be
Fig. 6. Finding the rotational part of the registration transformation as a combination of two rotations. (a) First rotation. (b) Second rotation.
G. Secretta, P.H. Gregson / Computers in Biology and Medicine 29 (1999) 333±359
347
computed. Vector e~ is simply their vector product. Thus, a quaternion can be de®ned using a pair of vectors. Now, a succession of rotations must be de®ned in terms of quaternions. Since the product of quaternions is also a quaternion then any number of rotations is equivalent to a single rotation. The succession of two rotations de®ned by quaternions p and q, respectively, corresponds to the operator r 0 qp rpÿ1 qÿ1 qp r qpÿ1
13
Thus, we can ®nd the quaternion which corresponds to the rotational part of the transformation that transforms one pair of vectors in one image set into the corresponding pair in the other image set. We de®ne our vectors in each set as vectors with origin at the reference vertex and with the orientation of each triangle side adjacent to that vertex. The estimates of the sides of the triangle are used to provide the pair of vectors in each image set. The problem is depicted in Fig. 6. There are two pairs of vectors V~1 , V~1 0 , V~2 and V~2 0 . Vector V~1 0 corresponds to V~1 , and vector V~2 0 corresponds to V~2 . The task at hand is to ®nd a normalized quaternion given two pairs of corresponding vectors. The solution is sought in the form of two successive rotations. Quaternions representing rotations can only be computed if corresponding vectors have the same length. To satisfy this requirement the longest vector is sought and the three other vectors are scaled up to have the same length. The ®rst rotation maps V~1 into V~1 0 . The axis of this rotation coincides with vector product V~1 V~1 0 . According to Eq. (12) the normalized quaternion that represents this rotation can be written as: q1 cos b
V~1 V~1 0 sin b jV~1 V~1 0 j
v !2 u V~1 V~1 V~1 V~1 u V~1 V~1 0 t 1ÿ q1 jV~1 jjV~1 0 j jV~1 V~1 0 j jV~1 jjV~1 0 j 0
0
14
15
The ®rst rotation results in V~2 going to V~2int . Now the second unit quaternion q2 that conically rotates V~2int into V~2 0 about axis V~1 0 can be calculated. As for the ®rst quaternion, two corresponding vectors that lie in the plane of rotation are used. The vectors P~2int and P~2 0 that are projections of V~2int and V~2 0 on the plane of rotation are used. Quaternion q2 that conically rotates V~2int into V~2 0 also rotates P~2int into P~2 0 and can be determined from them. Projection vectors P~2int and P~2 0 are found as V~2int V~1 0 ~ 0 P~2int V~2int ÿ V1 jV~1 0 j2
16
348
G. Secretta, P.H. Gregson / Computers in Biology and Medicine 29 (1999) 333±359
V~2 0 V~1 0 ~ 0 P~2 0 V~2 0 ÿ V1 jV~1 0 j2 Unit quaternion q2 is found according to Eq. (12) as v !2 u u 0 0 0 0 V~1 P~2 P~2int V~1 t P~2 P~2int sin g 1ÿ q2 cos g jV~1 0 j jP~2 0 jjP~2int j jV~1 0 j jPV~2 0 jjP~2int 0 j
17
18
Finally, according to Eq. (13), multiplication of quaternion q2 by quaternion q1 results in a quaternion q that represents rotation of vectors V~1 and V~2 into V~1 0 and V~2 0 respectively. Quaternion q is also a unit quaternion because it is the product of unit quaternions q1 and q2. Once this resultant quaternion q is found the rotation is de®ned. 3.4. Applying the registration transformation Since the transformation is carried out on a voxel-by-voxel basis, the low-resolution image is transformed to match the high-resolution image so as to retain as much image information as possible [11]. Thus, we transform the SPECT image set to match the MRI image set. The transformation is carried out by forming a vector from the reference vertex to the voxel under consideration in the MRI image set. If this vector is V~MRI , then by applying V~SPECT qV~MRI qÿ1 t~
19
where q is the quaternion representing the rotational part, t~ is the translational component of the transformation and vector V~SPECT in the SPECT image refers to the same anatomical location as V~MRI . The value found at this voxel in the SPECT image set is plotted at the MRI coordinates of the voxel under consideration in the transformed SPECT volume. Thus, the transformed SPECT image resulting from this procedure will have the SPECT data plotted with respect to the MRI coordinate axes, with the same anatomical points appearing at the same coordinates in each image set. The image sets are now in register. Two problems remain. 1. SPECT voxels are much larger than MRI voxels and so adjacent voxels in the transformed SPECT image have identical values, resulting in a strongly `pixelated' image. 2. Since the imaging planes of the SPECT data are not parallel to those of the MRI data, the transformed SPECT data exhibits visible edges due to data being taken from dierent planes. Both of these problems are evident in Fig. 8. To solve these problems, we use a multi-linear interpolation. Collections of new data points which all arise due to the same voxel in the original SPECT data are termed `supervoxels'. The value of the transformed voxel is computed as follows (Fig. 7): 1. Assume that the x and y coordinate axes lie in the slice plane and z is orthogonal to the slice plane. 2. Label the supervoxel containing the new voxel as K and the two nearest supervoxels on the
G. Secretta, P.H. Gregson / Computers in Biology and Medicine 29 (1999) 333±359
349
Fig. 7. Computing the grey scale value of a new voxel. (a) Labeling of the supervoxels. (b) Interpolation in the plane of the new voxel.
same plane as L and M. 3. Label supervoxels on the nearest plane to the new voxel and adjacent to K, L and M as K', L' and M ', respectively. 4. Using the x and y coordinates of the new transformed voxel, 4.1. Compute G as the linear interpolation of K and the average of L and M. 4.2. Compute G' as the linear interpolation of K'and the average of L' and M '. 5. Compute the new value of the transformed voxel as the linear interpolant of G and G'. This interpolation solves both pixelation and interplane transitions eects (Fig. 8).
350
G. Secretta, P.H. Gregson / Computers in Biology and Medicine 29 (1999) 333±359
Fig. 8. Results of interpolation. (a) Transformed SPECT image without interpolation illustrating `pixelation' and plane transitions. (b) Transformed SPECT image after multi-linear interpolation.
G. Secretta, P.H. Gregson / Computers in Biology and Medicine 29 (1999) 333±359
351
4. Experimental results 4.1. Head frame mounting accuracy Experiments were conducted to determine the registration accuracy of the new technique. Accurate registration requires that both the rotational and translational components of the transformation be correct and that the voxel values of the transformed image set are correctly computed. The ®rst step was to determine the accuracy with which the head frame could be repositioned. To perform this experiment, a volunteer's head was placed in a known position with the aid of a custom-®tted bite plate molded to ®t the volunteer's teeth. The bite plate was axed to a Plexiglas bar (5 cm wide, 5 mm thick, 15 cm long) which was clamped in a vise that was bolted to a table. The head frame was equipped with a forward-projecting laser. The volunteer positioned his head accurately by biting on the bite plate, moving his head until his teeth indexed positively on the bite plate. An assistant mounted the head frame on the volunteer, and the position of a small spot created by the laser on a tangent screen placed 1.53 m in front of the laser was recorded. The head frame was removed and remounted 30 times, with the location of the laser spot being recorded each time. After acquiring the data, the mean and standard deviation of the dot coordinates were computed. These values were computed only for the horizontal plane (motions due to trans-axial rotations) for two reasons. Firstly, the geometry of the head frame suggests, and the volunteer con®rms, that re-positioning accuracy should be lower with respect to transaxial misalignment than with respect to misalignments in the other two planes. Users report that the head frame mounting seems less ®rmly anchored to the head for coronal and transaxial rotations than for sagittal rotations. Secondly, the volunteer found it extremely dicult and very time consuming to bite down on the bite plate without causing its Plexiglas mounting to bend. For lateral (transaxial) rotations, the Plexiglas bar resisted bending very well as it was 5 cm wide. Further, lateral mis-alignments of the head were very obvious to the volunteer as they lead to discomfort in the jaws due to the teeth not ®tting the bite plate properly. In contrast, vertical mis-alignments resulted in signi®cant bending of the Plexiglas bar with relatively low forces applied because the bar was only 5 mm thick. Additionally, vertical motions did not result in loss of the ®t between teeth and the bite plate and so the volunteer had very little indication of vertical misalignment. Coronal rotations also lead to discomfort as under these conditions the Plexiglas bar is being twisted and the jaws feel as if they are being pried open. The results of the experiment are found in Table 1. The right-most column refers to the error distance measured on a cylinder of 200 mm radius, suciently large to enclose a head. Table 1 Results of the head frame remounting experiment
Standard deviation Maximum mag.
Distance at 1.53 m (mm)
Angular error (min, s)
Distance at 200 mm
2.41 5.97
5 '240 13 '230
0.31 0.78
352
G. Secretta, P.H. Gregson / Computers in Biology and Medicine 29 (1999) 333±359
We note that the maximum error magnitude is in the order of one voxel at this radius. Misalignment of the head frame will result in, at most, a one-voxel shift in the placement of the ®ducial reference marks in the images and, most of the time, a shift considerably less than this. Thus, the portion of the error due to head frame mis-alignment is generally greatly below one voxel. (Note that brain structures are more central and so the errors in alignment will be proportionally smaller.) We conducted MRI-to-MRI, SPECT-to-MRI and SPECT-to-SPECT studies to determine the accuracy of transformation. For assessing MRI-to-MRI registration, the reference triangles mounted on the head frame were ®lled with cod liver oil to generate bright spots in MRI images. To permit accurate re-positioning, we used a phantom consisting of a U-shaped piece of ABS pipe of 2-inch diameter, ®lled with water and rigidly attached to the head frame. This phantom shape allowed us to position the head frame precisely with the aid of a laser, while simultaneously providing and image-able object with sucient spatial complexity that we could verify our test data. For SPECT-to-SPECT and SPECT-to-MRI comparisons, we used a more physiologicallyrelevant phantom consisting of a plastic canister reinforced with plastic rods and containing ®llable compartments to emulate structures in the brain. The main part of the canister was ®lled with radioactive technetium (Tc). A phantom was used because ethical considerations precluded the use of a volunteer in these studies. For the SPECT-related studies, the tubes of the two reference triangles were ®lled with a radioactive solution with a higher concentration than the solution in the canister. Three SPECT data sets at diering orientations were acquired. Finally, two MRI image sets of a volunteer were acquired to illustrate algorithm performance with real images. 4.2. Estimating accuracy and precision To assess accuracy, we need two sets of images of one object in dierent poses. We call these set 1 and set 2. We wish to evaluate the dierences between set 1 and set 2 after set 2 has been transformed to match the pose of set 1. Ideally, there would be no dierence. The ®rst study assessed the accuracy of MRI-to-MRI registration. Our MRI scanner produced in-plane resolution of 0.87 mm/pixel and interslice spacing of 5.2 mm. For this study, the phantom was mounted on a test jig that permitted only out-of-image-plane rotation (the worst case) of the phantom. V-shaped notches were made in the test jig to denote angles to which the phantom with the attached headframe was to be rotated. After imaging in the MRI scanner, the actual angles were determined by ®tting a laser to the head frame to provide an optical lever, reproducing the angular positions of the head frame, noting the positions of the resulting laser spots on a tangent screen and computing the resulting angles. Repeated repositioning trials produced a positioning error of 21 '300. Four MRI scans of the phantom in four dierent angular positions were performed. Rotations were computed using each of the two triangles for each of the six pair-wise combinations of the four images (Table 2). From the data, it is clear that the rotations estimated by the two triangles are in close agreement with actual rotation angles and with each other. In Table 2, MR1 refers to MRI image 1. The columns headed `left triangle' and `right
G. Secretta, P.H. Gregson / Computers in Biology and Medicine 29 (1999) 333±359
353
Table 2 Angles of rotation in degrees for MR-to-MR registration Registration
Left triangle
Right triangle
Actual angle
Maximum error
MR1±MR2 MR1±MR3 MR1±MR4 MR2±MR3 MR2±MR4 MR3±MR4
11.81 8.59 11.25 7.42 19.46 12.47
11.78 8.38 11.39 7.37 19.10 11.98
11.75 8.22 11.30 6.98 19.72 12.93
ÿ 0.06 ÿ 0.37 ÿ 0.09 ÿ 0.44 0.62 0.95
triangle' present the angular dierence between the poses of the object (the phantom) in the two images as computed using the left triangle and the right triangle, respectively. The actual rotation angle and the maximum error are also given. The small rotation angles used represent physiologically-relevant positioning dierences between imaging sessions and so are realistic. The worst-case error on the surface of a cylinder with a radius of 200 mm is thus 3.3 mm or 3.8 pixels. The absolute average of the worst-case errors was 0.428 which is 1.47 mm or 1.69 pixels on the cylinder. The angular standard deviation of all estimates was 0.448 which is 1.52 mm or 1.7 pixels on the 200 mm cylinder. While this is not strictly the precision of the system (it does not look at the `spread' of the errors for all pixels within one pair of image sets but rather the spread for six imageset pairs), it does provide some indication of repeatability. These results are clinically acceptable as changes that might be overlooked due to errors of this size are not considered clinically signi®cant. The second study looked at SPECT-to-MRI registration using the canister phantom. One MRI image set and three SPECT image sets were acquired. The SPECT system captures 64 64 pixel images, with pixels about 4 mm square. The consistency with which the SPECT image sets were brought into register with the MRI set are shown in Table 3. We were not able to determine accuracy because to do so we need to accurately determine the inherent coordinate systems of both scanners. The coordinate system of the MRI scanner may be modi®ed by changes to pulse sequences and it is also aected by minor disturbances in its magnetic ®eld and so the reliability with which we could determine its coordinate frame was considered dubious.
Table 3 Angles of rotation for the SPECT-to-MR registration Registration
Left triangle
Right triangle
Angle error
SP1±MR SP2±MR SP3±MR
5.2 5.3 11.3
5.4 5.6 11.9
0.2 0.3 0.6
354
G. Secretta, P.H. Gregson / Computers in Biology and Medicine 29 (1999) 333±359
Table 4 Angles of rotation for the SPECT-to-SPECT registration Registration
Left triangle
Right triangle
Angle error
SP1±SP2 SP1±SP3 SP2±SP3
7.8 11.6 7.2
7.5 11.1 7.9
0.3 0.5 0.7
Comparison of the transformation angles estimated with right and left triangles shows that they agree well. This shows that the estimates are consistent. For this study, the worst case discrepancy is 0.68 which is 2.1 mm or about 0.5 pixels on the surface of the 200 mm cylinder. The consistency data for the SPECT-to-SPECT registration study appear in Table 4. The worst-case dierence in angular position is 0.78. This corresponds to 2.4 mm on the 200 mm cylinder or 0.6 pixels. This is clinically acceptable precision. Finally, real data were analyzed to illustrate the eectiveness of registration. Two MRI image sets of a volunteer were acquired and registered. Real data are useful because they have higher spatial frequencies than our phantom and so any mis-registration is clearly visible. Fig. 9(a) shows one image from the ®rst image set of the pair. Fig. 9(b) shows the corresponding plane, after the registration transformation, taken from the second image set of the pair. The algorithm's ability to retain relatively high spatial frequencies is clear. The sulci between the gyri are retained to a large degree and the distinction between grey matter and white matter is also clear. There is some loss of resolution; this is unavoidable if we are to prevent `pixelation' and inter-plane discontinuities in the image slices. The small black wedge at the bottom of the image represents the amount in the x±y plane by which the registration algorithm was required to transform the second image set. Not shown is the out-of-plane rotation of the head for which compensation was also accomplished. Fig. 9(c) shows the dierence or error image. The grey levels in this image are not meaningful as our display stretches contrast. The degree of mis-registration is visible in Fig. 9(c), appearing to be a sort of `double exposure' of the head. Note that should registration be perfect, the image would be a uniform grey level. While not perfect, the registration error is very small, certainly compared to the error without registration. Typical maximum deviation from the mean grey level inside the region of interest (the brain) in the error image is 230 grey levels, or 212% of the whole range of 256 grey levels, while the standard deviation for the same region equals 5 grey levels or 2%.
5. Conclusions and future work A new reference frame and a new algorithm have been developed to enable rapid, clinicallyuseful registration of image volumes both between imaging modalities and within a given modality for sequential scanning. Applications of the these developments include integrating
G. Secretta, P.H. Gregson / Computers in Biology and Medicine 29 (1999) 333±359
Fig. 9. MR-to-MR registration. (a) Benchmark image. (b) Registered image. (c) Result of subtraction.
355
356
G. Secretta, P.H. Gregson / Computers in Biology and Medicine 29 (1999) 333±359
Fig. 9 (continued)
functional and anatomical data, monitoring disease progression and planning invasive therapies. The new approach has proved to be an accurate tool for registration of unimodal and multimodal three-dimensional data sets. The head frame provides sucient adjustment range while preserving an acceptable degree of positional repeatability when remounted. Performance of the system, consisting of the head frame and the algorithm, is sucient for many applications where 3D registration is required. MRI-to-MRI registration error was under 2 mm for the region of interest (the head, as exempli®ed by the 200 mm cylinder). These errors are not considered clinically meaningful. The maximum discrepancy between estimates of angular transformation for the SPECT-toSPECT registration was 2.4 mm for our SPECT scanner. SPECT-to-MRI registration exhibited similar performance. Among the approaches to the registration problem that are reported in current literature there are two that reach similar accuracy. One of them proposes use of a combination of ®ducial markers and structured light to achieve registration [11] and has RMS error 1.1720.45 mm. This method implies use of certain auxilliary equipment such as special light projector and videocamera to capture the pattern of light on the patient's body. The other method does not require any additional equipment [12]. The algorithm minimizes a mutual information metric in an iterative way and provides suboxel accuracy. One drawback of this approach is the considerable amount of processing time needed for the algorithm to converge, which typically
G. Secretta, P.H. Gregson / Computers in Biology and Medicine 29 (1999) 333±359
357
is between 30 and 60 min. The other problem with this method is that it requires that images of both modalities have high resolution. Our algorithm uses minimal auxilliary equipment and is independent of the image resolutions in the data sets. However, further work is required to more accurately determine the angle of rotation, as we suspect that much of our error derives from noise and mislocalization of reference dots in cases when most of them fall close to the vertex of the triangle. One way to overcome this drawback is to increase the dimensions of the triangles and mount them closer to the patient's head. The system should also be modi®ed to permit image acquisition in all three planes. Finally, the transformation should be estimated using both triangles instead of only one. This would provide more measurements. It would also permit making measurements over greater displacements (between triangles) thereby increasing accuracy.
6. Summary We present a new method of registering three dimensional image volumes of the brain of a given patient acquired at dierent times or with dierent imaging modalities (MRI and SPECT), or both. Registration is an essential requirement for fusing the data from the two image sets so as to either increase the available information by exploiting complementary imaging modalities, or to measure small changes over time for prognostication, disease assessment, etc. The new technique exploits an external, removable, re-mountable reference frame which is attached to the head. The head frame is suciently adjustable to ®t most patients, the positional accuracy with which it may be remounted indicates that it contributes negligibly to the registration errors. The head frame is equipped with two removable triangles which have perimeter tubes containing contrast material (MRI) or radioactive material (SPECT). These triangles generate four dots in each image which intersects them. The position and spacing of these dots are dependent on the orientation of the head frame, and therefore the head, with respect to the imaging plane. These dots form the ®ducial references. Using computer vision techniques, the dots are automatically located in each image. The locations of the triangle sides are estimated from the set of dots displayed in all the images of the image set. The location of a reference vertex of the triangle is determined. This procedure is followed for both image sets. Quaternion mathematics is used to determine the rotational transformation which carries one head pose into the other. This transformation is then applied to each pixel in the image set. Bilinear interpolation is used to remove `pixelation' eects due to non-parallel imaging planes in the two image sets and to dierent spatial resolutions of imaging modalities. The new algorithm is tested with MRI-to-MRI, SPECT-to-MRI, and SPECT-to-SPECT registration tasks. Performance is deemed to be clinically acceptable, with registration errors below 2 mm throughout the head volume. SPECT-based registration is consistent to within 6 mm.
358
G. Secretta, P.H. Gregson / Computers in Biology and Medicine 29 (1999) 333±359
Finally, MRI-to-MRI registration was performed on the images of a volunteer. These results show that the method preserves much of the anatomical information present in these images. References [1] M. Singh, Towards proton MR spectroscopic imaging of stimulated brain function, IEEE Trans. Nucl. Sci. 39 (4) (1992) 1161. [2] G.K. Matsopoulos, S. Marshall, J.N.H. Brunt, Multiresolution morphological fusion of MR and CT images of the human brain, IEEE Proc.-Vis. Image Signal Process. 141 (3) (1994) 137±143. [3] R. Sauter, M. Shneider, K. Wicklow, H. Kolem, Current status of clinically relevant techniques in magnetic resonance spectroscopy, Siemens Electro-Medica 60 (29-56) (1992) 32. [4] D. Boulfelfel, R.M. Rangayan, L.J. Hann, R. Kloiber, Prereconstruction restoration of myocardial single photon emission computed tomography images, IEEE Trans. Med. Imaging 11 (3) (1992) 336±399. [5] W.A. Kalender, Quo Vadis CT? CT in the year 2000, Siemens Electro-Medica 61 (2) (1993) 30. [6] O. Kuler, G. Gerig, Computer-analysis of 3D medical images, Siemens Rev. (1991) 4. [7] P.J. Besl, N.D. McKay, A method for registration of 3D shapes, IEEE Trans. Pattern Anal. Mach. Intell. 14 (2) (1994) 239±256. [8] W.R. Fright, A.D. Linney, Registration of 3D head surfaces using multiple landmarks, IEEE Trans. Med. Imaging 12 (3) (1993) 515. [9] L. Brand, Vector and Tensor Analysis, in: John Wiley and Sons, NY, 1948, pp. 403±421. [10] R.J. Moorhead II, Z. Zhu, Signal processing aspects of scienti®c visualization, IEEE Signal Process. Mag. 12 (5) (1995) 20. [11] R.B. Backman, A temporal 3D-registration framework for computer-integrated surgery, Ph.D. thesis, Department of Computer Science, Robotics and Vision Group, University of Western Australia, 1998. [12] C.R. Meyer, J.L. BoeCR, B. Kim1, P.H. Bland, K.R. Zasadny, P.V. Kison, K. Koral, K.A. Frey, R.L. Wahl, Demonstration of accuracy and clinical versatility of mutual information for automatic multimodality image fusion using ane and thin-plate spline warped geometric deformations, Medical Image Analysis 1(3) (1996± 1997) 195±206. Dr. Peter Gregson (B.Eng. 74, M.Eng. 77, Ph.D. 88) is a Professor in the Department of Electrical Engineering at DalTech and is Director of the Computer Vision and Image Processing Laboratory. He has founded three companies designing and manufacturing data acquisition and control systems and was employed by the Defense Research Establishment. Dr. Gregson's research is focused on the development of new theory, architectures and algorithms for `early' computer vision and image processing and on applications of this technology to medical imaging, robot navigation and industrial control and measurement. He is developing theory, algorithms and analog VLSI circuitry for real-time, vision based short-range robot navigation, for automated surveillance monitoring and for motion estimation. Current medical imaging projects include automated assessment of diabetic retinopathy; automated determination of refractive error in pre-verbal children; automated correction of motion artifact in CT images; detection and quanti®cation of prostrate tumors and benign prostate hyperplasia in 3D ultrasound; multi-modal image registration and detection, localization and quanti®cation of multiple sclerosis lesions in 3D MRI. Dr. Gregson's laboratory has several ongoing industrial collaborations. Dr. Gregson has a major interest in eective, relevant teaching. He was the recipient of the Wighton Fellowship in 1996. This annual award is made to one professor of Engineering and Applied Science in Canada for innovative undergraduate laboratory presentation.
Mr. Gleb Secretta is a Ph.D. student in the department of Electrical Engineering at DalTech. He is a member of the Computer Vision and Image Processing Laborator. Mr. Secretta was born in 1967 in Moscow, Russia. He graduated with Bachelor of Engineering and Master of Engineering degrees in Systems and Control from the Moscow State Technical University in 1991. His Russian degree thesis proposed a new type of proportional electro-magnetic actuator for the fuel injection pump of a diesel engine. Before coming to the laboratory, he was an engineer with a Moscow-based company designing and producing embedded systems for inductrial applications. Mr. Secretta completed his M.A.Sc. degree at DalTech. His thesis proposed a new method for registration of multi-modal and time-sequential three-dimensional images. Using this method, two or more medical images acquired with dierent modalities are merged to construct a single image with higher information content. The method also allows
G. Secretta, P.H. Gregson / Computers in Biology and Medicine 29 (1999) 333±359
359
uni-modal three-dimensional registration of time-sequential images for detecting changes in pathology. These changes are used in both research and clinically to assess disease progression for prognostication and treatmen planning. Two other projects with which Gleb is heavily involved are the development of theory and algorithms for near-real-time counting and classi®cation of zooplankton and for counting and classifying droplets and bubbles in fast-¯owing, dynamic oil-water-air mixtures.