Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

3D spherical harmonic invariant features for sensitive and robust quantitative shape and function analysis.. 2007

You don't seem to have a PDF reader installed, try download the pdf

Item Metadata

Download

Media
ubc_2008_spring_uthama_ashish.pdf
ubc_2008_spring_uthama_ashish.pdf
ubc_2008_spring_uthama_ashish.pdf [ 7.08MB ]
ubc_2008_spring_uthama_ashish.pdf
Metadata
JSON: 1.0066294.json
JSON-LD: 1.0066294+ld.json
RDF/XML (Pretty): 1.0066294.xml
RDF/JSON: 1.0066294+rdf.json
Turtle: 1.0066294+rdf-turtle.txt
N-Triples: 1.0066294+rdf-ntriples.txt
Citation
1.0066294.ris

Full Text

3D Spherical Harmonic Invariant Features for Sensitive and Robust Quantitative Shape and Function Analysis in Brain MRI by Ashish Uthama B. E. (Electronics and Communication Engineering), Visvesvaraya Technological University, 2003 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF Master of Applied Science in THE FACULTY OF GRADUATE STUDIES (Electrical & Computer Engineering) THE UNIVERSITY OF BRITISH COLUMBIA October 2007 Ashish Uthama, 2007 Abstract A novel framework for quantitative analysis of shape and function in magnetic resonance imaging (MRI) of the brain is proposed. First, an efficient method to compute invariant spherical harmonics (SPHARM) based feature representation for real valued 3D functions was developed. This method addressed previous limitations of obtaining unique feature representations using a radial transform. The scale, rotation and translation invariance of these features enables direct comparisons across subjects. This eliminates need for spatial normalization or manually placed landmarks required in most conventional methods [1-6], thereby simplifying the analysis procedure while avoiding potential errors due to misregistration. The proposed approach was tested on synthetic data to evaluate its improved sensitivity. Application on real clinical data showed that this method was able to detect clinically relevant shape changes in the thalami and brain ventricles of Parkinson's disease patients. This framework was then extended to generate functional features that characterize 3D spatial activation patterns within ROIs in functional magnetic resonance imaging (fMRI). To tackle the issue of intersubject structural variability while performing group studies in functional data, current state-of- the-art methods use spatial normalization techniques to warp the brain to a common atlas, a practice criticized for its accuracy and reliability, especially when pathological or aged brains are involved [7-11]. To circumvent these issues, a novel principal component subspace was developed to reduce the influence of anatomical variations on the functional features. Synthetic data tests demonstrate the improved sensitivity of this approach over the conventional normalization approach in the presence of intersubject variability. Furthermore, application to real fMRI data collected from Parkinson's disease ii patients revealed significant differences in patterns of activation in regions undetected by conventional means. This heightened sensitivity of the proposed features would be very beneficial in performing group analysis in functional data, since potential false negatives can significantly alter the medical inference. The proposed framework for reducing effects of intersubject anatomical variations is not limited to functional analysis and can be extended to any quantitative observation in ROIs such as diffusion anisotropy in diffusion tensor imaging (DTI), thus providing researchers with a robust alternative to the controversial normalization approach. 111 Table of Contents Abstract^ ...ii Table of Contents  ^ iv List of Tables  ^ viii List of Figures^ ix Statement of Co-authorship^ xii 1 Thesis Introduction and Summary^ 1 ^1.1^Thesis Motivations and Objectives  2 ^ 1.2^Principles of Magnetic Resonance Imaging (MRI)^ 3 1.3^Structural Brain MRI^ 8 1.3.1 Summary of Structural Data Acquired^  10 1.3.2 MR Image Preprocessing^  10 1.3.3 Segmentation of Regions of Interest in the Brain^  12 1.3.4 Quantifying Structural Changes in Brain ROIs  15 1.4 Functional Magnetic Resonance Imaging^ 17 1.4.1 Principles of BOLD imaging^ 17 1.4.2 Hemodynamic Response 18 1.4.3 Functional Experiment Setup^ 19 1.4.4 Summary of Functional Data Acquired^ 22 1.4.5 fMRI Data Preprocessing ^ 23 1.4.6 Functional Activation Map Generation for Individual Subjects^ 27 1.4.7 Differences in Activation Maps at the Group Level^ 29 1.5^Thesis Contributions^ 32 1.5.1^Structural ROI Analysis  33 1.5.2 Functional ROI Analysis ^  34 1.6^Thesis Overview^ 35 iv 1.7^References^ 37 3 2 Unique Invariant SPHARM Features for ROI Based Analysis ^ 47 ^2.1^Introduction^ 47 ^ 2.2^Methods 48 2.2.1 Proposed Invariant SPHARM features^ 49 2.2.2^Statistical Analysis of Features 52 2.3^Data and Imaging^ 53 2.3.1^Synthetic Structural Data Generation^  53 2.3.2 Real Structural MR Data Acquisition  54 ^2.4^Results and Discussion^  55 2.4.1^Validation on Synthetic Data 55 2.4.2 Analysis of Real MR Structural Data^ 57 ^2.5^Conclusions^ 59 2.6^References 61 Detection and Analysis of Anatomical Shape Changes in MRI (PD thalamus) ^ 63 3.1 Background^ 63 3.2 Results 67 3.3 Discussion^ 70 3.4 Conclusions 72 3.5 Methods^ 73 3.5.1 MR Imaging at the University of British Columbia^ 73 3.5.2 MR Imaging at the University of North Carolina 74 3.5.3 Thalami Shape Analysis^ 74 3.5.4 Shape Analysis -- Technical Aspects^ 75 ^3.6^Authors' contributions^  79 3,7^Acknowledgements 79 v 3.8^References^ 80 4 Detection and Analysis of Anatomical Shape Changes in MRI (PD brain ventricles)^ 87 4.1^Introduction^  87 4.2^Methods 89 4.2.1 Proposed Invariant SPHARM features^ 89 4.2.2^Statistical Analysis of Features 93 4.2.3 Automated segmentation of the lateral ventricles^  94 4.3^Data and Imaging^ 95 4.3.1^Synthetic Structural Data Generation^ 95 4.3.2 Real Structural MR Data Acquisition 97 4.4^Results and Discussion^ 97 4.4.1 Validation on Synthetic Data 97 4.4.2 Analysis of Real MR Structural Data^ 99 4.5^Conclusions^  101 4.6^References  102 5 Extending SPHARM Features to Functional Activation Analysis in fMRI ^ 105 5.1^Introduction^  105 5.2^Methods  108 5.2.1 Invariant SPHARM-Based Spatial Features^  109 5.2.2 fMRI Analysis Using SPHARM Features  113 5.2.3^Statistical Group Analysis ^  115 5.3^Data and Imaging^  116 5.3.1^Synthetic fMRI Data  116 5.3.2 fMRI Data of PD Patients^  120 5.4^Results and Discussion  120 v i 5.4.1 Validation on Synthetic Data^ 121 5.4.2 Application to Clinical fMRI Data of PD Patients ^ 123 5.5 Conclusions^ 124 5.6 References 126 6 Validation and Comparisons to ConventionalflWRI Analysis Approach 129 6.1^Introduction^ 129 6.2 Methods 133 6.2.1 Invariant SPHARM-Based Spatial Features^ 133 6.2.2 fMRI Group Analysis Using SPHARM Features 138 6.2.3 Statistical Group Analysis ^ 140 6.2.4 ROI Normalization Approach (ROI-N)^ 141 6.3 Data and Imaging^ 142 6.3.1 Synthetic fMRI Data 142 6.3.2 Real fMRI Data^ 144 6.4 Results and Discussion 147 6.4.1 Validation on Synthetic Data^ 147 6.4.2 Application to Clinical fMRI Data 152 6.5 Conclusions^ 155 6.6 References 156 7 Conclusions and Future Work^ 161 7.1 Technical Contributions 161 7.2 Impact on Neuroscience Research^ 164 7.3 Future Work^ 165 7.4 References 166 8 APPENDIX^ 169 vii List of Tables Table 1.1 MR relaxation times for commonly occurring brain tissue^ 9 Table 1.2 MR relaxation time for Oxygenated and Deoxygenated blood^ 18 Table 2.1 Results of pre-drug vs. post-drug analysis^ 58 Table 2.2 Results of PD patients vs. controls analysis 59 Table 3.1 Results of volumetric and shape analysis. Numbers indicate the p-values obtained from the permutation test^ 68 Table 4.1 Result of pre-drug vs. post-drug analysis ^  100 Table 4.2 Result of comparison between control subjects and PD patients. Values indicated with a * are statistically significant at an alpha value of 0.05^ 100 Table 5.1 Permutation test results (p value) of fMRI data analysis from the Parkinson's disease study (the two groups being 'before drug' and 'after drug'). IG: Internally guided, EG: Externally guided, L: Left, R: Right, CRB: Cerebellar hemisphere, CAU: Caudate. Numbers in bold are statistically significant (alpha = 0.05).   124 Table 6.1 Maximum radius of ROIs and the corresponding bandwidth used ^ 152 Table 6.2 SPHARM analysis of the two frequency tasks. Only ROIS with at least one significant result are shown^  153 viii List of Figures Figure 1.1 Gradient fields along the X and Y axes^ 7 Figure 1.2 Whole brain images obtained using different relaxation time weighted protocols showing the different tissue contrasts. Left: T i weighted image. Right: T2 weighted image. 9 Figure 1.3 Left: A raw T i weighted image (axial slice), notice the noise near the nasal cavity and the eyes (upper part of image). Right: The extracted brain using the combined expectation maximization and geodesic active contour method proposed by Huang et al [33].   11 Figure 1.4 Result of anisotropic diffusion filtering on a T1 brain volume (one slice shown).^  12 Figure 1.5 Result of the automatic segmentation of the right brain ventricle. The left venctricle's binary ROI mask is rendered on orthogonal slices of the subjects Ti image.14 Figure 1.6 The MR-compatible rubber squeeze bulb used in the functional experiment. 20 Figure 1.7 Visual feedback given to subject performing the bulb squeezing task. The width of the black horizontal bar reflected the amount of 'squeezing' that the bulb recorded. Subjects had to modulate the squeezing action to keep the width within the undulating white pathway. 21 Figure 1.8 The stimulus design used for the functional study. The tunnel width shown on the Y axis indicates the width of the undulating pathway shown in Figure 1.7. ^ 21 Figure 1.9 Left: A slice from a T 1 weighted anatomical scan of a subject. Right: Corresponding slices of 12 * scans from functional volumes acquried the duration of the experiment (130 such volumes were collected over a duration of 260 seconds) 23 Figure 1.10 The process of brain extraction, followed be co-registration of the high resolution T1 anatomical scan to the first functional T2 * voume^ 24 Figure 1.11 Motion correction re-aligns the functional volumes in case the subject moves during the time course of acquisition^ 26 Figure 1.12 The process of obtaining subject level statistical parameter maps. The time course of each voxel is statistically compared to an expected response based on the task perfored to obtain a parameter which describes the activation at that voxel. This process is reapeated for all voxels in the brain resulting in the subject level SPM. 28 Figure 1.13 The functional ROIs studied in this thesis. All ROIs shown occur in pairs (in the left and righ brain hemishpheres), only regions in left hemisphere are shown. Indicated ROIs are Prefrontal cortex (PMC), Supplementary motor cortex (SMA), ix Primary motor cortex (Ml), Caudate (CAU), Putamen (PUT), Thalamus (THA), Cerebellar hemisphere (CER).^ 32 Figure 2.1 (a) The human brain ventricles. The two ventricles would have to be analyzed jointly since the relative orientations might be of importance. One of the spherical shells at r =35 is also shown. (b) The parameterization of the surface at r = 16 (top) and 35 (b ttom) 52 Figure 2.2 (a) Surface rendering of a real left thalamus (top) and right thalamus (bottom) from a PD patient, note the rough edges. (b) Group A synthetic samples (see text) generated with parameter set (2, 5, 10). (c) Group B synthetic sample generated with parameter set (2, 5, 8). Note the realistic looking edges in the synthetic data 54 Figure 2.3 Plot of the Euclidian distance between the two groups A and B for different values of cb. Distances resulting in a significant difference (alpha = 0.05) between the groups as obtained by the permutation test are indicated with a star. Non-Significant distances are indicated with a dot.  56 Figure 3.1 The SPHARM-based method for shape determination. The shape to be specified (the thalamus) and two concentric spherical shells are shown. On the right is the intersection between the thalamus and shells as a function of rotation (0) and azimuth (T). The rotation angle spans from 0 to 2it radians, and the azimuth angle is from 0 to IC r dians. 67 Figure 3.2 Two sets typical thalami registered and shown on brain from a PD subject. Note that the registration here is solely for visualization purposes, and is not required for the calculation of shape differences. Also, although the thalami here were first smoothed with a 12mm FWHM Gaussian kernel for visualization purposes, no smoothing was performed for the shape analysis and group comparison 69 Figure 4.1 Intersection of a left brain ventricle with 2 of the 128 shells used to obtain the SPHARM representation. The corresponding sampled surfaces of the two shells are also sh w . 93 Figure 4.2 (a) Surface rendering of a real left thalamus (top) and right thalamus (bottom) from a PD patient, note the rough edges. (b) Group A synthetic samples (see text) generated with parameter set (2, 5, 10). (c) Group B synthetic sample generated with parameter set (2, 5, 8). Note the realistic looking edges in the synthetic data 96 Figure 4.3 . Plot of the Euclidian distance between the two groups A and B for different values of cb, (ca was fixed at 10 for all points). Distances resulting in a significant difference (alpha 0.05) between the groups as obtained by the permutation test are indicated with a star. Non-Significant distances are indicated with a dot.  98 Figure 5.1 A cross section of fMRI activation statistics from the right cerebellar hemisphere of a real subject^  118 Figure 5.2 A cross section of synthetic fMRI activation statistics for group B from different sets. A pattern consisting of two activation levels corrupted by noise is present, - delta values closer to the functional centre and +delta away from the centre. (a) and (c) represent extreme cases where the pattern present is visually discernible in the presence of noise. (b) and (d) represent very low signal strengths but are still detected by the proposed features. The delta values correspond directly to those in Figure 5.4. A Gaussian noise source of zero mean and o of one was used in the above figures  119 Figure 5.3 A cross section of fMRI activation statistics for group A synthetic data. Note the absence of a pattern  119 Figure 5.4 Results of the synthetic data experiment. The normalized Euclidian feature distance between groups A and B for each set is plotted against delta for the 3DMI, SPHARM-fs and SPHARM-f features. The black diamond markers represent the distances which yielded significant difference using 3DMI features, the red triangles represent SPHARM-fs features while the blue pentagrams represent SPHARM-f features. The closer the markers are to the zero delta mark, the smaller the failure range, hence better the performance. 3DMI has a failure range of (-1.5, 1.25), SPHARM-fs has a range (-1.1, 1.3), while SPHARM-f features return the best results by being able to discriminate the groups even at every low values of delta (-0.2, 0.7)  121 Figure 6.1 Representing data as a set of spherical functions. Real fMRI data from the right cerebellar hemisphere of a normal subject performing the maximum frequency task (Section III-B) is shown (t>2 for clarity). The shells were obtained using a bandwidth of L=128. The physical values of both angles are non-uniformly spaced (Chebyshev nodes) [19]. The threshold, bandwidth and the depicted transparent ROI boundary are for illustration purposes only.  135 Figure 6.2 (a) Cross section of a synthetic pattern injected into the binary mask of a real right cerebellar hemisphere. SNR = 1.204 dB. (b) Cross section of a real fMRI activation statistic map mapped to the same color map  144 Figure 6.3 False positive rates using SPHARM-f features for various values of the parameter L (Bandwidth). a = 0.05.  148 Figure 6.4 The average SNR plotted against increasing delta value.  149 Figure 6.5 Sensitivity of the SPHARM approach. The mean Euclidean distance between the feature vectors of a group with only noise and of another group with increasing SNR is plotted against the SNR. The star and triangle markers denote distance which yielded a p value less 0.05, signifying that the method detected differences between the groups. 150 Figure 6.6 Sensitivity of the ROI-normalization approach. The number of voxels surviving various height thresholds is plotted against the SNR. A cluster size threshold of 4 connected voxels was used in each case. These results were computed from the exact same data generated for Figure 6.5.   151 xi Statement of Co-authorship All technical content in this thesis were written by Ashish Uthama and edited by his supervisor, Dr. Rafeef Abugharbieh. Dr. Martin J. McKeown provided the clinical explanations relevant to the results presented in this thesis. The experimental data were collected by Samantha J. Palmer, Drew Smith, Cara Slagle, Mechelle Lewis and Xuemei Huang. All data analysis methods were designed and developed by Ashish Uthama under the supervision of his supervisor. xii 1 Thesis Introduction and Summary Magnetic resonance imaging (MRI) is the most widespread brain imaging technique in use today. In its basic form, MRI is routinely used in clinical neuroscience to non- invasively capture anatomical images of the brain for use in detection, diagnosis and tracking of neurological disorders such as Multiple sclerosis (MS), Parkinson's, (PD) and Alzheimer's disease. Active research in exploring other uses of MR imaging technology has resulted in exciting new avenues which are distinct research domains in their own right. One such area is the study of brain function using functional magnetic resonance imaging (fMRI). The mainstream MRI technique, referred to as structural or anatomical MRI, provide high resolution 3D images of the brain. This facilitates the study of (geometric) shape of cortical and subcortical brain structures, e.g. the thalamus or the ventricles. On the other hand, fMRI, which is currently one of the most actively researched MR imaging modalities, aims to identify regions of the brain involved in performing specific actions such as motor or visual tasks. Such knowledge not only increases our understanding of brain functionality, but also helps neurologists better track the progress of neurodegenerative diseases and uncover whether, and if so, how the human brain reacts to disease by employing various compensatory mechanisms [12]. 1 1.1 Thesis Motivations and Objectives One of the main challenges in MR image analysis, particularly functional imaging, is the need to meaningfully combine results from a subject pool and deducing group differences. A representation of the observed data for each subject in a common space is typically required to perform such analysis. Currently, the most widespread approach used in functional imaging is to spatially normalize the images of subject brains to a common template or atlas. The availability of easy to use software tools (e.g. The SPM toolbox [13]) tends to encourage the indiscriminant use of such an approach without due diligence to the validity of the underlying assumptions involved. Several shortcomings of this approach have been highlighted in recent studies; the major concern is the validity of normalization for pathological brains [7-10]. The core idea behind these normalization approaches, the voxel based morphometry method (VBM) was initially criticized for its susceptibility to misregistration [11], and this was later corroborated with evidence from functional studies [14, 15]. These limitations intensified the interest in alternative methods to analyze the data directly in the subject space without introducing errors due to the normalization process. Researchers have adopted simple measures of observed data such as percentage of activation, mean value etc. within specific region of interests (ROIs) to perform their analysis [16, 17]. While these methods have the advantage of being simple to compute, they clearly do not incorporate any spatial information of activation within the ROIs. More recent studies proposed a feature based approach to characterize and discriminate functional data based on spatial activation patterns [18, 19]. However, these approaches lacked a definite means to account for the structural variations in the ROI masks. 2 The primary aim of this thesis was to develop a feature based approach to sensitively discriminate functional patterns of activation in specific regions of interests despite structural intersubject variations. First, in order to characterize the structural variations, we developed a spherical harmonics based feature representation for binary ROI masks. This, in itself, is a significant contribution to ROI shape analysis area since it enabled direct comparisons of the shape across subjects without a need for mutual registration. We then extended our feature representation to functional data in order to obtain an invariant feature representation of the actual spatial distribution of activation within ROIs. Furthermore, to reduce the effect of structural variations of the ROI shape on these features, we developed a principal component subspace projection technique that further increased our method's sensitivity. 1.2 Principles of Magnetic Resonance Imaging (MRI) Magnetic resonance imaging was originally referred to as Nuclear Magnetic Resonance (NMR), based on the underlying physical phenomenon used in image capture and generation. The core of the technology relies on the fact that protons (nuclear particles) have a spin, and since protons also carry a positive charge, this results in a magnetic moment [20]. The net effect of this magnetic moment depends on the number of protons and neutrons in a nucleus, implying that different elements would have different magnetic properties. An even number of particles with spin can cancel each other out, thereby showing no outward manifestation of this effect. However, the nuclei of certain isotopes of elements like Hydrogen have unpaired particles, thus possessing a net spin. The average Hydrogen content in the human body is known to be around 63% by mass [21]. 3 This makes it an ideal target to manipulate using external magnetic fields to image the human body. A nucleus with a net magnetic moment, when placed in a magnetic field of strength B, will tend to align itself along or opposite to the field direction. Alignment along field lines is called the low energy state and alignment in the opposite direction is termed the high energy state. According to the laws of nuclear physics, a particle can switch states by either absorbing or emitting a photon with energy E equal to the difference between these two states (1.1) where h is Planck's constant (h = 6.626x10 -34 J s). This process of absorption and emission enables the MR phenomenon to be indirectly observed. E = hyB^(1.1) y is called the gyromagnetic ratio, an intrinsic property. The frequency of this photon is given by (1.2), termed the Larmor frequency. E v — h (1.2) Consider an object with thousands of such nucleus, the resonance in MRI comes from the fact that there is continuous exchange of energy between the two states of the nuclear particles constituting the object. This phenomenon is also termed spin-spin interaction. In normal conditions, random alignments of particles cancel out and result in a zero net magnetic moment. However, under the influence of a constant external field B, the number of particles in the lower energy state (parallel to the field lines) will outnumber the ones in the higher state (anti-parallel) thereby resulting in the object possessing a net 4 magnetic moment M o along a direction parallel to B. According to accepted convention, the static field B is defined to be along the Z direction in 3D Cartesian coordinate system, B. Hence, under equilibrium, the net magnetic moment is given by (1.3) MZ = Mo^( 1 . 3 ) This equilibrium state can be disrupted by the application of external energy. Depending on the energy difference between the two states (1.1) for each particle, this would cause particles to flip and hence change the direction of the net magnetic moment vector. This external source of energy is usually in the form of a radio frequency (RF) pulse, whose precise energy is carefully computed to flip the net magnetic moment by a known amount. Since the energy is pulsed, this flip is temporary, and the particles return back to the equilibrium state by dissipating energy to their surroundings. According to terminology from initial applications of MR to study crystal structures, this interaction is termed spin-lattice interaction. Depending on the energy of the pulse and intrinsic characteristics of the particles, there is a certain time constant involved before equilibrium in the direction of magnetic moment is achieved. The time constant, also called the T1 relaxation time, is defined as the time taken for the magnetic field M z to return to 63% of its equilibrium value M o (1.4) / Mz = Mo 1— e rt^(1.4) The flipped magnetic moment, also termed the transverse magnetic moment, actually rotates around the Z axis. Initially, when the external energy is applied, most particles 5 rotate in phase. After this external source is removed, individual particles being to rotate differently based on their surroundings. In addition to returning back to equilibrium state in the direction of magnetism, the rotation also begins to go out of phase. Usually, before Mo is restored, the rotation settles down from a focused, in phase rotation, to the initial random phase. The time constant involved in this relaxation of phase is termed T2. Ideally, the main static field B is assumed to be uniform throughout the substance being studied, however, in reality the field is not perfectly homogenous. This results in different parts of the object having different energy differences (1.1), and hence different Larmor frequencies (1.2), the final result is that the relaxation time T2 is in practice influenced by both particle interaction and field inhomogeneity resulting in a shorter relaxation time termed T2 * . The T2 time constant can still be measured using repeated complex RF pulse sequences. Since the magnetic moments are rotating about the Z axis, the magnitude of moments can be estimated using pick up coils placed in close vicinity. The induced current in this coil is then observed to deduce the various relaxation times. To obtain an internal image of the substance being studied, the static field B is augmented with additional fields called the gradient fields. First, a gradient field Gz is applied along the Z direction. This changes the perceived B value and correspondingly changes the Larmor frequency (1.1) along the Z axis. Now an RF pulse with energy equal to the exact Larmor frequency at a given location on the Z axis is applied. Magnetic moment of particles along the X-Y plane at that location would be selectively flipped by this pulse. This process is called slice selection since only a thin slice of the object being imaged has been excited. To further isolate signals picked by the coils to specific (x, y) locations within this slice, two more gradient fields are used as shown in Figure 1.1. 6 These gradient fields introduce a controlled phase shift in the precession process of the spinning particles. The frequency and phase are both reflected in the current signals picked by the coils. These signals are said to be in K-space since they reflect the frequency component observed. A Fourier synthesis from this space yields the 2D image of the slice being studied. This process is then repeated along the volume to obtain a 3D representation. This is one of the most basic approaches to obtain an MR volume, more complex pulse sequences, like the echo-planar imaging (EPI) [22] are used in practice. ,ig.ro.v. ,, .;';',17...a--ii.....4:4.:.:9:.:4,41.° ::"....•'...:.-...;!::-.:'''''' G Y Figure 1.1 Gradient fields along the X and Y axes. 7 1.3 Structural Brain MRI Various combinations of pulse sequences and measured signals result in a wide array of possible MR imaging protocols [23]. Currently, the most widespread set of protocols are those which image the structure of the brain [24]. Structural MRI aims to obtain high resolution images of the different types of brain tissues. Brain tissue is mainly classified into: • Grey Matter: Mainly consists of unmyelinated neurons (Myelin is a protective sheath around axons of the neurons). Grey matter concentration has been known to be positively correlated with human intelligence. • White Matter: Mainly consists of myelinated neurons. Known to be involved in message transfer across the human central nervous system • Cerebrospinal fluid (CSF): A clear isotonic fluid which occupies the space inside and around the human brain. Mainly concentrated in the brain ventricles. Structural MRI plays a significant role in the diagnosis of neurological diseases like Multiple Sclerosis [25, 26]; where specific pathological conditions like lesions appear in the white matter. The high level of detail in imaged brain anatomy enables neuroscientists to selectively study the anatomy of individual brain structures such as the thalamus, ventricles etc [2, 3, 27, 28]. The most commonly used brain structure imaging protocols are the T1 and T2 weighted images (Figure 1.2). Table 1.1 shows the relaxation times for some of the brain tissue. 8 Table 1.1 MR relaxation times for commonly occurring brain tissue Tissue Type Ti relaxation time (ms) T2 relaxation time (ms) Adipose tissue (fat) 240-250 60-80 Cerebrospinal fluid 2200-2400 500-1400 White matter 780 90 Grey matter 920 100 The difference in relaxation times for different tissues shows up as difference in intensities in the obtained MR image. As the values in the table show, Ti images have a higher contrast between the various tissues and are more suited for segmentation of brain tissue. However, certain pathological tissues, like grey matter in multiple sclerosis affected brain [291 are better contrasted in T2 weighted images. Moreover, as explained later, T2 images have certain properties that make them sensitive to the amount of oxygen in the blood, a basic concept used in functional imaging. Figure 1.2 Whole brain images obtained using different relaxation time weighted protocols showing the different tissue contrasts. Left: T i weighted image. Right: T2 weighted image. 9 1.3.1 Summary of Structural Data Acquired The structural data used in this thesis was obtained from MRI scans of 21 control (normal) subjects and 21 PD patients. Two T 1 weighted high-resolution anatomical images (voxel dimensions 1.0x1.0x1.0mm, 176x256 voxels, 160 slices) were obtained with a gap of two hours between the two scans. This second observation helps to check the consistency of our analysis methods. The MR data were acquired on a 3.0 Tesla Siemens scanner (Siemens, Erlangen, Germany) with a birdcage type standard quadrature head coil and an advanced nuclear magnetic resonance echoplanar system. Foam padding was used to limit head motion within the coil. 1.3.2 MR Image Preprocessing The raw images obtained by MRI are influenced by a variety of noise sources [30]. Basic preprocessing steps are necessary to increase the SNR, most importantly when automated analysis of the images follows. The pre-processing stage typically involves a varying number of steps depending on the nature of the study; here we explain the steps carried out in processing our data within the context of the work presented in this thesis. 1.3.2.1 Brain Extraction MR images of the brain normally include signal from the brain tissues, spinal chord, eyes, nasal cavity etc. These regions are typically not of interest and need to be excluded from the image data that is to be analyzed, moreover regions like the nasal cavity actually degrade MR signal around it. The presence of these artifacts is seen as a confounding factor for most automatic registration methods, typically employed in multi-subject studies, hence the importance of skull stripping or brain extraction [31][32]. 10 Figure 1.3 Left: A raw T1 weighted image (axial slice), notice the noise near the nasal cavity and the eyes (upper part of image). Right: The extracted brain using the combined expectation maximization and geodesic active contour method proposed by Huang et al [33]. To process the data in this thesis, we use a recently proposed fully automated method by Huang et al [33] which was shown to provide quantitative advantages over existing methods. This technique begins by first pre-processing MRI scans to generate inputs to a deformable model — a geodesic active contour [70]. After the model converges to a stable solution, the mask is post-processed to further refine the perimeter of the brain cortex. The pre-processing and post-processing stages employ the expectation-maximization (EM) algorithm on a mixture of Gaussian models as well as mathematical morphology and connected component analysis. Figure 1.3 shows the result of using this technique on one of our subject's T 1 weighted brain MR volume. 1.3.2.2 Anisotropic Diffusion Filtering Most brain tissue and shape analysis techniques are automated and require a high SNR for effective results. One of the most common preprocessing steps is in reducing the impact of noise using spatial smoothing. Normal Gaussian kernel based spatial filters are 11 not well suited for smoothing structural MR data since the isotropic nature of the kernels would affect the signal integrity near tissue boundaries. As explained later, these boundaries guide automated segmentation of brain tissues and regions and thus need to be preserved. A form of antistrophic smoothing proposed by Perona et al [34] is widely used in MR image processing to preserve the edges while improving the SNR. Anisotropic filters have spatially varying diffusion coefficients which allow them to selectively perform intra-region smoothing rather than inter-region smoothing, thereby reducing the degradation of image boundaries. We used an ITK [35] implementation of a 3D version of this popular filter to process our structural data. Figure 1.4 shows the effect of using this filter on a noisy T 1 brain image. Figure 1.4 Result of anisotropic diffusion filtering on a Ti brain volume (one slice shown). 1.3.3 Segmentation of Regions of Interest in the Brain Segmentation of brain areas involves the process of localizing specific regions within the brain, either based on tissue class or based on the constituting brain structure. Segmentation methods can be broadly classified into manual, semi-automatic and 12 automatic segmentation techniques [36, 37]. Accurate segmentation of specific regions in the brain is very important for region of interest (ROI) based analysis approaches. In the study of shape or structure of brain regions like the thalamus for example, errors in segmentation can have a very significant impact on the analysis and subsequent medical inference. ROI segmentations are almost always done on high resolution anatomical scans, like the T 1 MR volumes. As explained in later sections, functional studies, which involve comparatively lower resolution scans make provision to acquire a high resolution anatomical (usually T i ) volume to enable segmentation of regions of interest. 1.3.3.1 Manual Segmentation of the Thalamus Most of the ROIs obtained for the data used in this thesis were segmented using manual techniques. Considerable experience is required to identify and delineate some of the ROIs in MR images (Figure 1.13). Even though this is a very time consuming process, manual segmentation by a trained expert is considered a gold standard in segmentation applications. A number of software packages exist to help trained neuroscientists accurately delineate regions of the brain in MR images. ROIs involved in the functional study (explained in section 1.4) were segmented manually using the Amira software package (Mercury Computer Systems, San Diego, USA) on the collected Ti images. Two of these ROIs, the left and the right thalamus were selected for structural analysis since it was shown that they undergo atrophy in Parkinson's disease patients [38]. 1.3.3.2 Automatic Segmentation of Brain Ventricles Manual segmentation processes are frequently limited by the availability of trained experts and the time involved in the process. Automatic segmentation methods are 13 preferred whenever the region of interest has a high SNR and when its accuracy and robustness have been thoroughly tested. We used the voxel based morphometry (VBM) [39] approach to segment the left and right brain ventricles from the T 1 brain volume scans. Despite its known limitation in aligning cortical features of the brain [11], it was expected to perform better for the ventricles due to their high-contrast and relatively more regular boundaries compared to the cortical gyri and sulci [40]. First, the left and the right ventricles were individually segmented on the tissue probability map (TPM) for the cerebrospinal fluid (CSF), then, using the VBM normalization approach, these maps were warped on to the subject brain. The areas in the subject brain corresponding to the segmented area on the TPM were labeled as the ventricles. A connected component analysis was then used to clean out the segmented ventricle of small disconnected voxels. Figure 1.5 displays the result of one such segmentation. Figure 1.5 Result of the automatic segmentation of the right brain ventricle. The left ventricle's binary ROI mask is rendered on orthogonal slices of the subjects T i image. 14 1.3.4 Quantifying Structural Changes in Brain ROIs It is now well known that the brain consists of different regions each associated with different functionality. It is also known that some neurological diseases affect certain parts of the brain; this has encouraged targeted analysis of anatomical regions in the hope of learning more about these diseases [5, 27, 28]. A region of interest (ROI) is a specific part of the brain delineated in 3D space either using automatic segmentation or manual methods for such a targeted study. While analyzing the structural features of an ROI, the data being analyzed is the binary ROI mask of the region in question. The simplest and most common approach to quantify the binary ROI mask is to use the volume or the voxel count (number of voxels contained in the ROI), e.g. Anstey et al used the hippocampal volume to investigate changes observed in mild cognitive impairment (MCI) [41]. However, studies using volumetric measures do not reflect shape changes that might be indicative of neurological change. Vetsa et al. used medial representations to characterize the surface of the caudate and identify shape changes in schizophrenia [42]. This approach requires accurate alignment of the ROIs across different subjects for valid analysis and hence residual errors in mutual alignment could adversely affect the final results. Alternate methods to quantify this mask rely on obtaining descriptive feature vectors which reflect the structural information in an ROI. These features may be invariant to rotations, translations and scale, thus enabling easy comparison across different subjects. This eliminates the need for mutual alignment, saving time and potential inconsistencies arising from registration errors. Zimring et al [43] used spherical harmonic (SPHARM) 15 based invariant indices to accurately estimate and study volume of brain lesions, reducing the errors due to subject positioning. Tootoonian et al [28] used scale, rotation and translation invariant features derived using a SPHARM representation of 3D ROIs to analyze shape changes in brain MRI data and demonstrated the advantages of using such an approach along with traditional volumetric analysis. However, the parameterization technique used there did not allow for shapes with non-convex topology. Gerig et al [44] proposed using a SPHARM surface representation for brain structures like the ventricles to better understand neurodevelopment and neurodegenerative changes. Their approach however requires explicit alignment of the structures along with volume normalization and is valid only for structures with spherical topology. 1.3.4.1 SPHARM features for Quantifying Structure in Brain ROIs This thesis presents a novel way to quantify 3D ROIs obtained from brain MR images using invariant spherical harmonics based features. Unlike previous implementations [68,69], this approach places no restriction on the allowable topology of the ROI masks while ensuring that the obtained features are unique to its shape. These features are translation, scale and rotational invariant which allow them to be directly compared between subject groups without a need for mutual registration. The sensitivity of these features was validated on synthetic data with realistic intersubject variability. Furthermore, application of this method to the ROI masks of the left and right thalami of real PD patient data revealed systematic difference in shape when compared to normal subjects [27]. Section 2 describes the method in detail including derivations and synthetic tests along with preliminary clinical results. Section 3 and Section 4 build on the results and presents the clinical relevance of our findings. 16 1.4 Functional Magnetic Resonance Imaging Functional magnetic resonance imaging (fMRI) is currently one of the most actively researched branches of the MR technique. Using MR pulse sequences which can detect change in the blood Oxygen content, this imaging modality can be used to localize regions in the brain responding to specific stimuli. Neuroscientists use this technique on normal human subjects to obtain functional maps of the brain describing the role of various regions in processing common cognitive tasks. These maps enable them to compare the functional responses of pathological brains, enabling a detailed study of the effects of neurological damage on functional response. Recently, using this technique, researchers have been able to show plastic reorganization of neurocognitive networks which would otherwise not be detected using outward manifestation of symptoms [12, 45]. 1.4.1 Principles of BOLD imaging The seminal work by Ogawa et al [46] laid the foundation for the study of function using MR imaging. It was noticed that oxygenated and deoxygenated blood have distinct T2 relaxation times due to their different magnetic properties (Table 1.2). In practice, due to magnetic field in-homogeneities, the observed relaxation time is much faster and is called the T2 * signal. 17 Table 1.2 MR relaxation time for Oxygenated and Deoxygenated blood Blood Type T1 relaxation time (ms) T2 relaxation time (ms) Oxygenated blood 1350 50 Deoxygenated blood 1350 200 This phenomenon, however, does not enable precise computation of the amount of Oxygen in the blood. Instead, neuroscientists use this to measure the change in blood Oxygen concentrations, hence the name blood oxygen level dependent (BOLD) imaging. 1.4.2 Hemodynamic Response The change in blood Oxygen content in response to activation in neurons is called the Hemodynamic response. Neurons are the computing cells of the brain, a basic unit which processes the various stimuli that the brain receives. In response to a specific stimulus, neurons associated with the region that processes this information are activated. Active neurons require energy derived from bio-chemical reactions which consume Oxygen. Hence an increase in activation levels requires a larger amount of Oxygen when compared to the rest state. This increase in Oxygen requirement triggers a Hemodynamic response which increases the blood flow into these regions. However, this process overcompensates and results in an increased oxygenated blood flow into the active regions [47]. From an imaging point of view, if this region was monitored before the stimulus was presented and after the Hemodynamic response was triggered, a change in the T2 * signal would be observed. In practice, this change is of a very small magnitude and hence hard to detect. Neuroscientists resort to the basic technique of obtaining 18 repeated observations to obtain a better SNR. A comprehensive review by Heeger et al outlines this process is more detail [48]. 1.4.3 Functional Experiment Setup The main challenge in functional imaging is to overcome the low SNR observed in the BOLD signal. Given the current limitations of MR scanners in terms of image resolution, SNR and the acquisition time, the most robust fMRI experiment used is the block design approach. This is also the simplest and easy to analyze technique unlike other more complex schemes like event based experiments. In the single task block design, a subject performs a given task repeatedly, in blocks of time (duration), while MR brain images are being continuously acquired. The tasks are interspersed with rest periods to let the Hemodynamic response settle down to baseline levels. The data so obtained is also called the time-series or the 4D functional data. Given this structure of a functional experiment, the volumes (3D brain images) need to be acquired with the sufficiently small acquisition time to accurately sample the BOLD response. This temporal resolution requirement limits the spatial resolution of fMRI volumes; often the functional volumes acquired have a very low spatial resolution compared to the structural images obtained (which have a longer acquisition time). Usually, a high resolution structural scan is separately obtained before or after the functional experiment. Functional activations observed in the fMRI data can then be examined for correspondence to physical structures using the anatomical scans. ROIs drawn on the anatomical scans can be mapped to the functional data to enable ROI based analysis of the time series. The precise nature of the applied stimulus or task performed varies with each fMRI study [49]. 19 1.4.3.1 Bulb Squeezing Task Experiment In neurological diseases like Parkinson's disease (PD), one of the primary external manifestations of the disease is motor task related disabilities. Hence, a motor task involving the squeezing of a bulb was used in our study. Subjects lay on their back in the functional magnetic resonance scanner viewing a computer screen via a projection-mirror system. All subjects used an in-house designed response device in their right hand, a custom-built MR-compatible rubber squeeze-bulb (Figure 1.6) connected to a pressure transducer outside the scanner room. Figure 1.6 The MR-compatible rubber squeeze bulb used in the functional experiment The subjects lay with their forearm resting down in a stable position, and were instructed to squeeze the bulb using an isometric hand grip and to keep their grip constant throughout the study. Using the squeeze bulb, subjects were required to control the width of an "inflatable ring" (shown as a black horizontal bar on the screen) in order to keep the ring within an undulating pathway without scraping the sides (Figure 1.7). 20 Figure 1.7 Visual feedback given to subject performing the bulb squeezing task. The width of the black horizontal bar reflected the amount of 'squeezing' that the bulb recorded. Subjects had to modulate the squeezing action to keep the width within the undulating white pathway. The pathway used a block design with sinusoidal sections in two different frequencies (0.25 and 0.75 Hz) in a pseudo-random order, and straight parts in between where the subjects had to keep a constant force (Figure 1.8). STIMULI DESIGN FOR FMRI TASK 302, 1500 ^ 2000 SAMPLE POINT Figure 1.8 The stimulus design used for the functional study. The tunnel width shown on the Y axis indicates the width of the undulating pathway shown in Figure 1.7. The frequencies were chosen based on prior findings and pilot studies were used to determine that PD subjects could comfortably perform the required task. Each block lasted 20 seconds, alternating a sinusoid, constant force, sinusoid and so on to a total of 4 21 minutes. Before the first scanning session, subjects practiced the task at each frequency until errors stabilized and they were familiar with the task requirements. Custom Matlab software (The Mathworks, Natick, MA) and the Psychtoolbox [50] were used to design the stimulus and present the visual feedback. 1.4.4 Summary of Functional Data Acquired Functional MRI for the experiment described in the section 1.4.3.1 was acquired on a Philips Achieva 3.0 T scanner (Philips, Best, the Netherlands) equipped with a head-coil. Echo-planar (EPI) T2 * weighted images with blood oxygenation level-dependent (BOLD) contrast were acquired. The 2D slice acquired had a matrix size of 64 x 64 with pixel size 3 x 3 mm. Thirty-six axial slices of 3mm thickness were collected in each volume, with a gap thickness of lmm. Slices were selected to cover the dorsal surface of the brain and include the cerebellum ventrally. A high resolution, 3-dimensional Ti-weighted image consisting of 170 axial slices was also acquired of the whole brain to facilitate anatomical localization of activation for each subject. Figure 1.9 shows a sample T1 slice against a time sequence of T2 * . 22 Figure 1.9 Left: A slice from a T i weighted anatomical scan of a subject. Right: Corresponding slices of T2 * scans from functional volumes acquried the duration of the experiment (130 such volumes were collected over a duration of 260 seconds). 1.4.5 fMRI Data Preprocessing The raw time course of the fMRI data has to be put through a complex image processing pipeline before the final results can be obtained. The following is a brief introduction to the various steps in this pipeline. While it includes most of the commonly used steps, it must be noted that there is considerable variation in the methods themselves and the order in which they are used [51]. Most fMRI studies collect the functional time-series as well as at least one high resolution anatomical image. The task performed and the image acquisition parameters are assumed to be the same for all subjects involved. While some studies are limited to a single subject, especially in cases of rare neurological conditions, most studies image a number of subjects. These subjects are usually age and sex matched between two different populations under study (for example, PD patients and normal subjects). The goal of the analysis then is to use the data from individual subjects and deduce the global group-wise difference in activation response. 23 1.4.5.1 Brain Extraction and Co-registration The first step is similar to the preprocessing involved in structural scans (Section 1.3.2.1). The brain is extracted from the raw MR image individually. The structural scan is then aligned to the first volume of the functional time series to compensate for any subject movement between the two scans. This registration is rigid; involving translation and rotation about the three axes. This process helps neuroscientists to identify regions in the brain which would otherwise have been difficult using the low spatial resolution fMRI images. Brain Extraction and Co-registration Figure 1.10 The process of brain extraction, followed be co-registration of the high resolution T 1 anatomical scan to the first functional T2 * voume 1.4.5.2 Slice Timing Correction for Individual Functional Volumes As explained earlier, the 3D volume acquisition of the brain proceeds with the sequential acquisition of single 2D slices. This implies that each slice of a volume is acquired with a definite time lag with respect to the previous slice. Slice timing correction is important in functional analysis since slice acquisition differences appear as a phase shift in the time- series of individual voxels. Since the BOLD signal being analyzed is temporal in nature, this artifact can confound further analysis. Slice timing correction procedures attempts to 24 correct for this time difference by correcting the phase of each acquired slice in the Fourier domain without affecting the magnitude of the observed signal. Slice timing correction was applied using Sinc interpolation in Brain Voyager (Brain Innovations, the Netherlands, www.brainvoyager.com) for the data collected for analysis in this thesis. 1.4.5.3 Motion Correction of Functional Volumes A complete fMRI study scan requires the subject to be inside the MRI scanner for up to an hour depending on the experiment design. Despite head rests and motion arrestors in the machine, it is inevitable that there is some amount of displacement in the subject's position during the course of the scan. Since the fMRI data analysis depends on the time series analysis of a voxel, there is an inherent assumption that a specific voxel always corresponds to the same physical volume throughout the session. Head motion correction steps attempt to rectify this situation by realigning individual scans. Typically, in cases of stronger forms of stimuli, subjects tend to move in synch with the change in stimulus (stimulus correlated motion) [52]. This causes a periodic change in the time course of a voxel which is in phase with the stimulus. These types of artifacts will have an adverse effect on the analysis and may as well render the entire data unfit for analysis. Researchers usually analyze the data sets manually and screen out data sets which show large head motions. A common practice is to align all the functional volumes to the first acquired volume. This is achieved using rigid registration processes, usually based on the intensity metric [53]. Accurate motion correction is extremely important, since task correlated head motion (common in motor tasks) can have adverse influence on the functional analysis. 25 There a number of practical approaches to correct for subject motion [54]. The data used in this thesis was processed with a motion corrected independent component analysis (MCICA) method proposed by Liao et al [55]. Motion Correction Figure 1.11 Motion correction re-aligns the functional volumes in case the subject moves during the time course of acquisition. 1.4.5.4 3D Spatial Smoothing of the Functional Volumes Spatial smoothing involves using a 3D kernel to remove high frequency component in the individual fMRI volumes [56]. Spatial smoothing is primarily done to increase the SNR. Smoothing is also assumed to reduce the effect of any residual error in the motion correction and registration process, by increasing the effective overlap between assumed homologous regions across subject brains. One branch of fMRI analysis also requires an estimate of the smoothness in the data for multiple comparison correction [57]. Often, since the true smoothness cannot be estimated, it is taken to be equal to the amount of smoothing performed. Given that fMRI resolution is low to begin with, smoothing can lead to signal degradation and reduced sensitivity in later analysis steps [58, 59]. 26 The functional data used in this thesis was not smoothed to enable the derived features to be more sensitive to spatial locations of activation statistics. The feature based characterization approach described in this thesis does not require an estimate of smoothness of the data; neither does it employ registration to warp the functional data. Moreover, the sensitivity of the proposed features despite high SNR further validates using the functional data as is, without any smoothing. 1.4.6 Functional Activation Map Generation for Individual Subjects The first post-processing step in fMRI analysis is to individually analyze a subject's functional response. The goal of this step is to obtain a measure of activation in each voxel of the subject's brain. There are numerous methods to obtain this measure, with a common approach being to study the time course of each voxel and compare it against an expected response in order to obtain a parametric measure of confidence in their similarity. A major hurdle is in defining the expected response. One of the most simple and widely used methods is to model the Hemodynamic response as a Volterra kernel [60] function and convolve this with a boxcar function representing the task-rest periods. Various approaches are then defined to obtain an activation metric ranging from simple correlation [61] to general linear models [13]. One drawback of these approaches is in the assumption that the Hemodynamic response is uniform in time and space. In this thesis, we use a more sophisticated approach based on independent component analysis (ICA) to extract the components involved in the observed time course. McKeown et al [62] showed that this approach could efficiently describe the response and isolate effects of non-task related signals, movements and other artifacts. The end 27 Task performed Observed Time course of i'th voxel result of a subject level analysis is what is called a statistical parameter map (SPM). This SPM is a 3D volume with the value at each voxel position representing a quantitative measure of confidence of that location being active during the task. Obtain parameter for i'th voxel Compute Measure of Activation for each voxel Assemble parameters in 3D to obtain subject SPM Figure 1.12 The process of obtaining subject level statistical parameter maps. The time course of each voxel is statistically compared to an expected response based on the task performed to obtain a parameter which describes the activation at that voxel. This process is repeated for all voxels in the brain resulting in the subject level SPM. 28 1.4.7 Differences in Activation Maps at the Group Level Once the subject level parameter maps have been obtained, the next step in functional analysis is to combine maps within each group and compare them against each other. This is a key step in functional analysis and there has been quite a few competing methods proposed. A main contribution of the thesis is in this step of functional analysis. By far the most popular approach is to warp each subject's brain to a common atlas, thereby spatially normalizing it [63]. Then the difference in means of the parameter at each voxel location between the two groups is statistically tested for significant difference. A final group level parameter map is then generated using the results of this statistical analysis. Further inferences are drawn directly from this map. This approach heavily relies on the non-rigid registration processes to perform the warping of the brain. A comparison of some of these methods was presented by Crivello et al. [64]. However, given the large amount of intersubject variations in the human brain, current spatial normalization methods may give an imperfect registration result [14]. This might result in signals from functionally distinct areas, especially small subcortical structures, to be inappropriately combined [15]. Spatial normalization may therefore lead to poor sensitivity in the fMRI data analysis due to reduced functional overlap across subjects, as observed by Stark and Okado [65]. An alternative approach more recently investigated is to align the subjects at the region of interest (ROI) level, as opposed to whole brain normalization. These studies are usually hypothesis driven and examine signals from specific parts of the brain. Miller et al. [66] extended the work of Stark and Okado [65] by utilizing large deformation diffeomorphic 29 metric mappings (LDDMM). Also, an approach that uses continuous medial representations (cm-rep) of the ROIs has been proposed [67]. However, these approaches rely on structural features and manually placed landmarks in the higher resolution scans to deform the functional maps into a common space. To circumvent these potential problems associated with normalization, many researchers employ feature based ROI analysis without using any normalization [16]. This results in enhanced sensitivity in activation detection in some cases [17]. However, to date, most methods ignore information encoded by the spatial distribution of the activation statistics and instead use simpler invariant (to pose) measures such as mean voxel statistics or percentage of active voxels within an ROI as features [16]. An attempt in incorporating spatial information was proposed in [19] using sums of activation statistics within spheres of increasing radii. However, these features have limited sensitivity to spatial patterns since they only capture changes in the radial direction. Ng et al. [18] proposed an alternate approach employing three dimensional moment invariant features (3DMI) which were shown to be more sensitive than conventional mean voxel statistics and percentage of activated voxels based approaches in detecting task-related activation differences. While a novel and interesting approach by itself, the method does not account for intersubject variability present in the ROI masks or its effect on the analysis. Also, as the authors note, higher order 3DMI features are more susceptible to noise and hence limit the maximum number of discriminative features that can be derived [18]. 30 SPHARM-based methods have previously been proposed in the context of 3D shape retrieval systems by Vranic et al [68] and Kazhdan et al [69]. Both authors proposed obtaining invariant SPHARM features by intersecting 3D shapes with shells of growing radii. This approach, however, could not detect independent rotations of a shape along the shells, thereby resulting in a non-unique representation [69]. In [28], the use of invariant SPHARM descriptors for analyzing anatomical structures (represented as binary 3D shapes) in MR was proposed. However, the shape parameterization employed was limited to convex topologies, which significantly limits its applicability. The three approaches described above were used to analyze the geometrical shape of an ROI surface, not the intensity distributions of an ROI volume, such as the spatial activation patterns in fMRI data. 1.4.7.1 SPHARM Features for Analyzing Spatial Activation Patterns This thesis proposes the use of SPHARM features to characterize spatial patterns of activation in functional MR data for the first time. These features are invariant to scale, translation and rotation and hence can be directly compared across subjects, thus avoiding the various problems associated with the normalization approach. The affect of anatomical variability between subjects on these functional features is accounted for by using a novel principal component subspace approach. Synthetic data experiments constructed with real brain ROIs and synthetically created spatial activation patterns showed significantly improved performance of the SPHARM approach against one of the latest spatial normalization approaches. Furthermore, when applied to real data, SPHARM features were able to detect changes in regions which were unnoticed using the 31 conventional normalization approach. Some of the ROIs analyzed in the real data set are shown in Figure 1.13. Figure 1.13 The functional ROIs studied in this thesis. All ROIs shown occur in pairs (in the left and righ brain hemishpheres), only regions in left hemisphere are shown. Indicated ROIs are Prefrontal cortex (PMC), Supplementary motor cortex (SMA), Primary motor cortex (M1), Caudate (CAU), Putamen (PUT), Thalamus (THA), Cerebellar hemisphere (CER). 1.5 Thesis Contributions This thesis makes three important contributions to region based analysis of MR images. The first contribution enables generating invariant SPHARM features to describe 3D functions uniquely. The second contribution involved extending this method to analyze spatial patterns of activation in functional data. Finally, the principal component subspace approach presented reduces the impact of anatomical variations in the functional features. 32 1.5.1 Structural ROI Analysis An automated feature based discrimination method which complements traditional volumetric analysis of anatomical ROIs using structural features was presented. These SPHARM features are invariant to translation, rotation and scale, thereby eliminating the need for mutual registration or spatial normalization as required by other methods. Moreover, unlike other methods, there is no need for manual intervention once the ROIs are obtained. This reduces any chance of bias or variability in landmark placement which might affect the results of mutual registration dependent methods. This approach was used to analyze real data and was able to discriminate shape changes in the thalami [27] and ventricles of PD patients. Group representation using shapes closest to group-mean vectors enabled neuroscientists to further study the nature of the change observed. With increasingly accurate and robust automatic segmentation algorithms being developed, this feature based approach has the potential to be completely automated once the raw MRI image is obtained. In summary, the contributions of this thesis to ROI based shape analysis are: • Addressed previous limitations of SPHARM based invariant feature representation using a radial transform to obtain unique features, thereby extending its use to more complex shapes including those with possible spatial disconnects. • Use of this SPHARM based unique invariant features for the analysis of shape in anatomical regions of interest. 33 1.5.2 Functional ROI Analysis Spatial patterns of activation have been shown to be more sensitive than traditional mean or percentage activation in the analysis of fMRI data [18]. In this thesis we presented a novel scheme to quantify and analyze spatial functional patterns using SPHARM features. For the first time, we proposed the use of a principal component subspace approach to quantify and reduce the influence of intersubject anatomical variability in functional pattern analysis using feature based representation. Using projections on this reduced subspace we demonstrated improved sensitivity by minimizing confounding effects of inter-subject anatomical variations seen in ROI binary masks. A significant improvement over a recently proposed ROI normalization approach was shown using synthetic data with real world intersubject anatomical variability modeled. The proposed method was then applied to real fMRI data collected from PD patients and normal subjects. In addition to the ROIs detected by the normalization approach, SPHARM features detected changes in activation pattern of clinically relevant ROIs. In summary, for the analysis of functional activation patterns in regions of interest in fMRI data, the contributions of this thesis are: • Invariant SPHARM based feature representation for ROI based fMRI data analysis • Principal component subspace approach to reduce the impact of intersubject structural variability in the functional analysis 34 1.6 Thesis Overview Sections in this thesis are based on the following manuscripts: • Section 2: Invariant SPHARM features for anatomical ROIs are introduced. Contributions in addressing previous limitations are presented along with a clinical application for shape analysis of the thalamus in PD patients. o Ashish Uthama, Rafeef Abugharbieh, Anthony Traboulsee and Martin J. McKeown , "Invariant SPHARM shape descriptors for complex geometry MR region of interest analysis in MR region of interest analysis," Engineering in Medicine and Biology, pp. 1322 - 1325, 2007. • Section 3: Clinical relevance of the SPHARM features are presented with shape changes observed in the PD thalamus discussed in greater detail. o Martin J. McKeown, Ashish Uthama, Rafeef Abugharbieh, Samantha Palmer, Mechelle Lewis and Xuemei Huang, "Shape (But Not Volume) Changes in the Thalami in Parkinson Disease", Submitted to BioMed Central Neuroscience, 2007. o Martin J. McKeown, Ashish Uthama, Rafeef Abugharbieh, Samantha Palmer, Mechelle Lewis and Xuemei Huang, "MRI in Parkinson's Disease identifies shape, but not volume changes in the thalamus", Abstract, Selected for poster presentation at XVII WFN World Congress on Parkinson's Disease and Related Disorders, Amsterdam, The Netherlands, 2007 35 • Section 4: Invariant SPHARM shape features were applied to the study of shape change in lateral brain ventricles in PD patients. o Ashish Uthama, Rafeef Abugharbieh and Martin J. McKeown, "Invariant SPHARM shape descriptors for the analysis of brain ventricles", manuscript being prepared for submission. • Section 5: SPHARM features are extended to characterize 3D spatial activation patterns in functional MR imaging. A novel principal component subspace approach to reduce influence of anatomical variability is presented. o Ashish Uthama, Rafeef Abugharbieh, Anthony Traboulsee, Martin J. McKeown, "3D Regional fMRI Activation Analysis Using SPHARM-Based Invariant Spatial Features and Intersubject Anatomical Variability Minimization", manuscript being prepared for submission. • Section 6: Functional MR SPHARM features are compared against a competing, state-of-art normalization approach. o Ashish Uthama, Rafeef Abugharbieh, Martin J. McKeown, Samantha J. Palmer, Anthony Traboulsee and Martin J. McKeown, "Characterizing fMRI Activations in ROIs while Minimizing Effects of Intersubject Anatomical Variability", manuscript being prepared for submission. • Section 7: Conclusion and future work. 36 1.7 References [1] D. Mamah, L. Wang, D. Barch, G. A. de Erausquin, M. Gado and J. G. Csernansky, "Structural analysis of the basal ganglia in schizophrenia," Schizophrenia. Research, vol. 89, pp. 59-71, 2007. [2] D. Jong, J. Kim, T. Chung, S. K. An, Y. C. Jung, J. Lee, J. Lee, I. Kim and S. I. Kim, "Shape deformation of the insula in schizophrenia," Neuroimage, vol. 32, pp. 220-227, 2006. [3] J. J. Levitt, C. F. Westin, P. G. Nestor, R. S. Estepar, C. C. Dickey, M. M. Voglmaier, L. J. Seidman, R. Kikinis, F. A. Jolesz, R. W. McCarley and M. E. Shenton, "Shape of caudate nucleus and its cognitive correlates in neuroleptic-naive schizotypal personality disorder," Biological Psychiatry, vol. 55, pp. 177-184, 2004. [4] L. Shen, "A surface-based approach for classification of 3D neuroanatomic structures," Intelligent Data Analysis, vol. 8, pp. 519-542, 2004. [5] M. Styner, J. A. Lieberman, D. Pantazis and G. Gerig, "Boundary and medial shape analysis of the hippocampus in schizophrenia," Medical Image Analysis, vol. 8, pp. 197- 203, 2004. [6] Bermel. Selective caudate atrophy in multiple sclerosis: A 3D MRI parcellation study. Neuroreport 14(3), pp. 335, 2003 37 [7] S. Krishnan, M. J. Slavin, T. T. Tran, P. M. Doraiswamy and J. R. Petrella, "Accuracy of spatial normalization of the hippocampus: implications for fMRI research in memory disorders," Neuroimage, vol. 31, pp. 560-571, 2006. [8] H. Juch, I. Zimine, M. Seghier, F. Lazeyras and J. D. Fasel, "Anatomical variability of the lateral frontal lobe surface: implication for intersubject variability in language neuroimaging," Neurolmage, vol. 24, pp. 504-514, 2005. [9] M. W. G. Vandenbroucke, R. Goekoop, E. J. J. Duschek, J. C. Netelenbos, J. P. A. Kuijer, F. Barkhof, P. Scheltens and S. A. R. B. Rombouts, "Interindividual differences of medial temporal lobe activation during encoding in an elderly population studied by fMRI," Neurolmage, vol. 21, pp. 173-180, 2004. [10] K. Swallow, T. Braver, A. Snyder, N. Speer and J. Zacks, "Reliability of functional localization using fMRI," Neurolmage, vol. 20, pp. 1561-1577, 2003. [11] F. L. Bookstein, ""Voxel-based morphometry" should not be used with imperfectly registered images," Neuroimage, vol. 14, pp. 1454-1462, 2001. [12] R. Mallol, A. BarrOs-Loscertales, M. Lopez, V. Belloch, M. Parcet and C. Avila, "Compensatory cortical mechanisms in Parkinson's disease evidenced with fMRI during the performance of pre-learned sequential movements," Brain Research, vol. 1147, pp. 265-271, 2007. [13] K. J. Friston, A. P. Holmes, K. J. Worsley, J. B. Poline, C. Frith and R. S. J. Frackowiak, "Statistical Parametric Maps in Functional Imaging: A General Linear Approach," Human Brain Mapping, vol. 2, pp. 189-210, 1995. 38 [14] 0. T. Carmichael, H. A. Aizenstein, S. W. Davis, J. T. Becker, P. M. Thompson, C. C. Meltzer and Y. Liu, "Atlas-based hippocampus segmentation in Alzheimer's disease and mild cognitive impairment," Neuroimage, vol. 27, pp. 979-990, 2005. [15] M. Ozcan, U. Baumgartner, G. Vucurevic, P. Stoeter and R. Treede, "Spatial resolution of fMRI in the human parasylvian cortex: Comparison of somatosensory and auditory activation," Neuroimage, vol. 25, pp. 877-887, 2005. [16] R. T. Constable, P. Skudlarski, E. Mencl, K. R. Pugh, R. K. Fulbright, C. Lacadie, S. E. Shaywitz and B. A. Shaywitz, "Quantifying and Comparing Region-of-Interest Activation Patterns in Functional Brain MR Imaging: Methodology Considerations," Magnetic Resonance Imaging, vol. 16, pp. 289-300, 1998. [17] A. Nieto-Castanon, S. S. Ghosh, J. A. Tourville and F. H. Guenther, "Region of interest based analysis of functional imaging data," Neuroimage, vol. 19, pp. 1303-1316, 2003. [18] B. Ng, R. Abugharbieh, X. Huang and M. J. McKeown, "Characterizing fMRI activations within regions of interest (ROIs) using 3D moment invariants," in Computer vision and Pattern Recognition Workshop, pp 63, 2006, [19] D. Kontos and V. Megalooikonomou, "Fast and effective characterization for classification and similarity searches of 2D and 3D spatial region data," Pattern Recognition, vol. 38, pp. 1831-1846, 2005. [20] K. Tanaka, "Basic principles of magnetic resonance imaging," Rinsho Byori, vol. 48, pp. 614-620, 2000. 39 [21] M. A. Foster, Magnetic Resonance in Medicine and Biology. Pergamon Press, 1984. [22] M. Poustchi-Amin, S. A. Mirowitz, J. J. Brown, R. C. McKinstry and T. Li, "Principles and applications of echo-planar imaging: a review for the general radiologist," Radiographics, vol. 21, pp. 767-779, 2001. [23] E. F. Jackson, L. E. Ginsberg, D. F. Schomer and N. E. Leeds, "A review of MRI pulse sequences and techniques in neuroimaging," Surgical Neurology, vol. 47, pp. 185- 199, 1997. [24] M. Symms, H. R. Jager, K. Schmierer and T. A. Yousry, "A review of structural magnetic resonance neuroimaging," Journal of Neurology Neurosurgery and Psychiatry., vol. 75, pp. 1235-1244, 2004. [25] A. Traboulsee, "MRI relapses have significant pathologic and clinical implications in multiple sclerosis," Journal of Neurological Science, vol. 256 Supplement 1, pp. S19- 22, 2007. [26] R. Bakshi, "Magnetic resonance imaging advances in multiple sclerosis," Journal of Neuroimaging, vol. 15, pp. 5S-9S, 2005. [27] A. Uthama, R. Abhugharbieh, A. Traboulsee and M. J. McKeown, "Invariant SPHARM shape descriptors for complex geometry MR region of interest analysis in MR region of interest analysis," Engineering in Medicine and Biology, pp. 1322-1325, 2007 40 [28] S. Tootoonian, R. Abugharbieh, X. Huang and M. J. McKeown, "Shape vs. volume: Invariant shape descriptors for 3D region of interest characterization in MRI," in International Symposium on Biomedical Imaging, pp. 754 -757, 2006. [29] R. Bakshi, R. H. B. Benedict, R. A. Bermel, S. D. Caruthers, S. R. Puli, C. W. Tjoa, A. J. Fabiano and L. Jacobs, "T2 Hypointensity in the Deep Gray Matter of Patients With Multiple Sclerosis: A Quantitative Magnetic Resonance Imaging Study," Archives of. Neurology, vol. 59, pp. 62-68, 2002. [30] A. Macovski, "Noise in MRI," Magnetic Resonance Medicine, vol. 36, pp. 494-497, 1996. [31] K. Boesen, K. Rehm, K. Schaper, S. Stoltzner, R. Woods, E. Luders and D. Rottenberg, "Quantitative comparison of four brain extraction algorithms," Neuroimage, vol. 22, pp. 1255-1261, 2004. [32] J. M. Lee, U. Yoon, S. H. Nam, J. H. Kim, I. Y. Kim and S. I. Kim, "Evaluation of automated and semi-automated skull-stripping algorithms using similarity index and segmentation error," Computational and Biological Medicine, vol. 33, pp. 495-507, 2003. [33] A. Huang, R. Abugharbieh, R. Tam and A. Traboulsee, "MRI Brain Extraction with Combined Expectation Maximization and Geodesic Active Contours," Signal Processing and Information Technology, 2006 IEEE International Symposium on, pp. 107- 111, 2006. [34] P. Perona and J. Malik, "Scale -space and edge detection using anisotropic diffusion," IEEE Trans. Pattern Analysis and Machine Intelligence, vol. 12, pp. 629 -639, 1990. 41 [35] L. Ibanez, W. Schroeder, L. Ng and J. Cates, "The ITK Software Guide," [36] L. P. Clarke, R. P. Velthuizen, M. A. Camacho, J. J. Heine, M. Vaidyanathan, L. O. Hall, R. W. Thatcher and M. L. Silbiger, "MRI segmentation: methods and applications," Magnetic Resonance Imaging, vol. 13, pp. 343-368, 1995. [37] D. L. Pham, C. Xu and J. L. Prince, "Current methods in medical image segmentation," Annual Review of Biomedical Engineering, vol. 2, pp. 315-337, 2000. [38] J. M. Henderson, K. Carpenter, H. Cartwright and G. M. Halliday, "Loss of thalamic intralaminar nuclei in progressive supranuclear palsy and Parkinson's disease: clinical and therapeutic implications," Brain, vol. 123 ( Pt 7), pp. 1410-1421, 2000. [39] J. Ashburner and K. J. Friston, "Voxel-based morphometry-the methods," Neuroimage, vol. 11, pp. 805-821, 2000. [40] J. Ashburner and K. J. Friston, "Comments and controversies. Why voxel-based morphometry should be used," Neuroimage, vol. 14, pp. 1238-1243, 2001. [41] K. J. Anstey and J. J. Mailer, "The role of volumetric MRI in understanding mild cognitive impairment and similar classifications," Aging and Mental. Health, vol. 7, pp. 238-250, 2003. [42] Y. S. K. Vetsa, M. Styner, S. M. Pizer, J. A. Lieberman and G. E. -. Gerig, "Caudate Shape Discrimination in Schizophrenia using Template-Free Non-Parametric Tests," MICCAI, vol. 2879, pp. 661-669, 2003 42 [43] D. Goldberg-Zimring, A. Achiron, S. K. Warfield, C. R. G. Guttmann and H. Azhari, "Application of spherical harmonics derived space rotation invariant indices to the analysis of multiple sclerosis lesions' geometry by MRI," Magnetic Resonance Imaging, vol. 22, pp. 815-825, 2004. [44] G. Gerig, M. Styner, D. Jones, D. Weinberger and J. Lieberman, "Shape analysis of brain ventricles using SPHARM," pp. 171-178, 2001 [45] R. Cabeza, N. D. Anderson, J. K. Locantore and A. R. McIntosh, "Aging gracefully: compensatory brain activity in high-performing older adults," Neuroimage, vol. 17, pp. 1394-1402, 2002. [46] S. Ogawa, T. M. Lee, A. R. Kay and D. W. Tank, "Brain Magnetic Resonance Imaging with Contrast Dependent on Blood Oxygenation," Proceedings of the National Academy of Sciences, vol. 87, pp. 9868-9872, 1990. [47] D. G. Norris, "Principles of magnetic resonance assessment of brain function," Journal Magnetic Resonance Imaging, vol. 23, pp. 794-807, 2006. [48] D. J. Heeger and D. Ress, "What does fMRI tell us about neuronal activity?" Nature Reviews Neuroscience, vol. 3, pp. 142-151, 2002. [49] E. Amaro Jr and G. J. Barker, "Study design in fMRI: basic principles," Brain Cogn., vol. 60, pp. 220-232, 2006. [50] D. H. Brainard, "The Psychophysics Toolbox," Spatial Visualization, vol. 10, pp. 433-436, 1997. 43 [51] S. C. Strother, "Evaluating fMRI preprocessing pipelines," IEEE Engineering in Medicine Biology and. Magnetic, vol. 25, pp. 27-41, 2006. [52] J. V. Hajnal, R. Myers, A. Oatridge, J. E. Schwieso, I. R. Young and G. M. Bydder, "Artifacts due to stimulus correlated motion in functional imaging of the brain," Magnetic Resonance Medicine, vol. 31, pp. 283-291, 1994. [53] S. C. Strother, J. R. Anderson, X. L. Xu, J. S. Liow, D. C. Bonar and D. A. Rottenberg, "Quantitative comparisons of image registration techniques based on high- resolution MRI of the brain," Journal of Computer Assisted Tomography, vol. 18, pp. 954-962, 1994. [54] T. R. Oakes, T. Johnstone, K. S. Ores Walsh, L. L. Greischar, A. L. Alexander, A. S. Fox and R. J. Davidson, "Comparison of fMRI motion correction software tools," Neuroimage, vol. 28, pp. 529-543, 2005. [55] R. Liao, J. L. Krolik and M. J. McKeown, "An information-theoretic criterion for intrasubject alignment of FMRI time series: motion corrected independent component analysis," IEEE Transaction on Medical Imaging, vol. 24, pp. 29-44, 2005. [56] H. J. Jo, J. Lee, J. Kim, Y. Shin, I. Kim, J. S. Kwon and S. I. Kim, "Spatial accuracy of fMRI activation influenced by volume- and surface-based spatial smoothing techniques," Neuroimage, vol. 34, pp 540-564, 2007. [57] K. J. Worsley, "An overview and some new developments in the statistical analysis of PET and fMRI data," Human Brain Mapping, vol. 5, pp. 254-258, 1997. 44 [58] M. Reimold, M. Slifstein, A. Heinz, W. Mueller-Schauenburg and R. Bares, "Effect of spatial smoothing on t-maps: arguments for going back from t-maps to masked contrast images," Journal of Cerebral Blood Flow Metabolism, vol. 26, pp. 751-759, 2006. [59] A. Geissler, R. Lanzenberger, M. Barth, A. R. Tahamtan, D. Milakara, A. Gartus and R. Beisteiner, "Influence of fMRI smoothing procedures on replicability of fine scale motor localization," Neuroimage, vol. 24, pp. 323-331, 2005. [60] K. J. Friston, A. Mechelli, R. Turner and C. J. Price, "Nonlinear responses in fMRI: the Balloon model, Volterra kernels, and other hemodynamics," Neuroimage, vol. 12, pp. 466-477, 2000. [61] 0. Friman, J. Cedefamn, P. Lundberg, M. Borga and H. Knutsson, "Detection of neural activity in functional MRI using canonical correlation analysis," Magnetic Resonance in Medicine, vol. 45, pp. 323-330, 2001. [62] M. J. Mckeown, S. Makeig, G. G. Brown, T. P. Jung, S. S. Kindermann, A. J. Bell and T. J. Sejnowski, "Analysis of fMRI data by blind separation into independent spatial components," Human Brain Mapping, vol. 6, pp. 160-188, 1998. [63] J. Talairach and P. Tournoux, Co-Planar Stereotaxic Atlas of the Human Brain: 3- Dimensional Proportional System: An Approach to Cerebral Imaging. Thieme, 1988, [64] F. Crivello, T. Schormann, N. Tzourio-Mazoyer, P. E. Roland, K. Zilles and B. M. Mazoyer, "Comparison of spatial normalization procedures and their impact on functional maps," Human Brain Mapping, vol. 16, pp. 228, 2002. 45 [65] C. E. Stark and Y. Okado, "Making memories without trying: medial temporal lobe activity associated with incidental memory formation during recognition," Journal of Neuroscience, vol. 23, pp. 6748-6753, 2003. [66] M. I. Miller, M. F. Beg, C. Ceritoglu and C. Stark, "Increasing the power of functional maps of the medial temporal lobe by using large deformation diffeomorphic metric mapping," Proceedings of the National Academy of Science U. S. A., vol. 102, pp. 9685-9690, 2005. [67] P. A. Yushkevich, J. A. Detre, D. Mechanic-Hamilton, M. A. Fernandez-Seara, K. Z. Tang, A. Hoang, M. Korczykowski, H. Zhang and J. C. Gee, "Hippocampus-specific fMRI group activation analysis using the continuous medial representation," Neuroimage, vol. 35, pp. 1516-1530, 2007. [68] D. V. Vranic, "An improvement of rotation invariant 3D-shape based on functions on concentric spheres," in 2003, pp. III-757-60 vol.2, 2003 [69] M. Kazhdan, T. Funkhouser and S. Rusinkiewicz, "Rotation invariant spherical harmonic representation of 3D shape descriptors," in Symposium on Geometry Processing, pp. 156-65, 2003. [70] V. Caselles, "Geodesic active contours," International Journal of Computer Vision (IJCV), vol. 22, no. 1, pp. 61-79, 1997. 46 2 Unique Invariant SPHARM Features for ROI Based Analysis 2.1 Introduction Magnetic resonance imaging (MRI) is increasingly recognized as a suitable modality for studying the human brain anatomy in vivo. Improvements in both MRI scanner hardware and acquisition sequences are yielding progressively higher resolution 3D images of internal brain structures. One of the main benefits of this increased resolution is in the study of structure of brain anatomy. Researchers use MR images to delineate specific brain areas, like the thalamus, to observe structural changes either over time or across two or more subject groups, e.g. patients with Parkinson's disease (PD) and normal volunteers. This process of delineating specific structures to perform targeted analysis is termed region of interest (ROI) analysis. Numerous methods have been proposed to quantify and analyze the shape of ROIs. The simplest and most common approach uses the overall volume within an ROI, e.g. Anstey et al used the hippocampal volume to investigate changes observed in mild cognitive impairment (MCI) [I J. However, studies using volumetric measures do not reflect actual shape changes that might occur which might be indicative of neurological change. Vetsa A version of this chapter has been published. Uthama A, Abugharbieh R, Traboulsee A, McKeown MJ, "Invariant SPHARM Shape Descriptors for Complex Geometry in MR Region of Interest Analysis" 29th Annual International Conference IEEE Engineering in Medicine and Biology Society 2007, pp. 1322-1325, France 2007 47 et al. used medial representations to characterize the surface of the caudate and identify shape changes in schizophrenia [2]. However, this approach requires accurate alignment of the ROIs across different subjects for valid analysis and hence residual errors in mutual alignment could adversely affect the final results. Alternate methods rely on obtaining descriptive feature vectors which may reflect the structural information in an ROI. These features may be invariant to rotations, translations and scale, thus enabling easy comparison across different subjects. This eliminates the need for mutual alignment, saving time and potential inconsistencies arising from registration errors. In an earlier work, we used scale, rotation and translation invariant features derived using a spherical harmonic (SPHARM) representation of 3D ROIs to analyze shape changes in brain MRI data [3] and demonstrated the advantages of using such an approach along with traditional volumetric analysis. However, the parameterization technique used in our previous work did not allow for shapes with non-convex topology. In this paper, we extend our method and propose new enhanced SPHARM shape descriptors that are unique to any ROI. These features also enable the characterization of complex shapes with complex topology, including those with multiple disjoint parts. We validate our method on synthetic data and demonstrate the effectiveness of the proposed technique even in the presence of inter-subject variability and misalignment between subjects. We also use the proposed technique to investigate shape changes in the thalamus in patients with Parkinson's disease. 2.2 Methods In this section, we briefly summarize our earlier work on SPHARM-based ROI descriptors in MRI. We then proceed to present the proposed enhancements used to 48 enable the study of complex ROI topologies including a novel radial transform used to obtain unique features for any ROI. 2.2.1 Proposed Invariant SPHARM features A 3D ROI with spherical topology can be expressed as P(0,0) , a function defined on the unit sphere with B and 0 as the zenithal and azimuthal angles respectively. The value of the function at a given (6,0) is equal to the distance to the surface point from the centroid at those angles [3]. The SPHARM representation of this function is given by 2.1. fr d0 fiC: 0m0,0)sin(0)ti9. 0^0 (2.1) Yr. (0 ,0) is the complex conjugate of the m th order spherical harmonic of degree 1, where / ranges from 0 to L [4]. The accuracy of this representation depends on the value of L, also called the bandwidth. Invariance of the representation to translation is obtained by shifting the origin of the 3D ROI to its centroid before computing the SPHARM coefficients. In our previous work [3], we employed rotationally invariant features from this representation using (2.2). Scale normalization was achieved by scaling the vector obtained in (2.2) with its first element, which gives a rough measure of total volume [3]. E m *rnC 1 . m=-1 (2.2) Non convex topologies, i.e. surfaces with self occlusions, or topologies with disconnected components (e.g. ventricles) cannot be represented as a function of two 49 variables, 'F(B, O) , since these topologies do not always have a single distance value for all (0 , 0 ). A third variable r, the distance from the origin can be introduced to extend the definition to include such shapes. The augmented function W(r09,0) is then binary valued, taking a non-zero value only in points inside the 3D ROI. The SPHARM representation of such a function is given by: 2" sin(zkr)c: --= fr 2dr PO 111-2^Y,04#1(r,(9,0)sin OB.^(2.3) 0^0^o^r k is an index introduced to account for possible degeneracy [4]. Since direct computation of (2.3) is highly inefficient, we use an alternate approach by representing the data as a set of spherical functions. These functions are obtained by intersecting the 3D binary ROI with R spherical shells of increasing radii r. This intersection is simulated by sampling the 3D ROI on a spherical grid of constant dimensions 2Lx2L at each value of r [5]. The SPHARM coefficients are then obtained at each value of r as: 2ff cm -= fd0 ,fYI , (9, 0)1*, 0, 0)sinHdt9.^(2.4) o^0 This results in a SPHARM representation for each shell. Invariant features for these shells can then be derived using (2.2), as used by Kazhdan et al [6]. However, as noted in [6], this representation is not unique, since independent rotations of the shape along any of the constituting shells will result in the same feature vector representation even though the topologies of the shapes are drastically different. To overcome this limitation, we propose adapting the radial transform in (2.3) across these features as in (2.5). This yields unique features even under independent rotations along the shells. To achieve scale 50 invariance across different ROIs, we propose to keep the number of shells, R, constant at 2R„,„,„ for all ROIs being analyzed. R„,,„, is the radius measure in voxels of the smallest common sphere encompassing each of the ROI in the entire set being analyzed. 2/?. n c;nd Er^2,) ^ ir-2 s n(gIcr) c,71 . r=1 k =^r [1,2,3,...,2R...A^(2.5) Since this radial transform is applied across each value of r, as shown in (2.5), the vector cH must be of the same length for each shell. Hence, we use a common value of L for all shells, whose minimum value is obtained by equating the surface area of the encompassing sphere to the equiangular sampling grid (2.6). ^ 471R.2. 2L x 2L, L =^(2.6) The final unique, scale, translation and rotation invariant features are defined as in (2.7). These values are then reshaped into a single row to provide the final feature vector. k=2Rnwx m=1 ^A (p, q) = E Ec:c:*,p 1...L,q 1...2Rm ^(2.7) k =1 m=-1 Figure 2.1 shows how this method could be used to analyze complex disconnected 3D ROIs like the human brain ventricles. 51 Figure 2.1 (a) The human brain ventricles. The two ventricles would have to be analyzed jointly since the relative orientations might be of importance. One of the spherical shells at r =35 is also shown. (b) The parameterization of the surface at r = 16 (top) and 35 (bottom) 2.2.2 Statistical Analysis of Features The feature vectors thus obtained are of a high dimension, dependent on the value of L chosen. Moreover, the generating probabilities of these features are unknown. We thus employ a non-parametric permutation test for the statistical analysis of these features, since no assumptions of the feature parameters are required. In permutation tests, a distance measure between two groups of shapes (represented by their feature vectors) is tested against the distance obtained by all possible permutations of the vectors across the two groups. We use the Euclidian distance between the feature vectors as the distance measure. Also, as reported by Vetsa et al [2] we use a Monte Carlo approach to obtain a sufficiently high number of permutations, since testing against all 52 possible permutations is not feasible. A final p value is obtained as the ratio of the number of permutations which yielded a Euclidian distance greater than the original ordering of the samples and the total number of permutations performed. If the obtained p value is less than a statistical standard threshold of 0.05, the two groups are declared to have significant differences between them. In our study, this is a strong indication that the shapes in the two groups have a systematic difference in structure. 2.3 Data and Imaging This section describes the generation and acquisition of data used to test our proposed approach, which include synthetic validation data as well as real structural MRI data of patients and controls. 2.3.1 Synthetic Structural Data Generation To verify the performance of the proposed shape description technique, we test on synthetic data with simulated, but nevertheless realistic inter-subject variability. Similar to [7], the basic shape used to generate these data is an ellipsoid as defined in (2.8) where (xc, yc, zc) is the centroid with (a, b, cx) as parameters defining the 3D boundary of the ellipsoid. 0(x , y , z) (x- xc- n (a - n a ,)2 2 - ye - ny,^(z_zc_nz , y 01_ nh,)2^(cx - nc,)2 (2.8) Keeping a and b constant, two values of cx, cc, and cb, were used to generate two groups of shapes, A and B. Twenty samples of these shapes were created in each group. A common centroid of (64, 64, 64) was chosen in a grid of (128, 128, 128) for both groups, 53 keeping with the resolution of real data. Each sample was independently created by perturbing the centroid using (n„„ ny ,, nz,), generated from a Gaussian noise source with zero mean and a standard deviation of 5, N(0,5). The parameters (a, b, cx) were perturbed by (na„ nb„ nc ,) generated from N (0, 0.2). The ellipsoid was then rotated about a random combination of x, y, and/or z axis by an angle generated by a uniform distribution from 0-45°. To create a realistic surface, comparable to the irregular edges obtained from manual ROI segmentation, we add Gaussian noise of N (0, 05) and use a set of morphological operations to obtain the final shape shown in Figure 2.2. (a) Figure 2.2 (a) Surface rendering of a real left thalamus (top) and right thalamus (bottom) from a PD patient, note the rough edges. (b) Group A synthetic samples (see text) generated with parameter set (2, 5, 10). (c) Group B synthetic sample generated with parameter set (2, 5, 8). Note the realistic looking edges in the synthetic data. 2.3.2 Real Structural MR Data Acquisition The data used in this study consisted of MRI scans of 21 control subjects and 19 PD patients. Two scans were performed with a gap of two hours, one before and one after the administration of the drug L-dopa. The MR data were acquired on a 3.0 Tesla Siemens 54 scanner (Siemens, Erlangen, Germany) with a birdcage type standard quadrature head coil and an advanced nuclear magnetic resonance echoplanar system. The participant's head was positioned along the canthomeatal line. Foam padding was used to limit head motion within the coil. High-resolution T1 weighted anatomical images (3D SPGR, TR=14ms, TE=7.7 ms, flip angle=25°, voxel dimensions 1.0x1.0x1.0mm, 176x256 voxels, 160 slices) were acquired. ROIs for the left and the right thalamus were manually drawn on both the scan data sets using the IRIS 2000 software (NeuroLib, University of North Carolina) by a trained research assistant using a computer mouse. The research assistant was blinded to the disease category. 2.4 Results and Discussion This section presents quantitative results of employing the proposed SPHARM shape features on the synthetic and real MR data outlined in Section 2.3. 2.4.1 Validation on Synthetic Data Synthetic data generated as described earlier consisted of forty one data sets each set having two groups, A and B, each group with twenty sample shapes. For the entire set, group A was re-generated independently each time with a baseline parameter set (2, 5, 10). Group B was generated using a different baseline parameter set each time, starting with (2, 5, 8) and ending with (2, 5, 12), with the value of cb being incremented by 0.1 for each set. The proposed SPHARM invariants were then derived for all samples in each set. Rmax was found to be22; hence 44 shells were used to parameterize the ROIs. A value of 40 was chosen for L as per (2.6), resulting in a feature vector length of 1760 elements. 55 Non Significant * Significantm 8 coco 4-a-3 463 OA ca :5 0.2 — 1.4 1.2 — 0.8 — 0.6 — 0 8 8.5^9^9.5^10 11.5^12 The Euclidian distance between the SPHARM features of the two groups, A and B, in each set, was recorded. Permutation tests were then performed for each set to determine the level of significance in the difference between the features of the two groups. Figure 2.3 summarizes the results of this analysis. Cb Figure 2.3 Plot of the Euclidian distance between the two groups A and B for different values of cb. Distances resulting in a significant difference (alpha = 0.05) between the groups as obtained by the permutation test are indicated with a star. Non-Significant distances are indicated with a dot. The graph in Figure 2.3 shows the expected trend in the Euclidian distance measure between the two groups. For values of cb smaller than 9.5 and larger than 10.6, the features are able to detect a systematic difference in shape between the two groups. However, as the value of cb approaches 10 (the value of ca), there is less difference in the sample shapes between the two groups. For values of c b in the range (9.5, 10.6), the combined effect of the various perturbations in the parameters involved, coupled with the 56 inter-subject variability introduced, leaves almost no difference in constituting shapes of the two groups. This is shown in Figure 2.3 with dots, indicating that the permutation test was unable to find any systematic difference in the features vectors. This experiment demonstrates the sensitivity of the proposed technique in discriminating subtle shape changes between two groups. 2.4.2 Analysis of Real MR Structural Data The proposed technique was next used to analyze the structure of the left and the right thalamus in patients with Parkinson's disease. SPHARM features were computed for the data of 21 control subjects and 19 PD patients. R„,„, was found to be 20 voxels; hence 40 shells were used to parameterize the ROIs. The bandwidth L was chosen to be 36 as per (2.6), resulting in a feature vector length of 1440. Since two scans were performed for each subject, two feature vectors were obtained for each ROI, for each subject. Also, since the duration between the scans was less than two hours, no systematic structural changes were expected between these two observations. The only changes expected were due to repositioning, manual tracing inconsistencies and other imaging artifacts. To confirm this assumption, each subject group (Patients and Controls) was further sub- divided into two categories: pre-drug and post-drug. Permutation test results of the SPHARM feature vectors between these two sub-groups are shown in Table 2.1. Results from the volumetric analysis are also shown. 57 Table 2.1 Results of pre-drug vs. post-drug analysis Subject group ROI SPHARM (p Value) Volumetric (p Value) Controls Left Thalamus 0.3863 0.3270 Controls Right Thalamus 0.2828 0.2720 PD Left Thalamus 0.4738 0.2430 PD Right Thalamus 0.2341 0.3208 The results in Table 2.1 indicate that, as expected, there are no statistically significant changes in the shape of the thalamus between the two scans. Having established this, the two feature vectors from the two sessions were then averaged to generate a single feature vector for each ROI for each subject. Direct averaging is valid since these features represent the same shape. Changes due to rotation and translations are annulled by the invariance property of the features. This averaging reduces the variability caused by manual tracing errors and other noise sources, thereby increasing the sensitivity of the features to actual shape changes. A permutation test on the averaged features was then performed between the two subject groups. Results are indicated in Table 2.2; results from the volumetric analysis are also shown. 58 Table 2.2 Results of PD patients vs. controls analysis ROI SPHARM (p value) Volumetric (p value) Left Thalamus 0.0012* 0.0102* Our earlier work [3] on a smaller number of subjects observed differences in the volume of both the thalamus in patients with PD, when compared against healthy volunteers. However, the previous technique was unable to detect significant shape changes between the groups. In this present study, in addition to the changes observed with volumetric measures, we detect a significant difference in the shape of both thalami. The shape analysis is independent of the volume changes observed since the features derived are scale invariant. The new parameterization approach coupled with the increased number of subjects probably contributed to this increased sensitivity. The results indicate that PD patients studied have a marked difference in the volume of their thalamus when compared to control subjects, indicating possible atrophy of this structure. Further work is required to determine if the consistent atrophy in PD subjects corresponds to specific thalamic nuclei (substructures within the thalamus). 2.5 Conclusions We enhanced our earlier work generating invariant SPHARM features for 3D ROI analysis in MR data by extending our characterization capability to ROIs with complex topologies, including those with possible disconnections. We also proposed the use of a novel radial transform to obtain unique feature. The proposed features were able to detect subtle shape changes in synthetically created data with realistic inter-subject variations 59 and detect shape changes in the thalami of PD patients in addition to the volumetric changes observed. The use of the proposed method to explore shape changes of different brain structures as a biomarker of disease progression warrants further study. 60 2.6 References [1] K. J. Anstey and J. J. Mailer, "The role of volumetric MRI in understanding mild cognitive impairment and similar classifications," Aging Mental Health, vol. 7, pp. 238-250, 2003. [2] Y. S. K. Vetsa, M. Styner, S. M. Pizer, J. A. Lieberman and G. E. -. Gerig, "Caudate Shape Discrimination in Schizophrenia using Template-Free Non-Parametric Tests" in MICCAI 2003, vol. 2879, pp. 661-669, 2003. [3] S. Tootoonian, R. Abugharbieh, X. Huang and M. J. McKeown, "Shape vs. volume: Invariant shape descriptors for 3D region of interest characterization in MRI," in ISBI, pp. 754-757, 2006. [4] G. Burel and H. Henocq, "Three-dimensional invariants and their application to object recognition," Signal Process, vol. 45, pp. 1-22, 1995. [5] D. M. Healey Jr, D. N. Rockmore and S. B. Moore, "An FFT for the 2-sphere and applications," in ICASSP 1996, pp. 1323-1326 vol. 3, 1996. [6] M. Kazhdan, T. Funkhouser and S. Rusinkiewicz, "Rotation invariant spherical harmonic representation of 3D shape descriptors," in Symposium on Geometry Processing, 23-25 June 2003, pp. 156-65, 2003 [7] D. Goldberg-Zimring, A. Achiron, S. K. Warfield, C. R. G. Guttmann and H. Azhari, "Application of spherical harmonics derived space rotation invariant indices to the 61 analysis of multiple sclerosis lesions' geometry by MRI," Magnetic Resonance Imaging, vol. 22, pp. 815-825, 2004. 62 3 Detection and Analysis of Anatomical Shape Changes in MRI (PD thalamus) 3.1 Background The thalamic changes seen in Parkinson Disease (PD) may represent selective non- opaminergic degeneration [1], as there is selective neuronal loss in the centre median- parafascicular (CM-Pf) complex in Parkinson's disease [2], yet relative preservation of neurons in the limbic (mediodorsal and anterior principal) thalamic nuclei. Henderson et al. examined the thalamic intranuclear nuclei in 10 normal controls and 9 patients with PD [3]. As expected, they found a-synuclein-positive Lewy bodies in these nuclei in the thalamus, but they also found a significant reduction (40-55%) in the total neuronal number in the caudal intralaminar (CM-Pf) nuclei, regions that receive glutaminergic innervation [3]. This contrasted with the 70% loss of pigmented nigral neurons. A factor analysis has demonstrated that the size of neurons in the motor cortex is negatively correlated with the size and number of neurons in its thalamic relay, the VLp. There is also a positive correlation between the number of ventral anterior (VA) neurons and the pre- supplementary motor area (SMA) [4]. A version of this chapter has been submitted for publication. Martin J. McKeown, Ashish Uthama, Rafeef Abugharbieh, Samantha Palmer, Mechelle Lewis and Xuemei Huang, "Shape (But Not Volume) Changes in the Thalami in Parkinson Disease", Submitted to BioMed Central Neuroscience, 2007 A version of this chapter has been accepted for a poster presentation. Martin J McKeown, Ashish Uthama, Rafeef Abugharbieh, Samantha Palmer, Mechelle Lewis and Xuemei Huang. "MRI in Parkinson's Disease identifies shape, but not volume changes in the thalamus", Abstract, Selected for poster presentation at XVII WFN World Congress on Parkinson's Disease and Related Disorders, Amsterdam, The Netherlands, 2007 63 Bacci et al. suggested that CM-Pf degeneration may partially counteract the consequences of dopamine neuronal loss, as thalamic and dopamine inputs have antagonistic influence on neurotransmitter-related gene expression [1]. Moreover, the CM-Pf degeneration may be a direct consequence of nigrostriatal denervation, as depleting the striatum of dopamine results in the remaining Pf neurons being particularly hyperactive [5]. The normal role of the CM-Pf complex is incompletely understood, although it is clearly related to basal ganglia function [6]. The Pf nuclei receive input from the spinal cord and project to the striatum [7]. These projections may carry specific temporally-patterned inputs to striatal targets [8]. While the CM-Pf complex has traditionally been considered part of the reticulo-thalamo-cortical activating system, a recent proposal suggests that the CM-Pf complex participates in sensory driven attention processes, particularly unexpected events [9]. The advent of modern imaging techniques has allowed the non-invasive in vivo assessment of brain structures, such as the thalamus, in disease states. Thalamic morphological changes have been detected chronically after cortical injury, such as middle cerebral artery (MCA) infarction after 1 yr [10, 11] and tumor resection after —2yrs [12]. The study by Hulshoff Pol [12] detected a 5% decrease in ipsilateral thalamus and, interestingly, a 4.5% increase in contralateral thalamic volume after unilateral tumor resection, presumably on the basis of a compensatory mechanism. In PD and related disorders, some studies of structural and functional imaging have detected thalamic changes. Thalamic grey matter changes have been detected 64 contralateral to unilateral Parkinsonian resting tremor [13]. In PD with dementia, the thalamus, in addition to the hippocampus and anterior cingulate, represent the regions most affected [14]. Functionally, there is a connection between a major component in F- DOPA uptake in the striatum and a component from fluoro-deoxy--glucose (FDG), which had positive loadings in the thalamus and the cerebellum [15]. Most morphological studies based on imaging involve a number of steps manipulating the brain images. A typical approach would be to warp ("spatially normalize") the brain images to a common space [13]. Further smoothing of the data (e.g. using an isotropic 12mm Gaussian kernel) to minimize the effects of misregistration between different normalized brains may affect the ability to make inferences about small, subcortical structures like the thalamus. In fact this "Voxel Based Morphometry" approach has thus been a controversial approach (e.g., see discussions in [16] and [17]). Recent approaches try to reduce errors due to misregistration by aligning the subjects at the region of interest (ROI) level, as opposed to the whole brain level [18] [19] [20]. However, these approaches are designed to deal with a different problem, namely that of summarizing fMRI activation from several subjects. To quantify differences in morphology, it would be necessary to examine the different transformations required to warp each subject's ROIs to the exemplar ROI shape -- a non-trivial task. An alternative approach to warping brains to the same space is to segment brain structures individually on unmanipulated (i.e. unregistered and unwarped) brains [21]. Because no registration of the brain images is done, this requires summarizing the individual brain structures in a way that they are invariant to positioning of the head in 65 the scanner. For example, simply estimating the volume of an ROI such as the thalamus has this property, as it is invariant to the individual coordinate frame used. A number of such invariants (e.g. spatial variance) have been proposed to summarize the shape of brain structures [22], or even characterize the distribution (texture) of activation maps in fMRI [23]. We have recently proposed a method based on spherical harmonics (SPHARM) which provided a unique representation of brain structures, including regions with possible topological disconnections, such as the lateral ventricles [24]. In brief, the method involves first finding a spherical shell which encompasses the entire ROI. Subsequent smaller concentric shells are then derived and the intersection between the progressively smaller spherical shells and the brain structure is computed (Figure 3.1). The results from the intersections are then combined into a unique feature vector containing approximately 100's or 1000's of elements. This feature vector provides a unique representation of the shape which is independent of the spatial orientation of the structure (see section 3.5). We examined the thalami from 18 PD subjects and 18 age-matched controls. Using the above technique, we found significant differences between the two groups in the shape of the thalami, but not in the volume. This suggests that significant thalamic changes can be assessed noninvasively in PD, suitable for longitudinal studies. 66 Concentric Spheres r=31 Thalamus ED Figure 3.1 The SPHARM-based method for shape determination. The shape to be specified (the thalamus) and two concentric spherical shells are shown. On the right is the intersection between the thalamus and shells as a function of rotation (0) and azimuth ((p). The rotation angle spans from 0 to 27t radians, and the azimuth angle is from 0 to rc radians. 3.2 Results There were no significant differences in volume between sides in either controls or PD subjects, nor between controls and PD subjects in either the left or right thalamus (Table 3.1). In contrast, the SPHARM-based method found significant shape differences between the left and right thalamus in PD but not in controls. Significant shape differences between PD and controls were detected in both the left and right thalami. 67 Table 3.1 Results of volumetric and shape analysis. Numbers indicate the p-values obtained from the permutation test Group Volume SPHARM Control, Left vs. Right 0.5630 0.1470 PD, Left vs. Right 0.5780 0.0060 Left thalamus, PD vs. Control 0.4150 0.0270 Right thalamus, PD vs. Control 0.1730 0.0290 Multivariate regression did not find any statistically significant (p < 0.05) effects of handedness, side of symptoms, or dominant side of tremor on the distance measure. In order to better visualize and intuitively assess the shape differences in thalami between PD subjects and controls, we took thalami that had "typical" feature vectors (i.e. feature vectors closest to the mean of each group) and assumed that they represented exemplar shapes. We then spatially aligned these exemplar shapes (Figure 3.2). There appeared to be greater differences in the left thalami between controls and PD subjects. The largest differences appeared to be along the dorsal surface. Note that the registration of the thalami in this instance was solely for visualization purposes and was not incorporated into the analysis. 68  Normal PD Posterior^ Posterior Rigid registration of the Thalami PDNormal Final rendering Figure 3.2 Two sets typical thalami registered and shown on brain from a PD subject. Note that the registration here is solely for visualization purposes, and is not required for the calculation of shape differences. Also, although the thalami here were first smoothed with a 12mm FWHM Gaussian kernel for visualization purposes, no smoothing was performed for the shape analysis and group comparison. 69 3.3 Discussion It is well known that changes in the thalamus can be seen chronically after cortical injury [25, 26]. This is not only due to direct effects of axonal damage, as axonal-sparing cortical lesions also result in thalamic degeneration [27]. Such thalamic degeneration probably involves both anterograde and retrograde processes [28] and may be mitigated by growth factors [29] [30]. Brain development has a critical role on the extent of thalamic changes after a cortical lesion. Animal models have determined that perinatal lesions are far less likely to induce thalamic changes, compared to when the cortical lesions are made prenatally [31] or in adulthood [32]. . In contrast to the secondary effects of thalamic changes from cortical lesions, the thalamic changes in PD are related to selective non-dopaminergic neurodegeneration [1]. Consistent with prior results, we found no significant differences in the volume of the thalami between PD subjects and controls [3]. However, for the first time, we have demonstrated that the shape of the thalami undergoes systematic changes in PD. The reason that shape may change but not the volume may be due to the fact that particular nuclei (e.g. CM-Pf) are involved, and thus, at the typical resolution of MRI, this does not result in significant changes in the overall volume. Alternately, since thalamic volume may actually increase as a compensatory mechanism [12], other regions of the thalamus may hypertrophy. The ability to non-invasively quantify subtle morphological shape changes appears to be a powerful technique. We utilized standard structural MR imaging techniques without 70 any special sequences or any special scanner resolution requirements. We obtained robust results from pooling data from two different centers using two different types of scanners. We used manual segmentation of the thalami in this paper. Automatic segmentation of subcortical structures is an area of ongoing research [33], and often requires the tuning of many parameters, especially when the images may be pooled from scanners from different centers. Since the person at each centre doing the segmentation was blinded to disease status, it would be unexpected that a systematic bias was introduced into our results. Even then, any misspecification of the same ROI across subjects would tend to increase inter-subject variability and presumably reduce discriminability across groups making the task harder for the shape analysis approach. We used SPHARM-based invariant descriptors to quantify the shape of the thalami. The main advantages of such a method is that it does not require that all brain images be warped to a common space, nor does it require that the brain images be aligned in any way. A drawback of these invariant features approach is that it is difficult to invert the feature vectors, i.e. once given all the values of the different invariants, it is impossible to reconstruct the original image which gave that feature vector -- analogous to the fact that given the volume of an object, it is impossible to uniquely reconstruct the original shape. We therefore cannot create a "typical" brain structure by averaging the feature vectors and create the image that would give this feature vector (e.g. to create an "average left thalamus"). However, it is possible to cluster structures in the feature space and find the brain structure whose feature vector is in the middle of the cluster so as to use it as an exemplar shape, which we have done (Figure 3.2). 71 It is difficult to determine if the shape differences we detected are attributable to any specific nuclei. However, based on prior pathological studies, it would be likely that the differences we detected were related to degeneration in the CM-Pf complex. Given that progressive supranuclear palsy (PSP) has even greater involvement of VLp than PD [4], it remains to be seen if thalamic shape is a discriminable feature between these two conditions. We did not detect any association between overall shape change and handedness, and dominant side presentations or presence/absence of tremor. This may be due to the relatively small sample size employed in this study. However, because the feature vectors consist of many different components, we don't discount that there may be a subset of components that are sensitive to these disease parameters. 3.4 Conclusions Our results suggest that systematic changes in thalamic shape can be non-invasively assessed in PD in vivo and that shape changes, in addition to volume changes, may represent a new avenue to assess the progress of neurodegenerative processes. Although we cannot state which parts of the thalamus are directly affected, previous pathological studies would suggest that the shape changes detected in this study represent degeneration in the centre median-parafascicular (CM-Pf) complex, an area known to represent selective non-dopaminergic degeneration in PD. 72 3.5 Methods The study was approved by the appropriate Institutional Review Boards and Ethics Boards of the University of British Columbia (UBC) and the University of North Carolina (UNC). All structural data were obtained as part of fMRI studies whose results are reported elsewhere (e.g., [34]). 3.5.1 MR Imaging at the University of British Columbia All subjects gave written informed consent prior to participating. Nine volunteers with clinically diagnosed PD participated in the study (5 men, 4 women, mean age 68.1 ± 6.8 years, 7 right-handed, 2 left-handed). All subjects had mild to moderate PD (Hoehn and Yahr stage 2-3) [35] with mean symptom duration of 3.6 ± 2.6 years. We recruited ten healthy, age-matched control subjects without active neurological disorders (3 men, 7 women, mean age 55 ± 12.4 years, 9 right-handed, 1 left-handed). Exclusion criteria included atypical Parkinsonism, presence of other neurological or psychiatric conditions and use of antidepressants, sleeping tablets, or dopamine blocking agents. MRI was conducted on a Philips Achieva 3.0 T scanner (Philips, Best, The Netherlands) equipped with a head-coil. A high resolution, three dimensional (3D) Tl-weighted image consisting of 170 axial slices was acquired of the whole brain to facilitate anatomical localization of activation for each subject. The thalami were one of eighteen specific regions of interest (ROIs) that were manually drawn on each unwarped, aligned structural scan using the Amira software (Mercury Computer Systems, San Diego, USA). Although the thalami were manually segmented 73 on the axial slices, they were carefully examined in the coronal and saggital planes to ensure accuracy. The trained technician performing the segmentation was blinded to the disease state. 3.5.2 MR Imaging at the University of North Carolina All subjects gave written informed consent prior to participating. Nine volunteers with clinically diagnosed mild to moderate PD (Hoehn and Yahr stage 2-3 -- mean symptom duration of 2.1 ± 2.0 years) participated in the study (5 men, 4 women, mean age 58 ± 12yrs, all right-handed). We also recruited eight healthy, age-matched control subjects without active neurological disorders (5 men, 3 women, mean age 49 ± 14yrs, 8 right- handed). Images were acquired on a 3.0 Tesla Siemens scanner (Siemens, Erlangen, Germany) with a birdcage-type standard quadrature head coil and an advanced nuclear magnetic resonance echoplanar system. The head was positioned along the canthomeatal line. Foam padding was used to limit head motion. High-resolution T1 weighted anatomical images were acquired (3D SPGR, TR=14 ms, TE=7700 ms, flip angle=25°, voxel dimensions 1.0 x 1.0 x 1.0 mm, 176x256 voxels, 160 slices). ROIs (including thalami) were drawn manually by the same trained research associate with assistance from multiple on-line and published atlases (e.g. [36]). 3.5.3 Thalami Shape Analysis As described in the technical appendix, the analysis of each shape results in a unique feature vector, of length n = 1440. The left and right thalami were analyzed separately. 74 For comparison, we examined for any differences in volume. The volume of each thalamus was estimated as the number of voxels that each ROI contained multiplied by the volume of a single voxel. To assess the significance of group differences between feature vectors, we used a permutation test to generate a null distribution of Euclidean distances between feature vectors. The permutation test does not require a priori assumptions about the data distribution, and is thus preferred over T-test and F-test [37]. We assessed the differences in left vs. right thalami in controls, left vs. right in PD subjects, PD vs. controls for the left thalamus, and PD vs. controls for the right thalamus. 3.5.4 Shape Analysis -- Technical Aspects Let '11(B4 O) be a function defined on the unit sphere with 0 and 0 as the zenithal and azimuthal angles, respectively. The SPHARM representation for this function is given by (3.1) where (B, 0) is the complex conjugate of the m th order spherical harmonic of degree 1. 1 ranges from 0 to L [16]. Increasing the value of L, also called the bandwidth, improves the representation accuracy at the cost of higher computation time. This definition can also be extended to real valued 3D distributions T(r,O,0) (3.2), where r is the distance from the origin to a given voxel. k is an index introduced to account for possible degeneracy due to the additional dimension [38]. 2n^n C im = 1d0 f Y„:, OWO, 0)sin 0^0 (3.1) 75 c: =1/- 2 cir 'ffid011i sin (ffkr) Y,,,(0,0)`F(r , 0 , Osin (19)c 1 0 0^0^0^r As explained later in this section, rotationally invariant features can be derived from this spherical harmonic representation. In our application, we also need the features to be invariant to any translation of the entire ROI in 3D space. To achieve this, we move the origin of the function T(r,0,0), to the centroid of Ts (r,0,0) , where Tjr,6),Ø) is given by (3.3). , {1 if Tfr,9,0) ^ 0 `I', fr, 19,0)' 0 if T(r,19,0).= 0 (3.3) Since direct computation of (3.2) is highly inefficient [39], we use an alternate approach by representing the data as a set of spherical functions obtained by intersecting the 3D data with spherical shells. Alternatively, for each value of r, T(r,t9,0) can be visualized as a spherical shell comprising the function values at a distance r from the origin. r can then be incremented in steps of t to encompass the entire ROI. If the initial representation of the function is in the form of a cubic grid (regular isotropic voxels in our case), volumetric interpolation is required to resample the ROI in the spherical coordinate space. When analyzing multiple subjects' ROIs simultaneously, we define the maximum radius, R,,„ as the minimum radial distance in voxel count that encompasses all non-zero values of all subjects' ROIs being analyzed. To represent the values from the cubic grid of all (3.2) 76 ROIs with sufficient accuracy, 2R„,„„ shells are used. To achieve scale invariance, the shells must be distributed evenly throughout the spatial extent of each ROI. Since the ROI size across subjects is not uniform, shell spacing t must be adjusted for each subject separately. This procedure ensures that each shell captures similar features from the 3D ROI irrespective of its scale. Surface sampling along each of these shells is performed on an equiangular spherical grid of dimensions 2Lx2L [39]. The common bandwidth L for all shells of all functions is chosen to satisfy the sampling criterion for the largest shell in this set of ROI, namely the one with radius R„,,,,,.. The surface area for this shell represents the maximum surface shell area that needs to be sampled by the equiangular grid; hence, any value of L satisfying the required equiangular sampling (2L x2L) at this shell will be sufficient to represent data from smaller radii shells. The minimum value for L is obtained by equating the surface area of this largest shell to the equiangular sampling grid (3.4). Higher values of L are not used, since it increase computation time with no added benefit. Also, this will result in longer feature vectors, complicating the analysis. Furthermore, when the represented object is a discrete array, higher values of I, resulting from a larger L, may correspond to sampling noise [38]. Recognizing that in applications pertaining to discrimination, high accuracy in the SPHARM representation is not a necessity, we chose to use the minimum value for L as that obtained by (3.4). 4•.1? 2 . =2Lx2L, L = R max,ITT (3.4) To obtain the SPHARM representation for all shells, a discrete SPHARM transform is performed at each value of r to obtain c r'n, (3.5). Features derived from this representation, 77 however, do not provide a unique function representation [40, 41]. For instance, rotating the inner and outer shells by different amount will result in different spatial distribution of function values. However, in this approach, the derived features are insensitive to these rotational transforms, thus resulting in the same feature values for dissimilar spatial distributions. 2r r C 'r7 = f d0 SY,,: i (19,0>iqr,t9,0)s in McIO r [1,2,3,...,2R.]^(3.5) Burel and Henocq's original equation (3.2) does not have this problem, since a part of the basis function is a function of r. However, since (3.2) is computationally intractable [17], we proposed an efficient approach that uses a radial transform (3.6), derived from (3.2), to obtain a unique function representation. The transform (3.6) retains the relative orientation information of the shells, thus the features derived will be sensitive to independent rotations of the different shells, thereby ensuring that unique feature representation is obtained. 2Rm,m ^2C = Li r r=1 sin(gkr) ,„c k =^ (3.6) The range of k could be changed to obtain different lengths of the final feature vector. However, to avoid unnecessarily increasing the feature vector length or losing any 78 important information caused by reducing the range, we choose to keep the range of k the same as that of r, i.e. From the obtained representation (3.6), we then compute similarity transform invariant features using (3.7) for each value of / and k [38] with p and q are used to index these features. Note we reshape I into a single row vector of dimensions D = L x 2R.a. for later analysis. Ic=2Rm.„ m=1 I(p,q)= E E (cr,)*, p =1...L, q 1...2Rma„ k=1 3.6 Authors' contributions MJM conceptualized the study, supervised the collection of the data from UBC, and supervised the application of the SPHARM technique to medical data. AU developed the SPHARM-based method and performed the calculations. RA supervised the development of the SPHARM-based method and assisted in the application to medical data. SP collected the data at UBC and performed the manual segmentations of the thalami. ML collected the data at UNC and manually segmented the data from UNC. XH supervised the collection of UNC data and assisted in biological interpretation of the results. 3.7 Acknowledgements This work was supported by a grant from NSERC/CIHR CHRP-(323602 - 06) (MJM). (3.7) 79 3.8 References [I] Bacci J-J, Kachidian P, Kerkerian-Le Goff L, Salin P: Intralaminar thalamic nuclei lesions: widespread impact on dopamine denervation-mediated cellular defects in the rat basal ganglia. Journal of Neuropathology and Experimental Neurology, 63(1):20-31, 2004 [2] Henderson JM, Carpenter K, Cartwright H, Halliday GM: Degeneration of the centre median-parafascicular complex in Parkinson's disease. Annals of Neurology 2000, 47(3):345-352, 2000. [3] Henderson JM, Carpenter K, Cartwright H, Halliday GM: Loss of thalamic intralaminar nuclei in progressive supranuclear palsy and Parkinson's disease: clinical and therapeutic implications. Brain 2000, 123(Pt 7):1410-1421, 2000. [4] Halliday GM, Macdonald V, Henderson JM: A comparison of degeneration in motor thalamus and cortex between progressive supranuclear palsy and Parkinson's disease. Brain 2005, 128(Pt 10):2272-2280, 2005 [5] Aymerich MS, Barroso-Chinea P, Perez-Manso M, Munoz-Patino AM, Moreno-Igoa M, Gonzalez-Hernandez T, Lanciego JL: Consequences of unilateral nigrostriatal denervation on the thalamostriatal pathway in rats. European Journal of Neuroscience 2006, 23(8):2099-2108, 2006. [6] Kerkerian-Le Goff L, Bacci J-J, Salin P, Aymerich MS, Barroso-Chinea P, Obeso JA, Lanciego JL: Intralaminar Thalamic Nuclei are Main Regulators of Basal Ganglia: 80 Possible involvement in the pathophysiology of Parkinson's disease. In The Basal Ganglia VIII. Edited by Bolam JP, Ingham CA, Magill PJ: Springer US;331-339, 2005. [7] Nakamura Y, Otake K, Tokuno H: The parafascicular nucleus relays spinal inputs to the striatum: an electron microscope study in the rat. Neuroscience Research 2006, 56(1):73-79,2006 [8] Lacey CJ, Bolam JP, Magill PJ: Novel and distinct operational principles of intralaminar thalamic neurons and their striatal projections. Journal of Neuroscience 2007, 27(16):4374-4384, 2007. [9] Kimura M, Minamimoto T, Matsumoto N, Hori Y: Monitoring and switching of cortico-basal ganglia loop functions by the thalamo-striatal system. Neuroscience Research 2004, 48(4):355-360, 2004. [10] Ogawa T, Yoshida Y, Okudera T, Noguchi K, Kado H, Uemura K: Secondary thalamic degeneration after cerebral infarction in the middle cerebral artery distribution: evaluation with MR imaging. Radiology 1997, 204(1):255-262, 1997. [11] Tamura A, Tahira Y, Nagashima H, Kirino T, Gotoh 0, Hojo S, Sano K: Thalamic atrophy following cerebral infarction in the territory of the middle cerebral artery. Stroke 1991, 22(5):615-618, 1991. [12] Hulshoff Poi 14E, van der Flier WM, Schnack HG, Tulleken CA, Ramos LM, van Ree JM, Kahn RS: Frontal lobe damage and thalamic volume changes. Neuroreport 2000, 11(13):3039-3041, 2000. 81 [13] Kassubek J, Juengling FD, Hellwig B, Spreer J, Lucking CH: Thalamic gray matter changes in unilateral Parkinsonian resting tremor: a voxel-based morphometric analysis of 3-dimensional magnetic resonance imaging. Neuroscience Letters 2002, 323(1):29-32, 2002. [14] Summerfield C, Junque C, Tolosa E, Salgado-Pineda P, Gomez-Anson B, Marti MJ, Pastor P, Ramirez-Ruiz B, Mercader J: Structural brain changes in Parkinson disease with dementia: a voxel-based morphometry study. Archives of Neurology 2005, 62(2):281- 285, 2005. [15] Kaasinen V, Maguire RP, Hundemer HP, Leenders KL: Corticostriatal covariance patterns of 6-[18F] luoro-L-dopa and [18F] luorodeoxyglucose PET in Parkinson's disease. Journal of Neurology 2006, 253(3):340-348, 2006. [16] Bookstein FL: "Voxel-based morphometry" should not be used with imperfectly registered images. Neuroimage 2001, 14(6):1454-1462, 2001. [17] Ashburner J, Friston KJ: Why voxel-based morphometry should be used.[comment] Neuroimage 2001, 14(6):1238-1243., 2001. [18] Miller MI, Beg MF, Ceritoglu C, Stark C: Increasing the power of functional maps of the medial temporal lobe by using large deformation diffeomorphic metric mapping. Proceedings of the National Academy of Sciences of the United States of America 2005, 102(27):9685-9690, 2005. 82 [19] Stark CEL, Okado Y: Making memories without trying: medial temporal lobe activity associated with incidental memory formation during recognition. Journal of Neuroscience 2003, 23(17):6748-6753, 2003. [20] Yushkevich PA, Detre JA, Mechanic-Hamilton D, Fernandez-Seara MA, Tang KZ, Hoang A, Korczykowski M, Zhang H, Gee JC: Hippocampus-specific fMRI group activation analysis using the continuous medial representation. Neuroimage 2007, 35(4):1516-1530, 2007. [21] McKeown M, Li J, Huang X, Wang ZJ: Local Linear Discriminant Analysis (LLDA) For Inference Of Multisubject fMRI Data. . In IEEE Acoustics, Speech and Signal Processing: 2007; Honolulu, Hawaii; 2007:305 — 308, 2007. [22] Mangin JF, Poupon F, Duchesnay E, Riviere D, Cachia A, Collins DL, Evans AC, Regis J: Brain morphometry using 3D moment invariants. Medical Image Analysis 2004, 8(3):187-196, 2004 [23] Ng B, Abu-Gharbieh R, Huang X, McKeown MJ: Characterizing fMRI Activations within Regions of Interest (ROIs) Using 3D Moment Invariants, . In Proceedings of Conference on Computer Vision and Pattern Recognition: 2006; New York; 2006:63,2006 [24] Uthama A, Abugharbieh R, Traboulsee A, McKeown MJ: Invariant SPHARM Shape Descriptors for Complex Geometry in MR Region of Interest Analysis. In 29th Annual International Conference IEEE EMBS: 2007; Lyon, France; 2007. 83 [25] Matthews MA: Death of the central neuron: an electron microscopic study of thalamic retrograde degeneration following cortical ablation. Journal of Neurocytology 1973, 2(3):265-288, 1973 [26] Vicedomini JP, Corwin JV, Nonneman AJ: Role of residual anterior neocortex in recovery from neonatal prefrontal lesions in the rat. Physiology & Behavior 1982, 28(5):797-806, 1982. [27] Ross DT, Ebner FF: Thalamic retrograde degeneration following cortical injury: an excitotoxic process? Neuroscience 1990, 35(3):525-550, 1990. [28] Sorensen JC, Dalmau I, Zimmer J, Finsen B: Microglial reactions to retrograde degeneration of tracer-identified thalamic neurons after frontal sensorimotor cortex lesions in adult rats. Experimental Brain Research 1996, 112(2):203-212, 1996. [29] Kumon Y, Sakaki S, Watanabe H, Nakano K, Ohta S, Matsuda S, Yoshimura H, Sakanaka M: Ciliary neurotrophic factor attenuates spatial cognition impairment, cortical infarction and thalamic degeneration in spontaneously hypertensive rats with focal cerebral ischemia. Neuroscience Letters 1996, 206(2-3):141-144, 1996. [30] Yamada K, Kinoshita A, Kohmura E, Sakaguchi T, Taguchi J, Kataoka K, Hayakawa T: Basic fibroblast growth factor prevents thalamic degeneration after cortical infarction. Journal of Cerebellar Blood Flow and Metabolism 1991, 11(3):472-478, 1991. [31] Jackson GF, 3rd, Villablanca JR, Loopuijt LD: Few neocortial and thalamic morphological changes after a neonatal frontal cortical ablation contrast with the effects 84 of a similar lesion in fetal cats. Brain Research 1995, Developmental Brain Research. 90(1-2):62-72. [32] Villablanca JR, Burgess JW, Benedetti F: There is less thalamic degeneration in neonatal-lesioned than in adult-lesioned cats after cerebral hemispherectomy. Brain Research 1986, 368(2):211-225, 1986. [33] Amini L, Soltanian-Zadeh H, Lucas C, Gity M: Automatic segmentation of thalamus from brain MRI integrating fuzzy clustering and dynamic contours. IEEE Transactions on Biomedical Engineering 2004, 51(5):800-811, 2004. [34] Lewis M, Slagle C, Smith A, Truong Y, Bai P, McKeown M, Mailman R, Belger A, Huang X: Task specific influences of Parkinson's disease on the striato-thalamo-cortical and cerebello-thalamo-cortical motor circuitries. Neuroscience Letters 2007, 147(1):224- 235, 2007. [35] Hoehn MM, Yahr MD: Parkinsonism: onset, progression and mortality. Neurology 1967, 17(5):427-442, 1967. [36] Damasio H: Human Brain Anatomy in Computerized Images; 2005. [37] Ludbrook J: Advantages of permutation (randomization) tests in clinical and experimental pharmacology and physiology. Clinical and Experimental Pharmacology and Physiology 1994, 21(9):673-686, 1994. [38] Burel G, Henocq H: Three-dimensional invariants and their application to object recognition. Signal Process 1995, 45(1):1-22, 1995. 85 [39] Healey DM, Jr., Rockmore DN, Moore SB: An FFT for the 2-sphere and applications. In: 1996; 1996:1323-1326 vol. 1323, 1996. [40] Vranic DV: An improvement of rotation invariant 3D shape descriptor based on functions on concentric spheres. In IEEE International Conference on Image Processing: 2003; Barcelona; 2003:757-760, 2003. [41] Kazhdan M, Funkhouser T, Rusinkiewicz S: Rotation invariant spherical harmonic representation of 3D shape descriptors. In Proceedings of the 2003 Eurographics/ACM SIGGRAPH symposium on Geometry processing. Aachen, Germany: Eurographics Association; 2003. 86 4 Detection and Analysis of Anatomical Shape Changes in MRI (PD brain ventricles) 4.1 Introduction Magnetic resonance imaging (MRI) is increasingly recognized as a suitable modality for studying the human brain anatomy in vivo. Improvements in both MRI scanner hardware and acquisition sequences are yielding progressively higher resolution 3D images of internal brain structures. One of the main benefits of this increased resolution is in the study of structure of brain anatomy. Researchers use MR images to delineate specific brain areas, like the thalamus, to observe structural changes either over time or across two or more subject groups, e.g. patients with Parkinson's disease (PD) and normal volunteers. This process of delineating specific structures to perform targeted analysis is termed region of interest (ROI) analysis. Progress in automated segmentation is increasingly yielding robust and consistent results which are slowly replacing the manual delineation process [1, 2]. However, current methods are suitable mainly to segment regions with clear anatomical boundaries, like the ventricles [3]. Numerous methods have been proposed to quantify and analyze the ROIs. The simplest and most common approach uses the overall volume within an ROI, e.g. Anstey et al A version of this chapter will be submitted for publication. Uthama A, Abugharbieh R, Traboulsee A, McKeown M. J, "Invariant SPHARM Shape Descriptors for Complex Geometry in MR Region of Interest Analysis" 87 used the hippocampal volume to investigate changes observed in mild cognitive impairment (MCI) [4]. However, studies using volumetric measures do not reflect actual shape changes that might be indicative of neurological change. Vetsa et al. used medial representations to characterize the surface of the caudate and identify shape changes in schizophrenia [5]. However, this approach requires accurate alignment of the ROIs across different subjects for valid analysis and hence residual errors in mutual alignment could adversely affect the final results. Alternate methods rely on obtaining descriptive feature vectors which may reflect the structural information in an ROI. These features may be invariant to rotations, translations and scale, thus enabling easy comparison across different subjects. This eliminates the need for mutual alignment, saving time and potential inconsistencies arising from registration errors. Zimring et al [6] used spherical harmonic (SPHARM) based invariant indices to accurately estimate and study volume of brain lesions, reducing the errors due to subject positioning. In an earlier work, we used scale, rotation and translation invariant features derived using a SPHARM representation of 3D ROIs to analyze shape changes in brain MRI data [7] and demonstrated the advantages of using such an approach along with traditional volumetric analysis. However, the parameterization technique used in our previous work did not allow for shapes with non-convex topology. In this paper, we extend our method and propose new enhanced SPHARM shape descriptors that are unique to any ROI. These features also enable the characterization of complex shapes with complex topology, including those with multiple disjoint parts. We validate our method on synthetic data and demonstrate the effectiveness of the proposed technique even in the presence of inter-subject variability and misalignment between subjects. 88 Gerig et al [8] proposed using a SPHARM surface representation for shape of brain structures like the ventricles to better understand neurodevelopment and neurodegenerative changes. Their approach however requires explicit alignment of the structures along with volume normalization and is valid only for structures with spherical topology. The SPHARM features proposed here is not limited to spherical topology and natively incorporates alignment and scale issues as part of its invariant nature. Further studies [9, 10] have also shown the importance of ventricular shape and volume in understanding neurological disorders. In this study, we analyze the shape of lateral ventricles in population of subjects with Parkinson's disease and detect significant changes not detected by the conventional volumetric method. 4.2 Methods In this section, we briefly summarize our earlier work on SPHARM-based ROI descriptors in MRI. We then proceed to present the proposed enhancements used to enable the study of complex ROI topologies including a novel radial transform used to obtain unique features for any ROI. 4.2.1 Proposed Invariant SPHARM features A 3D ROI with spherical topology can be expressed as T(6 ) ,Ø) , a function defined on the unit sphere with 0 and 0 as the zenithal and azimuthal angles respectively. The value of the function at a given (0,0 ) is equal to the distance to the surface point from the centroid at those angles [7]. The SPHARM representation of this function is given by: 89 = fdofY,:(8,0)*,0)sin(e)de.^(4.1) 0^0 V. (0, 0) is the complex conjugate of the mth order spherical harmonic of degree 1, where / ranges from 0 to L [11]. The accuracy of this representation depends on the value of L, also called the bandwidth. Invariance of the representation to translation is obtained by shifting the origin of the 3D ROI to its centroid before computing the SPHARM coefficients. In our previous work [7], we employed rotationally invariant features from this representation using (4.2). Scale normalization was achieved by scaling the vector obtained in (4.2) with its first element, which gives a rough measure of total volume [7]. No= Ecrc*T.^(4.2) Non convex topologies, i.e. surfaces with self occlusions, or topologies with disconnected components cannot be represented as a function of two variables, W(0, 0) , since these topologies do not always have a single distance value for all (0, 0 ). A third variable r, the distance from the origin can be introduced to extend the definition to include such shapes. The augmented function P(r, 0,0) is then binary valued, taking a non-zero value only in points inside the 3D ROI. The SPHARM representation of such a function is given by: r 2 2i^s in (gkr )= Jr dr dO 2^17,77(9,0)T(r,0,0)sin(0)d6). 0^0^o (4.3) 90 k is an index introduced to account for possible degeneracy [11]. Since direct computation of (3) is highly inefficient, we use an alternate approach by representing the data as a set of spherical functions. These functions are obtained by intersecting the 3D binary ROI with R spherical shells of increasing radii r. This intersection is simulated by sampling the 3D ROI on a spherical grid of constant dimensions 2Lx2L at each value of r [12]. The SPHARM coefficients are then obtained at each value of r as: 2ff^fr C71 = PO NO ,0)41(r, 19 9 0)sin(9)cit9. 0^0 (4.4) This results in a SPHARM representation for each shell. Invariant features for these shells can then be derived using (2), as used by Kazhdan et al [13]. However, as the authors themselves note, this representation is not unique, since independent rotations of the shape along any of the constituting shells will result in the same feature vector representation even though the topologies of the shapes are drastically different. To overcome this limitation, we propose adapting the radial transform in (4.3) across these features as in (4.5). This yields unique features even under independent rotations along the shells. To achieve scale invariance across different ROIs, we propose to keep the number of shells, S, constant at 2R„,,,, for all ROIs being analyzed. R„,„, is the radius measure in voxels of the smallest common sphere encompassing each of the ROI in the entire set being analyzed. 91 2 Rrm,„ —r sin gkrcm^,/ 2 ^( )̂kl = Er2 rl • r=1 k^=^ (4.5) Since this radial transform is applied across each value of r, as shown in (4.5), the vector cm must be of the same length for each shell. Hence, we use a common value of L for all shells, whose minimum value is obtained by equating the surface area of the encompassing sphere to the equiangular sampling grid (4.6). 47rR„2,„„ = 2L x 2L, L =^(4.6) The final unique, scale, translation and rotation invariant features are defined as in (4.7). These values are then reshaped into a single row to provide the final feature vector. k=2Rmax m=/ A qp, q) = E^= 1...L,q 1...2Rmax^(4.7) k =1 m=-1 Figure 4.1 shows how this method could be used to analyze 3D ROIs like the human brain ventricles. 92 •Figure 4.1 Intersection of a left brain ventricle with 2 of the 128 shells used to obtain the SPHARM representation. The corresponding sampled surfaces of the two shells are also shown. 4.2.2 Statistical Analysis of Features The feature vectors thus obtained are of a high dimension, dependent on the value of L chosen. Moreover, the generating probabilities of these features are unknown. We thus employ a non-parametric permutation test for the statistical analysis of these features, since no assumptions of the feature parameters are required. In permutation tests, a distance measure between two groups of shapes (represented by their feature vectors) is tested against the distance obtained by all possible permutations 93 of the vectors across the two groups. We use the Euclidian distance between the feature vectors as the distance measure. Also, as reported by Vetsa et al [5] we use a Monte Carlo approach to obtain a sufficiently high number of permutations, since testing against all possible permutations is not feasible. A final p value is obtained as the ratio of the number of permutations which yielded a Euclidian distance greater than the original ordering of the samples and the total number of permutations performed. If the obtained p value is less than a statistical standard threshold of 0.05, the two groups are declared to have significant differences between them. In our study, this is a strong indication that the shapes in the two groups have a systematic difference in structure. 4.2.3 Automated segmentation of the lateral ventricles The segmentation of the lateral ventricles was obtained using the voxel-based morphometry approach [14]. The ventricles were segmented in the subject space without spatial normalization. The left and right ventricles were segmented on the tissue probability map (TPM) of the cerebrospinal fluid. This TPM, along with those for white and grey matter were registered to the T1 anatomical scan of each subject. These registration parameters were then used to obtain the ventricles segmentation in the subject brain space. Automatic 3D connectivity methods and morphological operators were then used to trim the results of segmentation. Figure 4.1 shows the result of once such segmentation. 94 4.3 Data and Imaging This section describes the generation and acquisition of data used to test our proposed approach, which include synthetic validation data as well as real structural MRI data of patients and controls. 4.3.1 Synthetic Structural Data Generation To verify the performance of the proposed shape description technique, we test the SPHARM features on synthetic data with simulated, but nevertheless realistic inter- subject variability. Similar to [6], the basic shape used to generate these data is an ellipsoid as defined in (4.8), where (xc, yc, zc) is the centroid with (a, b, cG) as parameters defining the 3D volume of the ellipsoid (4.8). The advantage of such an experimental setup is in ability to quantitatively control the shape of the object; this lets us gauge the performance of SPHARM features with increasing noise and variability. Keeping the equatorial radii, a and b constant, the polar radius cG can be systematically varied to create shapes with simple differences. The ability of the SPHARM features in discriminating two groups of ellipsoids generated with a small difference in their respective polar radii can be used as a measure of its sensitivity. X 2 y 2^Z 2qqx, y, z)=^+ 2 + 2̂a^b^cG (4.8) To ensure that most sources of variations observed in real life are considered in this experiment, each ellipsoid, generated in a 128x128x128 voxel grid, was put through some pre-processing. First, to create a more realistic irregular surface seen in ROI masks, uniformly distributed random noise was added to the generated ellipsoidal binary mask. 95 (a) (b) This was followed by a 3D Gaussian smoothing with a 5x5x5 voxel kernel with a standard deviation of 4. Finally, the result was converted back into a binary mask using a threshold of 5. These values were obtained heuristically from the observation of real ROIs of the thalami in normal subjects. To simulate the effects of intersubject variability in both size and shape. All the three radii values were perturbed independently by up to 0.2 units, resulting in a maximum voxel count (volume) change of up to 25%. This change is comparable to the variations observed in real data (Section IV-B). To simulate the effects of subject position within the scanner, the centroid for each ellipsoid was randomly perturbed by up to 5 voxels about the common origin fixed at the centre of the voxel grid (64, 64, 64). Further, the ellipsoid was then rotated about a random combination of x, y, and/or z axis by an angle generated from a uniform distribution between 0-45°. Figure 4.2 (a) Surface rendering of a real left thalamus (top) and right thalamus (bottom) from a PD patient, note the rough edges. (b) Group A synthetic samples (see text) generated with parameter set (2, 5, 10). (c) Group B synthetic sample generated with parameter set (2, 5, 8). Note the realistic looking edges in the synthetic data. 96 4.3.2 Real Structural MR Data Acquisition The data used in this study consisted of MRI scans of 23 control subjects and 23 PD patients. Two scans were performed with a gap of two hours, one before and one after the administration of the drug L-dopa. The MR data were acquired on a 3.0 Tesla Siemens scanner (Siemens, Erlangen, Germany) with a birdcage type standard quadrature head coil and an advanced nuclear magnetic resonance echoplanar system. The participant's head was positioned along the canthomeatal line. Foam padding was used to limit head motion within the coil. High-resolution T1 weighted anatomical images (3D SPGR, TR=14ms, TE=7.7 ms, flip angle=25°, voxel dimensions 1.0x1.0x1.0mm, 176x256 voxels, 160 slices) were acquired. 4.4 Results and Discussion This section presents quantitative results of employing the proposed SPHARM shape features on the synthetic and real MR data outlined before. 4.4.1 Validation on Synthetic Data Synthetic data generated as described earlier consisted of forty one data sets each set having two groups, A and B, each group with twenty sample shapes. Group A was re- generated independently each time with a baseline parameter set (2, 5, 10). Group B was generated using a different baseline parameter, starting with (2, 5, 8) and ending with (2, 5, 12), with the value of cb being incremented by 0.1. The proposed SPHARM invariants were then derived for all samples. Rma, was found to be 22; hence 44 shells were used to parameterize the ROIs. A value of 40 was chosen for L as per (4.6), resulting in a feature vector length of 1760 elements. The Euclidian distance between the SPHARM features of 97 0.2 — § 0.8 — .o 0.6 — (.3 co .r11 0.4 — 73 Non Significant * Significant m • ** ** 671 1.2 — 0 8^8.5 9 10^10.5^11^11.5^12 Cb the two groups, A and B, was recorded. Permutation tests were then performed to determine the level of significance in the difference between the features of the two groups. Figure 3 summarizes the results of this analysis. Figure 4.3 . Plot of the Euclidian distance between the two groups A and B for different values of cb, (cc, was fixed at 10 for all points). Distances resulting in a significant difference (alpha = 0.05) between the groups as obtained by the permutation test are indicated with a star. Non-Significant distances are indicated with a dot. The graph in Figure 4.3 shows the expected trend in the Euclidian distance measure between the two groups. For values of cb smaller than 9.5 and larger than 10.6, the features are able to detect a systematic difference in shape between the two groups. However, as the value of cb approaches 10 (the value of ca), there is less difference in the sample shapes between the two groups. For values of c b in the range (9.5, 10.6), the combined effect of the various perturbations in the parameters involved, coupled with the inter-subject variability introduced, leaves almost no difference in constituting shapes of 98 the two groups. This is shown in Figure 3 with dots, indicating that the permutation test was unable to find any systematic difference in the features vectors. This experiment demonstrates the sensitivity of the proposed technique in discriminating subtle shape changes between two groups. 4.4.2 Analysis of Real MR Structural Data The proposed technique was next used to analyze the structure of the left and the right ventricle in patients with Parkinson's disease. Both ventricles were segmented using the automated segmentation method from each scan. Since the subject was scanned twice within a short duration, two independently segmented ROIs were available for each ventricle in each subject. SPHARM features were computed for each of the ROI masks so obtained. R„,,,, was found to be 64 voxels; hence 128 shells were used to parameterize the ROIs. The bandwidth L was set at 114 as per (4.6), resulting in a feature vector length of 14592. Since two scans were performed for each subject, two feature vectors were obtained for each ROI, for each subject. Also, since the duration between the scans was less than two hours, no systematic structural changes were expected between these two observations. The only changes expected were due to repositioning, manual tracing inconsistencies and other imaging artifacts. To confirm this assumption, each subject group (Patients and Controls) was further sub-divided into two categories: pre-drug and post-drug. Permutation test results of the SPHARM feature vectors between these two sub-groups are shown in Table 4.1. Results from the volumetric analysis are also shown. 99 Table 4.1 Result of pre-drug vs. post-drug analysis Subject group ROI SPHARMp value Volumetric p value Controls Left Ventricle 0.9803 0.9563 Controls Right Ventricle 0.9610 0.9710 PD Left Ventricle 0.9838 0.9638 PD Right Ventricle 0.9800 0.9781 The results in Table 4.1 indicate that, as expected, there are no statistically significant changes in the shape of the ventricle between the two scans. Having established this, the two feature vectors from the two sessions were then averaged to generate a single feature vector for each ROI for each subject. Direct averaging is valid since these features represent the same shape. Changes due to rotation and translations are annulled by the invariance property of the features. This averaging reduces any potential variability introduced by the automatic segmentation and other noise sources, thereby increasing the sensitivity of the features to actual shape changes. A permutation test on the averaged features was then performed between the two subject groups. Results are indicated in Table 4.2. Results from the volumetric analysis are also shown. Table 4.2 Result of comparison between control subjects and PD patients. Values indicated with a * are statistically significant at an alpha value of 0.05. ROI SPHARMp value Volumetric p value Left Ventricle 0.0462* 0.2711 Right Ventricle 0.0090* 0.2837 100 As the ventricles may increase in size as a result of brain atrophy (so called "hydrocephalus ex vacuo"), it remains to be seen if different shape changes correspond to different aspects of PD such as dementia. 4.5 Conclusions We enhanced our earlier work generating invariant SPHARM features for 3D ROI analysis in MR data by extending our characterization capability to ROIs with complex topologies, including those with possible disconnections. We also proposed the use of a novel radial transform to obtain unique feature. The proposed features were able to detect subtle shape changes in synthetically created data with realistic inter-subject variations. Given the fact that these measures can be obtained non-invasively and serially this method provides another means by which clinicians may monitor progression of disease, and possibly test whether or not a given treatment affects the rate of disease progression. Hence, the use of the proposed method to explore shape changes of different brain structures as a biomarker of disease progression warrants further study. 101 4.6 References [1] J. S. Suri, S. Singh and L. Reden, "Computer Vision and Pattern Recognition Techniques for 2-D and 3-D MR Cerebral Cortical Segmentation (Part I): A State-of-the- Art Review," Pattern Analysis & Applications, vol. 5, pp. 46-76, 2002. [2] J. S. Suri, S. Singh and L. Reden, "Fusion of Region and Boundary/Surface-Based Computer Vision and Pattern Recognition Techniques for 2-D and 3-D MR Cerebral Cortical Segmentation (Part-II): A State-of-the-Art Review," Pattern Analysis & Applications, vol. 5, pp. 77-98, 2002. [3] M. G. Linguraru, M. Ballester and N. Ayache, "Deformable Atlases for the Segmentation of Internal Brain Nuclei in Magnetic Resonance Imaging," International Journal of Computers, Communications & Control, vol. 2, pp. 26-36, 2007. [4] K. J. Anstey and J. J. Mailer, "The role of volumetric MRI in understanding mild cognitive impairment and similar classifications," Aging Ment. Health., vol. 7, pp. 238- 250, Jul. 2003. [5] Y. S. K. Vetsa, M. Styner, S. M. Pizer, J. A. Lieberman and G. E. -. Gerig, Caudate Shape Discrimination in Schizophrenia using Template-Free Non-Parametric Tests. , vol. 2879, 2003, pp. 661-669. [6] D. Goldberg-Zimring, A. Achiron, S. K. Warfield, C. R. G. Guttmann and H. Azhari, "Application of spherical harmonics derived space rotation invariant indices to the 102 analysis of multiple sclerosis lesions' geometry by MRI," Magn. Reson. Imaging, vol. 22, pp. 815-825, 7. 2004. [7] S. Tootoonian, R. Abugharbieh, X. Huang and M. J. McKeown, "Shape vs. volume: Invariant shape descriptors for 3D region of interest characterization in MRI," in ISBI, 2006, pp. 754-757. [8] G. Gerig, M. Styner, D. Jones, D. Weinberger and J. Lieberman, "Shape analysis of brain ventricles using SPHARM," in 2001, pp. 171-178. [9] M. Styner, J. A. Lieberman, R. K. McClure, D. R. Weinberger, D. W. Jones and G. Gerig, "Morphometric analysis of lateral ventricles in schizophrenia and healthy controls regarding genetic and disease-specific factors," Proceedings of the National Academy of Sciences, vol. 102, pp. 4872-4877, 2005. [10] X. Huang, Y. Z. Lee, M. McKeown, G. Gerig, H. Gu, W. Lin, M. M. Lewis, S. Ford, A. I. TrOster, D. R. Weinberger and M. Styner, "Asymmetrical ventricular enlargement in Parkinson's disease," Movement Disorders, vol. 22, pp. 1657-1660, 2007. [11] G. Burel and H. Henocq, "Three-dimensional invariants and their application to object recognition," Signal Process, vol. 45, pp. 1-22, 1995. [12] D. M. Healey Jr, D. N. Rockmore and S. B. Moore, "An FFT for the 2-sphere and applications," in 1996, pp. 1323-1326 vol. 3. 103 [13] M. Kazhdan, T. Funkhouser and S. Rusinkiewicz, "Rotation invariant spherical harmonic representation of 3D shape descriptors," in Symposium on Geometry Processing, 23-25 June 2003, 2003, pp. 156-65. [14] J. Ashburner and K. J. Friston, "Voxel-based morphometry-the methods," Neuroimage, vol. 11, pp. 805-821, 2000. 104 5 Extending SPHARM Features to Functional Activation Analysis in fMRI 5.1 Introduction Functional magnetic resonance imaging (fMRI) is one of the most commonly used modalities for human brain mapping. Medical researchers often use this non-invasive imaging technology to investigate normal brain functions as well as the effects of neurological disease such as Parkinson's (PD) [1]. Typically, the functional response obtained in an fMRI experiment is analyzed on a voxel-by-voxel basis with the resultant activation statistics assembled into a statistical parameter map (SPM). The most common approach in generating an SPM is based on the general linear model (GLM) [2]. A key challenge in functional neuroimaging is determining how to meaningfully combine results across subjects with varying brain sizes, brain shapes, and possibly even disease- induced brain shape changes. A typical approach is to register each subject's brain to a common atlas and compute the SPMs in the template space. A comparison of some of the spatial normalization methods was presented by Crivello et al. [3]. However, current spatial normalization methods may give an imperfect registration result [4], resulting in signals from functionally distinct areas, especially small subcortical structures, to be A version of this chapter will be submitted for publication. Uthama A, Abugharbieh R, Traboulsee A, McKeown MJ, "3D Regional fMRI Activation Analysis Using SPHARM-Based Invariant Spatial Features and Intersubject Anatomical Variability Minimization", 105 inappropriately combined [5]. Spatial normalization may therefore lead to poor sensitivity in the fMRI data analysis due to reduced functional overlap across subjects as was observed by Stark and Okado [6]. Extensive spatial smoothing is often performed to increase the functional overlap across possibly misaligned subjects, but this procedure degrades the overall spatial resolution. An alternative approach that was recently investigated is to align the subjects at the region of interest (ROI) level, as opposed to whole brain normalization. Miller et al. [7] extended the work of Stark and Okado [6] by utilizing large deformation diffeomorphic metric mappings (LDDMM). Also, an approach that uses continuous medial representations (cm-rep) of the ROIs has also been proposed [8]. However, these approaches still rely on structural features in the higher resolution scans to deform the functional maps into a common space. This inherently assumes a high level of correspondence between the structural and functional neuroanatomy across subjects. In the presence of brain pathologies like PD, which are known to induce compensatory changes affecting, this assumption may become problematic. To circumvent these potential problems associated with normalization, many researchers employ direct ROI analysis without using any normalization [9]. This results in enhanced sensitivity to detection of activation [10]. However, to date, most direct ROI analysis methods ignore information encoded by the spatial distribution of the activation statistics and instead use simpler measures such as mean voxel statistics or percentage of active voxels within an ROI as features [9], which ignore information encoded by the spatial distribution of the activation statistics. 106 An attempt in incorporating spatial information has been proposed [11], where the sums of activation statistics within spheres of increasing radii were considered. However, these features have limited sensitivity to spatial patterns since they only capture changes in the radial direction. Ng et al. [12] proposed an alternate approach employing three dimensional moment invariant features (3DMI) which were shown to be more sensitive than conventional mean voxel statistics and percentage of activated voxels based approaches in detecting task-related activation differences. The fact that the 3DMI approach detected significant task-dependent features (e.g. spatial variance) across non- normalized brain images is noteworthy, as it implies that the spatial features are sensitive to task-related effects and robust to inter-subject anatomical differences. Nevertheless in its current form, the 3DMI method will have reduced sensitivity to task-related functional changes if inter-subject anatomical variability is large, especially when analyzing older subjects or subjects with neurodegenerative disease. Also, as the authors noted, higher order 3DMI features are more susceptible to noise and hence limit the maximum number of discriminative features that can be derived [12]. SPHARM-based methods have previously been proposed in the context of 3D shape retrieval systems by Vranic et al [13] and Kazhdan et al [14]. Both authors proposed obtaining invariant SPHARM features by intersecting 3D shapes with shells of growing radii. This approach, however, could not detect independent rotations of a shape along the shells, thereby resulting in a non-unique representation [14]. In [15], the use of invariant SPHARM descriptors for analyzing anatomical structures (represented as binary 3D shapes) in MR was proposed. However, the shape parameterization employed was limited to convex topologies, which significantly limits its applicability. However, these 107 approaches were proposed to analyze the geometrical shape of an ROI, but not the intensity distributions of an ROI, such as the spatial activation patterns in fMRI data. In [20], we proposed a radial transform adapted from [16] to address the limitations of [13, 14]. This resulted in a novel SPHARM-based structural feature representation that is unique to the ROI and can be applied to any arbitrarily complex shape. We then extend this application of 3D SPHARM-based features to characterize spatial activation patterns in fMRI data. Importantly, these proposed features are unique to the ROI and invariant to similarity transformations, so alignment of brains (or ROIs) across subjects is not needed. Our key contribution is a novel approach to account for the underlying structural (anatomical) inter-subject variability based on subspace projection methods. We validate our proposed technique on synthetic data and demonstrate improved sensitivity as compared to the 3DMI technique. We also demonstrate the effectiveness of the proposed method on real fMRI data by discriminating medication (L-dopa) induced differences in activation patterns of Parkinson's disease (PD) patients. 5.2 Methods In this section we first describe our proposed invariant SPI-IARM-based spatial features. We then present the PCA subspace approach for minimizing effects of structural ROI variability on functional analysis. Finally, we describe the statistical analysis method used to generate the results. 108 5.2.1 Invariant SPHARM-Based Spatial Features Let P(9,0) be a function defined on the unit sphere with 8 and 0 as the zenithal and azimuthal angles, respectively. The SPHARM representation for this function is given by (4.1) where 17 : (0, 0) is the complex conjugate of the M th order spherical harmonic of degree 1. 1 ranges from 0 to L [16]. Increasing the value of L, also called the bandwidth, improves the representation accuracy at the cost of higher computation time. Burel and Henocq [16] also extended this definition to real valued 3D distributions '1+,0,0) (5.2), where r is the distance from the origin to a given voxel. k is an index introduced to account for possible degeneracy due to the additional dimension [16]. c; = 'do fY,:i (0, OW, 0)sin(0)d0 0^0 ^ (5.1) C kmi =ir^2ff^fr^. ( 7ricrl 'dr c10 fli sink I Yi* 0, OS* , 6) ,O)sinHc10 0^0^0 r^in (5.2) As explained later in this section, rotationally invariant features can be derived from this spherical harmonic representation. In our application, we also need the features to be invariant to any translation of the entire ROI in 3D space. To achieve this, we move the origin of the function '1+,0,0), to the centroid of Ts (r,0,0), where Tjr,9,0) is given by (5.3). 109 fr , {1 if Tfr,t9,0) ^ 0, 9,01= o if kilr,9,0)=-- 0 (5.3) Since direct computation of (5.2) is highly inefficient [17], we use an alternate approach by representing the data as a set of spherical functions obtained by intersecting the 3D data with spherical shells. Alternatively, for each value of r, '1+,0,0) can be visualized as a spherical shell comprising the function values at a distance r from the origin. r can then be incremented in steps of t to encompass the entire ROI. If the initial representation of the function is in the form of a cubic grid (regular isotropic voxels in our case), volumetric interpolation is required to resample the ROI in the spherical coordinate space. When analyzing multiple subjects' functions simultaneously (e.g. fMRI activation data for an ROI across a subject group), we define the maximum radius, Rmax, as the minimum radial distance in voxel count that encompasses all non-zero values of all subjects' functions being analyzed. To represent the values from the cubic grid of all functions with sufficient accuracy, 2R„,,,, shells are used. To achieve scale invariance, the shells must be distributed evenly throughout the spatial extent of each function. Since the ROI size across subjects is not uniform, shell spacing t must be adjusted for each subject separately. This procedure ensures that each shell captures similar features from the 3D function irrespective of its scale. Surface sampling along each of these shells is performed on an equiangular spherical grid of dimensions 2Lx2L [17]. The common bandwidth L for all shells of all functions is chosen to satisfy the sampling criterion for the largest shell in this set of functions, 110 namely the one with radius 1?,,,,,x. The surface area for this shell represents the maximum surface shell area that needs to be sampled by the equiangular grid; hence, any value of L satisfying the required equiangular sampling (2L x 2L) at this shell will be sufficient to represent data from smaller radii shells. The minimum value for L is obtained by equating the surface area of this largest shell to the equiangular sampling grid (5.4). Higher values of L are not used, since the computation time will increase. Also, this will result in longer feature vectors, complicating the analysis. Furthermore, Burel and Henocq [16] noted that when the represented object is a discrete array, higher values of 1, resulting from a larger L, may correspond to sampling noise. Recognizing that in applications pertaining to discrimination, high accuracy in the SPHARM representation is not a necessity, we chose to use the minimum value for L as that obtained by (5.4). 42 -RL„ = 2L x 2L, L -= Rn,a„.Fr (5.4) To obtain the SPHARM representation for all shells, a discrete SPHARM transform is performed at each value of r to obtain cm, (5.5). Features derived from this representation, however, do not provide a unique function representation as discussed in [13] and [14]. For instance, rotating the inner and outer shells by different amount will result in different spatial distribution of function values. However, in this approach, the derived features are insensitive to these rotational transforms, thus resulting in the same feature values for dissimilar spatial distributions. 111 2ir cm^id0 SY1,;(0 0) - (r It 9 ,0)Sin(e)d0 0^0 r =[1 , 2 , 3 , ... , 2Rm ]^(5.5) Burel and Henocq's original equation (5.2) does not have this problem, since a part of the basis function is a function of r. However, since (5.2) is computationally intractable [17], we earlier proposed an efficient approach that uses a radial transform (5.6), derived from (5.2), to obtain a unique function representation [20]. The transform (5.6) retains the relative orientation information of the shells, thus the features derived will be sensitive to independent rotations of the different shells, thereby ensuring that unique feature representation is obtained. m =^r 2^sin(gkr) c m r c kl rl r=1 k = [1,2,3,...,2R.]^(5.6) The range of k could be changed to obtain different lengths of the final feature vector. However, to avoid unnecessarily increasing the feature vector length or losing any important information caused by reducing the range, we choose to keep the range of k the same as that of r, i.e. 2R,nar From the obtained representation (5.6), we then compute similarity transform invariant features using (5.7) for each value of I and k as outlined in [16]. p and q are used to index 112 these features. Note we reshape I into a single row vector of dimensions D = L x 2R,,. for later analysis. k=2R,„a, Ap,q)= E E (c:)*, p =1...L,q = /...2Rmax k=1 m=-1 (5.7) 5.2.2 fMRI Analysis Using SPHARM Features Features extracted from fMRI activation statistics are influenced by two factors. The first factor, which is the one of interest, corresponds to the functional aspect of the data as exhibited by the spatial pattern of the activation statistics within the ROI. The second factor is the shape of the anatomically-defined ROIs. In functional studies, the structural information reflects pure inter-subject variability, which adversely affects the primary aim of the analysis, namely resolving actual functional pattern changes across subjects. In this section, we propose a novel approach for fMRI ROI characterization that uses SPHARM-based features described in the previous section in a manner that minimizes the effects of the underlying structural (inter-subject) variability. In our approach, invariant features describing the shape of an ROI, SPHARM-s, are first computed without incorporating the fMRI activation statistics within the ROI (5.8). We then compute invariant features comprising of both structural and functional information, SPHARM-fs, as in (5.9), where txyz is the real valued activation statistic obtained at the voxel position (x, y, z). The difference between SPHARM-s and SPHARM-fs is that SPHARM-s contains only shape information, whereas SPHARM-fs contain both structural and functional information. We note that variations in ROI shapes thus affect both vectors, SPHARM-f:s. and SP HARM-s. 113 {1 Within the ROI R01 (x, Y z) =^ SPHARM—s0 Elsewhere (5.8) 'I' nal (X, y,z) , xy, WithintheROI SPHARM — 0 Elsewhere (5.9) To quantify the effects of inter-subject structural variability, we first derive the SPHARM- s features for all subjects for a given ROI. Let I S be a matrix whose ith row corresponds to the SPHARM-s features of the subject. Each of these feature vectors will have a length D as explained in section 5.2.1. We apply principal component analysis (PCA) on the pooled structural information (of all subjects) in IS to determine the projection directions of maximal variance. Means of each column of is are first removed. A singular value decomposition of IS is then applied to obtain a DxD matrix V and a diagonal A, where the columns of V (eigenvectors) {vv} represents the principal component directions, and the diagonal elements of A {ii } (eigenvalues) represents the variance of the corresponding eigenvector. V spans a linear subspace where the first several projection directions {vi } reflect maximum variations in the row elements of L. Since i s contains the pooled structural information for all subjects, these directions correspond to maximum inter- subject structural variations. For the data used in this study, we observed that the first d (3 to 4) components of V account for up to 98% of the inter-subject variation. To mitigate inter-subject variability, we first define a new set of projections V' (5.10) by removing the first d components for V. The resulting subspace V' thus represents the set of projections minimally effected by the inter-subject variation in ROI shape observed in the SPHARM-s features across both groups. 114 V', {1, j d,...,D (5.10) To minimize the effects of the structural variations on functional analysis, let /f, be the matrix obtained by deriving the SPHARM-fs features for all subjects. The columns of If, , {fSj , are first mean centered and then projected onto the linear subspace V' (5.11). If = Ifir (5.11) We denote features resulting from this projection as SPHARM-f features. The obtained SPHARM-f features are minimally affected by structural variations, and thus ensure that subtle functional pattern changes will not be obscured by structural inter-subject variability. Thus, discriminability in detecting activation pattern differences is enhanced as will be demonstrated in Section 5.4.1. This approach is similar to methods used in the face recognition application, where weighted PCA subspaces reduce the sensitivity of different facial expressions within the same subject [18]. We validate our approach in section 5.4.1 on synthetic data, and show that the proposed SPHARM f features are more sensitive to functional pattern changes. Results indicate that for a given level of inter-subject structural variations, SPHARM-f features perform remarkably well even under high levels of noise. We demonstrate the practical use of our approach by analyzing the effects of drug on the functional activation patterns of PD patients. 5.2.3 Statistical Group Analysis To efficiently discriminate two groups of 3D fMRI distributions, we first derive SPHARM-f features for each subject's ROI as explained in section 5.2.2. Then, to 115 determine the level of statistical significance in the difference between the groups, we use the non-parametric permutation test as reported by Vesta et al [19]. This approach uses the mean Euclidian distance between the feature vectors of the two groups to compute a p value. A p value below the statistical standard threshold of alpha=0.05 is considered as an indication of significant difference existing between the classes. This test is well suited for our long feature vectors whose generating probability distributions are not known. 5.3 Data and Imaging This section describes the synthetic fMRI data we generated to validate the proposed fMRI spatial analysis technique. It also outlines the imaging procedure used to obtain real fMRI data from PD patients. 5.3.1 Synthetic fMRI Data To demonstrate the validity of the proposed spatial analysis method, synthetic fMRI data sets were generated using ROIs obtained from real MR data with actual inter-subject anatomical variability. Eighteen binary ROI masks of the right cerebellar hemisphere were obtained during two MR sessions of nine PD patients (imaging details presented in Section 5.3.2). The same subjects took part in the two sessions, which were an hour apart hence no significant structural changes were expected between them. However, since the ROIs were manually delineated twice, i.e. on each subject's image data pair from the two sessions, inter-scan variability in the delineated shape were also expected to be present. These 18 ROIs were then randomly assigned to two groups, A and B (9 each) to create a data set with realistic inter-subject variability in ROI structure (shape). To simulate affects of subject positioning in the scanner, these ROIs were then randomly rotated up to 116 an angle of ten degrees about a random combination of the three spatial axes. In an approach similar to that used in [11] we generated synthetic functional patterns within the binary ROI masks resulting in 3D SPMs. Since the aim of this experiment is to demonstrate the sensitivity of the proposed features in detecting functional patterns, we create group A with no patterns (pure noise) and group B with a simple, known, quantifiable pattern. This pattern, for group B ROIs, was created using two levels of baseline activation; -delta and +delta were used as activation values depending on the distance from the functional centroid. The functional centroid for each subject was obtained by a random perturbation of the ROI centroid by up to three voxels in all three spatial directions. This was to simulate possible inter-subject variability in the functional centroid in real data. —delta valued voxels were generated up to a distance of 4 voxels from the functional centroid while for distances greater than 4 voxels, +delta valued voxels were used. Considerable amount of research has gone into examining the noise model in the fMRI time series, yet no clear noise model exists for the spatial domain in the final parameter maps. Kontos and Megalooikonomou choose to use a probabilistic model to generate the pattern, but they do not add additional noise to the data [11]. Here, we choose to generalize all sources of noise to a Gaussian model. To quantify of the level of noise when compared to the signal amplitude present, we use the contrast to noise ratio (CNR) defined by (5.12). 117 CNR ^Asignal^21delta^ (5.12) o The aggregation of these two groups of functional data thus created is termed a set. Multiple such sets were created for the actual experiment as explained further in Section 5.1 by varying the delta value. It is important to note that creation of each set begins from the random assignment of the ROI pool into the two groups; hence the sets would be considerably different from each other. Figure 5.1 shows the cross section of a real fMRI volume taken from the right cerebellar hemisphere of a PD patient (data details explained in Section 5.3.2). Figure 5.2 shows some of the generated patterns in group B for different values of delta. Figure 5.3 shows a cross section of the generated synthetic fMRI activation values for a sample from group A. As explained earlier, group A samples do not have any patterns in their activation values and represent only noise. Slice number 6 15 20 25 15 20^25 30x Figure 5.1 A cross section of fMRI activation statistics from the right cerebellar hemisphere of a real subject. 118 all subsequent movements were scaled to this, they had to squeeze at 5-15% of MVC to accomplish the task. Using the squeeze bulb, subjects were required to control the width of an "inflatable ring" (shown as a black horizontal bar on the screen) in order to keep the ring within an undulating pathway without scraping the sides. The pathway used a block design with sinusoidal sections in two different frequencies (0.25 and 0.75 Hz) in a pseudo-random order, and straight parts in between where the subjects had to keep a constant force of 10% of MVC. The frequencies were chosen based on prior findings and pilot studies were used to determine that PD subjects could comfortably perform the required task. Each block lasted 20 seconds (exactly 10 TR intervals), alternating a sinusoid, constant force, sinusoid and so on to a total of 4 minutes. Before the first scanning session, subjects practiced the task at each frequency until errors stabilized and they were familiar with the task requirements. Custom Matlab software (The Mathworks, Natick, MA) and the Psychtoolbox [24] was used to design, present stimuli, and to collect behavioral data from the response devices. 6.3.2.2 Data Acquisition Functional MRI was conducted on a Philips Achieva 3.0 T scanner (Philips, Best, the Netherlands) equipped with a head-coil. Echo-planar (EPI) T2*-weighted images with blood oxygenation level-dependent (BOLD) contrast were acquired. Scanning parameters were: repetition time 2000 ms, echo time 3.7, flip angle 90°, field of view (FOV) 240.00 mm, matrix size = 128 x 128, pixel size 1.9 x 1.9 mm. Each functional run lasted 260 seconds. Thirty-six axial slices of 3mm thickness were collected in each volume, with a 145 gap thickness of 1 mm. Slices were selected to cover the dorsal surface of the brain and include the cerebellum ventrally. A high resolution, 3-dimensional T1-weighted image consisting of 170 axial slices was acquired of the whole brain to facilitate anatomical localization of activation for each subject. 6.3.2.3 fMRI Data Analysis The functional MRI data were pre-processed for each subject, using Brain Voyager (Brain Innovations, the Netherlands, www.brainvoyager.com ) trilinear interpolation for 3D motion correction and Sinc interpolation for slice time correction. No temporal or spatial smoothing was performed on the data. The data were then further motion corrected with MCICA, a computationally expensive but highly accurate method for motion correction [25]. The Brain Extraction Tool in MRIcro (http://www.sph.sc.edu/comd/rorden/mricro.html) was used to strip the skull from the anatomical scan and first functional image from each run to enable a more accurate alignment of the functional and anatomical images. Custom scripts in Amira software (Mercury Computer Systems, San Diego, USA) were used to co-register the anatomical and functional images. Manual segmentations to obtain ROIs were based on anatomical landmarks and guided by a neurological atlas [26]. Sixteen ROI's, hypothesised to be involved in motor tasks, were drawn separately in each hemisphere. They were outlined on the unwarped, aligned structural scan for each subject using Amira software. The ROIS are: primary motor cortex (M1) (Brodman Area 4), supplementary motor cortex (SMA) (Brodman Area 6), prefrontal cortex (PFC) (Brodman Area 9 and 10), caudate (CAU), putamen (PUT), 146 thalamus (THA), cerebellum (CER) and anterior cingulate cortex (ACC) (Brodman Area 28 and 32). The labels on the segmented anatomical scans were resliced at the fMRI resolution. The raw time courses of the voxels within each ROI were then extracted. Since the task frequencies were too fast to be directly measured with the sluggish Hemodynamic response inherent in BOLD IMRI signals, a block design was used for analysis. A hybrid Independent Component Analysis (ICA) / General Linear Model scheme was used to contrast each of the 2-frequency blocks with the static force blocks [27] and create statistical parametric maps (SPMs). 6.4 Results and Discussion This section presents the results of applying our proposed fMRI spatial analysis technique on the synthetic and real data described in Section 6.3. 6.4.1 Validation on Synthetic Data 6.4.1.1 False Positives as a Function of Bandwidth The bandwidth, L, is the main parameter of the proposed SPHARM approach. The first synthetic example was designed to analyze effects of varying it on the false positive rate of the proposed method. Two samples of TsYnft,m/ were obtained with no signal added (only noise). SPHARM-s features were derived for each of the ten binary ROI masks in PRO/ SPHARM-fs and SPHARM-f features were then derived for each of the 10 synthetic ROI fMRI activation statistics for the two samples of 'V" jA,m/ . These features were then analyzed as explained in section 6.2 to obtain a p value. Since R„,„„ for these sets was found to be 19 voxels, the number of shells used was set at 40 (6.7). The bandwidth used 147 ■ 6^12^18^24^30 36^42^48 LIT 4- 0 C 0 7 3 to obtain all the SPHARM features was varied over a range and at each value, 100 such experiments were performed. Knowing that no difference exists between the two groups, we were able to decide if the resulting p value from the permutation test was a false positive or not. a was set at 0.05 and any p value below this threshold was declared a false positive. The graph so obtained between the false positive count and the bandwidth used is shown in Figure 6.3. Number of trials = 100 25 coa) 17- 0 Figure 6.3 False positive rates using SPHARM-f features for various values of the parameter L (Bandwidth). a = 0.05 . Figure 6.3 clearly shows the stabilization of the false positive rate (FPR) well before the theoretical bandwidth of 36 is reached. A FPR of 3 in a 100 trials is below the expected 5 counts at the set alpha value, confirming FPR control. 148 6.4.1.2 Comparison with ROI-N For this example, two samples of Ts''"ft,,R/ were created, One with no signal (only noise) and another with a known non-zero signal value, SNR was controlled by the magnitude of delta (Section 6.3.1). Figure 6.4 plots the average SNR value in the signal group against the delta value. Results for both the SPHARM and the ROI-N approach were recorded for increasing SNR values. The results are depicted in Figure 6.5 and Figure 6.6. 0.5^ .5 ^ 25 delta Figure 6.4 The average SNR plotted against increasing delta value. Figure 6.5 summarizes the performance of the two SPHARM approaches. SPHARM-fs features are able to discern a difference in the two groups at -3.4402 dB and higher SNR values. SP HARM-f features perform better, picking up the difference at an SNR value of - 10.7674 dB. 149 -10.7674^-3.4402 7.7180 1 , 6X iv * SPHARM-f SPHARM-fs SNR (dB) Figure 6.5 Sensitivity of the SPHARM approach. The mean Euclidean distance between the feature vectors of a group with only noise and of another group with increasing SNR is plotted against the SNR. The star and triangle markers denote distance which yielded a p value less 0.05, signifying that the method detected differences between the groups. Figure 6.6 presents the results of the ROI normalization approach for different height thresholds for the absolute value of the t-statistic. Hayasaka et al [23] have shown that results are sensitive to the height threshold and cluster size threshold used, both of which are dependent on the estimated smoothness of the data. Based on their simulation results, we use a very liberal cluster size threshold of 4 connected voxels (in a 26 voxel neighborhood) corresponding to zero smoothness. A height threshold of 2 is also considered liberal since the corresponding uncorrected (for multiple comparisons) p value is 0.0304 with 18 degrees of freedom. 150 6.4.1.3 Discussion This example clearly shows that for the realistic pattern used, the sensitivity of both the SPHARM approaches outperforms the ROI-N method. SPHARM-f features perform better than SPHARM-fs features indicating the advantage of the principal component subspace approach in reducing intersubject structural variations. The graphs for both SPHARM methods show a consistent trend emphasizing the robustness of the method to noise. The ROI-N method on the other hand, shows a higher susceptibility to noise, especially by the t>131 graph in Figure 6.6. 200-  lOt >121* t >13Ia t >141, 150 0 100 0 50- 25 -3.4402 ^ 0.4522 ^ 4.2688 ^ 6.8817 SNR (dB) Figure 6.6 Sensitivity of the ROI-normalization approach. The number of voxels surviving various height thresholds is plotted against the SNR. A cluster size threshold of 4 connected voxels was used in each case. These results were computed from the exact same data generated for Figure 6.5. 151 6.4.2 Application to Clinical fMRI Data Section 6.3.2 presented the real data that was analyzed using both the SPHARM approaches (SPHARM-fs and SPHAR11/1-fi and the ROI normalization approach. 6.4.2.1 Results of SPHARM analysis SPHARM-fs and SPHARM-f features were derived for all participants as outlined in section 6.2. Table 6.1 presents the maximum radius and bandwidth used for each pair (left and right) of the Sixteen ROIs analyzed. Values for each ROI were obtained considering all twenty participants. All SPHARM-f features were obtained using a subspace with d=3 (10) of the SPHARM-s PC subspace. Table 6.1 Maximum radius of ROIs and the corresponding bandwidth used ROI RmAx Bandwidth ACC 16 30 CAU 11 20 CER 19 36 M1 15 28 PFC 16 30 PUT 6 12 SMA 16 30 THA 8 16 Table 6.2 displays the result of the SPHARM analysis. For each ROI, permutation tests were performed for both SPHARM features comparing the two tasks performed at different frequencies. Each p value indicates the result of the permutation test (N=10,000) with an alpha value a = 0.05 . The table only shows ROIs for which at least one test was 152 found to be significant. ROIs in the left hemisphere have a `L: prefix, while those in the right hemisphere have a `R_' prefix. Table 6.2 SPHARM analysis of the two frequency tasks. Only ROIS with at least one significant result are shown ROI SPHARM-fs Normal SPHARM-f Normal SPHARM-fs PD SPHARM-f PD L CER 0.0613 0.4694 0.0239* 0.0202* L_M1 0.1767 0.1761 0.0863 0.0104* L_PFC 0.5302 0.9497 0.1045 0.0378* L THA 0.1200 0.0425* 0.0651 0.2513 R CER 0.0546 0.1599 0.0478* 0.0309* R_PFC 0.6938 0.5273 0.0150* 0.0166* R_THA 0.0914 0.1057 0.1602 0.7946 6.4.2.2 Results from ROI Normalization The Sixteen ROIs from the twenty subjects were also analyzed using the ROI- normalization approach. Both tasks were performed in the same session enabling the same anatomical Ti ROI to be used to obtain corresponding functional maps. Hence the normalization template obtained was used to warp both tasks fMRI activation statistics (Section 6.2.4). Due to the uncertainty regarding the best height and cluster size threshold to be used, we used a range of height thresholds varying from 2 to 5 in steps of 0.5 in t- statistic magnitudes, while using a cluster size threshold of 4. Among the thirty two combinations of tests performed (2 subject groups with 16 ROIs each), only the left and the right cerebellar hemispheres in the PD subject group had at least one cluster left after the thresholding. 153 6.4.2.3 Discussion These results demonstrate that using the standard ROI normalization approach, increased movement rate was associated with an increase in activity of the cerebellum bilaterally in PD subjects, but not in age-matched, healthy controls. This is consistent with an increased reliance on visual feedback in PD subjects mediated by the cerebellum [28]. Using the SPHARM-fs features also identified a change in the shape of activation within the bilateral cerebellar hemispheres of PD subjects, and in addition was able to detect changes in the activation of the ipsilateral prefrontal cortex of PD subjects. This region is associated with performance monitoring [29], thus a change in the activity of this area may reflect an increased attentional demand in PD subjects as the movement becomes more difficult. The SPHARM-f approach, in which structural variations in ROI anatomy across subjects are minimized, sowed greater sensitivity to changes in BOLD activation features during increased movement speeds. Specifically, control subjects showed a change in the activation of contralateral thalamus with increasing movement rate, consistent with an increased output from the basal ganglia. PD subjects were apparently not able to adjust the output through the thalamus, and instead showed adjustments in the output from bilateral cerebellum, bilateral prefrontal cortex, and contralateral M1 in response to increased task demands. Increased activity of MI and cerebellum in PD has previously been suggested to be a compensatory mechanism [30, 31]. It is thus possible that compensatory mechanisms are recruited during tasks which make greater demands on the motor system, and rely on changes to the shape as well as the level of activation. 154 6.5 Conclusions In this paper, we proposed a new technique for analyzing spatial activation patterns in fMRI data using 3D SPI-IARM features. These features were obtained by intersecting 3D data distributions with concentric shells of increasing radii followed by a radial transform. Our other main contribution was a novel approach to account for anatomical shape variations across subjects while discriminating subtle changes in spatial distribution of fMRI activation patterns within an ROI. The effectiveness of the proposed approach in mitigating structural variability, robustness to noise and comparative sensitivity were validated on synthetic data. We also applied this approach to discriminate functional pattern changes within anatomical ROIs in fMRI data collected. Using this approach, we were able to demonstrate differences in the way that PD subjects and healthy controls respond to an increased task demand, reflecting failure of PD subjects to increase basal ganglia output, and a reliance on cerebellar and cortical activity to enable successful performance. This adjustment may reflect a compensatory mechanism in PD subjects. An interesting application of this approach would be to compare activation patterns between two subject groups that are expected to have systematic changes in both structural and functional aspects. Avenues to decouple these aspects in the feature domain would yield promising new directions in fMRI data analysis. 155 6.6 References [1] K. J. Friston, A. P. Holmes, K. J. Worsley, J. B. Poline, C. Frith and R. S. J. Frackowiak, "Statistical Parametric Maps in Functional Imaging: A General Linear Approach," Human Brain Mapping, vol. 2, pp. 189-210, 1995. [2] F. Crivello, T. Schormann, N. Tzourio-Mazoyer, P. E. Roland, K. Zilles and B. M. Mazoyer, "Comparison of spatial normalization procedures and their impact on functional maps," Hum. Brain Mapp. (USA), vol. 16, pp. 228, 2002. [3] O. T. Carmichael, H. A. Aizenstein, S. W. Davis, J. T. Becker, P. M. Thompson, C. C. Meltzer and Y. Liu, "Atlas-based hippocampus segmentation in Alzheimer's disease and mild cognitive impairment," Neuroimage, vol. 27, pp. 979-990, 2005. [4] M. Ozcan, U. Baumgartner, G. Vucurevic, P. Stoeter and R. Treede, "Spatial resolution of fMRI in the human parasylvian cortex: Comparison of somatosensory and auditory activation," Neuroimage, vol. 25, pp. 877-887, 4/15. 2005. [5] C. E. Stark and Y. Okado, "Making memories without trying: medial temporal lobe activity associated with incidental memory formation during recognition," J. Neurosci., vol. 23, pp. 6748-6753, Jul 30. 2003. [6] A. Geissler, R. Lanzenberger, M. Barth, A. R. Tahamtan, D. Milakara, A. Gartus and R. Beisteiner, "Influence of fMRI smoothing procedures on replicability of fine scale motor localization," Neuroimage, vol. 24, pp. 323-331, Jan 15. 2005. 156 [7] B. Thirion, P. Pinel, A. Tucholka, A. Roche, P. Ciuciu, J. F. Mangin and J. B. Poline, "Structural Analysis of fMRI Data Revisited: Improving the Sensitivity and Reliability of fMRI Group Studies," IEEE Trans. Med. Imaging, vol. 26, 2007. [8] M. I. Miller, M. F. Beg, C. Ceritoglu and C. Stark, "Increasing the power of functional maps of the medial temporal lobe by using large deformation diffeomorphic metric mapping," Proc. Natl. Acad. Sci. U. S. A., vol. 102, pp. 9685-9690, Jul 5. 2005. [9] P. A. Yushkevich, J. A. Detre, D. Mechanic-Hamilton, M. A. Fernandez-Seara, K. Z. Tang, A. Hoang, M. Korczykowski, H. Zhang and J. C. Gee, "Hippocampus-specific fMRI group activation analysis using the continuous medial representation," Neuroimage, vol. 35, pp. 1516-1530, May 1. 2007. [10] R. T. Constable, P. Skudlarski, E. Mencl, K. R. Pugh, R. K. Fulbright, C. Lacadie, S. E. Shaywitz and B. A. Shaywitz, "Quantifying and Comparing Region-of-Interest Activation Patterns in Functional Brain MR Imaging: Methodology Considerations," Magn. Reson. Imaging, vol. 16, pp. 289-300, 4. 1998. [11] A. Nieto-Castanon, S. S. Ghosh, J. A. Tourville and F. H. Guenther, "Region of interest based analysis of functional imaging data," Neuroimage, vol. 19, pp. 1303-1316, 2003. [12] D. Kontos and V. Megalooikonomou, "Fast and effective characterization for classification and similarity searches of 2D and 3D spatial region data," Pattern Recognit, vol. 38, pp. 1831-1846, 11. 2005. 157 [13] B. Ng, R. Abugharbieh, X. Huang and M. J. McKeown, "Characterizing fMRI activations within regions of interest(ROIs) using 3D moment invariants," in MMBIA, 2006. [14] D. V. Vranic, "An improvement of rotation invariant 3D-shape based on functions on concentric spheres," in 2003, pp. 111-757-60 vol.2, 2003 [15] M. Kazhdan, T. Funkhouser and S. Rusinkiewicz, "Rotation invariant spherical harmonic representation of 3D shape descriptors," in Symposium on Geometry Processing, 23-25 June 2003, pp. 156-65, 2003. [16] S. Tootoonian, R. Abugharbieh, X. Huang and M. J. McKeown, "Shape vs. volume: Invariant shape descriptors for 3D region of interest characterization in MRI," in ISBI, 2006, pp. 754-757, 2006. [17] A. Uthama, R. Abhugharbieh, A. Traboulsee and M. J. McKeown, "Invariant SPHARM shape descriptors for complex geometry MR region of interest analysis in MR region of interest analysis," in 2007, pp. 1322-1325, 2007. [18] G. Burel and H. Henocq, "Three-dimensional invariants and their application to object recognition," Signal Process, vol. 45, pp. 1-22, 1995. [19] D. M. Healy, D. N. Rockmore, P. J. Kostelec and S. Moore, "FFTs for the 2-Sphere- Improvements and Variations," Journal of Fourier Analysis and Applications, vol. 9, pp. 341-385, 2003. 158 [20] Y. Zhang and A. M. Martinez, "Recognition of expression variant faces using weighted subspaces," in Proceedings of the 17th International Conference on Pattern Recognition, ICPR 2004, 2004, pp. 149-152. [21] Y. S. K. Vetsa, M. Styner, S. M. Pizer, J. A. Lieberman and G. E. -. Gerig, Caudate Shape Discrimination in Schizophrenia using Template-Free Non-Parametric Tests. , vol. 2879, 2003, pp. 661-669. [22] J. Ashburner and K. J. Friston, "Nonlinear spatial normalization using basis functions," Hum. Brain Mapp., vol. 7, pp. 254-266, 1999. [23] S. Hayasaka and T. E. Nichols, "Validating cluster size inference: random field and permutation methods," Neuroimage, vol. 20, pp. 2343-2356, 2003. [24] D. H. Brainard, "The Psychophysics Toolbox," Spat. Vis., vol. 10, pp. 433-436, 1997. [25] R. Liao, J. L. Krolik and M. J. McKeown, "An information-theoretic criterion for intrasubject alignment of FMRI time series: motion corrected independent component analysis," IEEE Trans. Med. Imaging, vol. 24, pp. 29-44, Jan. 2005. [26] J. Talairach and P. Tournoux, Co-Planar Stereotaxic Atlas of the Human Brain: 3- Dimensional Proportional System: An Approach to Cerebral Imaging. Thieme, 1988, [27] M. J. McKeown, L. K. Hansen and T. J. Sejnowski, "Independent component analysis of functional MRI: what is signal and what is noise," Curr. Opin. Neurobiol., vol. 13, pp. 620-629, 2003. 159 [28] J. F. Stein and M. Glickstein, "Role of the cerebellum in visual guidance of movement," Physiol. Rev., vol. 72, pp. 967-1017, Oct. 1992. [29] A. M. Owen, A. C. Evans and M. Petrides, "Evidence for a two-stage model of spatial working memory processing within the lateral frontal cortex: a positron emission tomography study," Cereb. Cortex, vol. 6, pp. 31-38, Jan-Feb. 1996. [30] S. Thobois, P. Dominey, J. Decety, P. Pollak, M. C. Gregoire and E. Broussolle, "Overactivation of primary motor cortex is asymmetrical in hemiparkinsonian patients," Neuroreport, vol. 11, pp. 785-789, Mar 20. 2000. [31] U. Sabatini, K. Boulanouar, N. Fabre, F. Martin, C. Carel, C. Colonnese, L. Bozzao, I. Berry, J. L. Montastruc, F. Chollet and O. Rascol, "Cortical motor reorganization in akinetic patients with Parkinson's disease: a functional MRI study," Brain, vol. 123 ( Pt 2), pp. 394-403, Feb. 2000. 160 7 Conclusions and Future Work 7.1 Technical Contributions This thesis presented a novel spherical harmonics based framework for region of interest shape and functional analysis of neurological MRI data. We first introduced an invariant spherical harmonic representation to characterize 3D real valued functions. Unlike previous methods [1, 2], which were invariant to rotation, scale and translation of the 3D function, the method proposed here ensure that the features are also unique. We achieved this additional uniqueness property by proposing the use of a radial transform on a shell based representation of 3D real valued functions. This ensures a one to one correspondence between the 3D real valued function and the SPHARM representation, an important requirement for any feature based representation. We then proposed the use of these features to discriminate shape changes in anatomical regions of interest in brain magnetic resonance imaging data. Unlike previous approaches which studied the 3D surface based representation [6, 7], we used these features to study shape as represented by the 3D volume representation of the regions. The surface based representation methods have a serious limitation of not being able to characterize regions with concavities [7] or regions without spherical topology [6]. Our features do not place any restrictions on the allowable topology, vastly increasing their potential. Synthetic data experiments, closely modeled after real world observations were used to validate the sensitivity of these features to subtle shape changes. 161 We then extended the use of these invariant features to study the spatial distribution of activation values within a 3D region of interest in functional magnetic resonance imaging data. This area of research is dominated by spatial normalization based voxel-by-voxel analysis approach [9] which has been criticized for its various shortcomings [10, 11, 12]. Feature based analysis of these regions has been proposed earlier [13] but an important aspect, namely, accounting for intersubject variability in the shape of the ROI has not been addressed. To enable the proposed SPHARM approach to sensitively detect changes in the functional pattern of activation in 3D space, while being relatively robust to structural changes in the shape of the ROI, we proposed a novel principal component subspace approach. Using synthetic experiments modeled after real world observations, we were able to show the superior sensitivity of this approach when compared to a recently proposed alternate feature based approach [13] and also to the current state-of-art spatial normalization approach in analyzing functional MR data. In summary, the contributions of this thesis are: 1. Translation, rotation and scale invariant properties of spherical harmonics were augmented to include the important property of uniqueness, a limitation of previous approaches. • Uthama A, Abugharbieh R, Traboulsee A, McKeown M.J, "Invariant SPHARM Shape Descriptors for Complex Geometry in MR Region of Interest Analysis. In 29th Annual International Conference IEEE EMBS: 2007; Lyon, France; 2007. 162 2. These features were used to characterize 3D volumes of region of interests (as opposed to 3D surfaces) in the brain exposing interesting, clinically relevant changes in the PD brain. Presented in: • Martin J. McKeown, Ashish Uthama, Rafeef Abugharbieh, Samantha Palmer, Mechelle Lewis and Xuemei Huang, "Shape (But Not Volume) Changes in the Thalami in Parkinson Disease", Submitted to BioMed Central Neuroscience, 2007. • Uthama A, Abugharbieh R, Traboulsee A, McKeown M. J, "Invariant SPHARM Shape Descriptors for Complex Geometry in MR Region of Interest Analysis" , in preparation. 3. Spherical harmonic features were used to characterize spatial distribution properties of entire 3D volumetric fMRI data in region of interests. First presented in: • Uthama A, Abugharbieh R, Traboulsee A, McKeown MJ, "3D Regional fMRI Activation Analysis Using SPHARM-Based Invariant Spatial Features and Intersubject Anatomical Variability Minimization", in preparation 4. A principal component subspace approach to reduce influence of anatomical variations in feature based representation of fMRI data was proposed. Presented along with clinical results in: 163 • Uthama A, Abugharbieh R, Palmer S.J, Traboulsee A and McKeown M. J, "Characterizing fMRI Activations in ROIs while Minimizing Effects of Intersubject Anatomical Variability", in preparation. 7.2 Impact on Neuroscience Research The spherical harmonics features developed were first used to study the shape of 3D volumes of specific regions of the brain. Structural MR image were obtained for both normal subjects and subjects afflicted with Parkinson's disease. The use of spherical harmonic features revealed shape changes in the thalamus and the brain ventricles in subjects with Parkinson's disease. These changes correlate well with previous pathological (post-mortem) observations [3]. This demonstrates that our proposed approach is a promising and powerful alternative means for neuroscientists to non- invasively observe such changes in vivo. The SPHARM framework developed for functional analysis was shown to have higher sensitivity than the leading (normalization-based) competing approach; it was able to detect clinically relevant changes in the functional response of PD patients. The improved sensitivity of the spherical harmonics achieved using principal component subspace projection enabled the detection of what is hypothesized to be a compensatory mechanism [4, 5]. These changes would not have been detected using the current state-of- art approach, thus underscoring the impact of the improved sensitivity demonstrated by our approach on the clinical inferences drawn from functional MR data. 164 7.3 Future Work The SPHARM based MR analysis framework of ROI data presented here has potential for further investigation. Some key directions include: • Given that the SPHARM features are valid for any 3D real valued function; this framework can be easily extended to analyze other MR data. Some key modalities where this could be beneficial include characterizing and discriminating scalar measures extracted from diffusion tensor (DTI) imaging data such as fractional anisotropy (FA). • Further study into the core of the SPHARM representation might help localize the global changes now observed, to a more local description. This will help make better clinical inferences by being able to isolate parts of the ROI responsible for the observed change. • Since the proposed features can handle disjoint ROIs, further investigations can be carried out by analyzing sets of ROIs simultaneously. This could potentially shed more light on the functional connectivity between different regions. The main advantage of the proposed framework is its generality due to the ease with which this method can be used to perform group discrimination analysis on any real valued 3D data. We expect this method to reveal more interesting clinical observations from MR data, which would have a significant impact in neuroscience. 165 7.4 References [1] M. Kazhdan, T. Funkhouser and S. Rusinkiewicz, "Rotation invariant spherical harmonic representation of 3D shape descriptors," in Symposium on Geometry Processing, 23 -25 June 2003, pp. 156 -65, 2003. [2] D. V. Vranic, "An improvement of rotation invariant 3D-shape based on functions on concentric spheres," in 2003, pp. III-757-60 vol.2, 2003. [3] J. M. Henderson, K. Carpenter, H. Cartwright and G. M. Halliday, "Loss of thalamic intralaminar nuclei in progressive supranuclear palsy and Parkinson's disease: clinical and therapeutic implications," Brain, vol. 123 ( Pt 7), pp. 1410-1421, Jul. 2000. [4] S. Thobois, P. Dominey, J. Decety, P. Pollak, M. C. Gregoire and E. Broussolle, "Overactivation of primary motor cortex is asymmetrical in hemiparkinsonian patients," Neuroreport, vol. 11, pp. 785-789, Mar 20. 2000. [5] U. Sabatini, K. Boulanouar, N. Fabre, F. Martin, C. Carel, C. Colonnese, L. Bozzao, I. Berry, J. L. Montastruc, F. Chollet and O. Rascol, "Cortical motor reorganization in akinetic patients with Parkinson's disease: a functional MRI study," Brain, vol. 123 ( Pt 2), pp. 394-403, Feb. 2000. [6] D. Goldberg-Zimring, A. Achiron, S. K. Warfield, C. R. G. Guttmann and H. Azhari, "Application of spherical harmonics derived space rotation invariant indices to the analysis of multiple sclerosis lesions' geometry by MRI," Magnetic Resonance Imaging, vol. 22, pp. 815-825, 2004. 166 [7] S. Tootoonian, R. Abugharbieh, X. Huang and M. J. McKeown, "Shape vs. volume: Invariant shape descriptors for 3D region of interest characterization in MRI," in ISBI, pp. 754-757, 2006. [8] Uthama A, Abugharbieh R, Traboulsee A, McKeown MJ: Invariant SPHARM Shape Descriptors for Complex Geometry in MR Region of Interest Analysis. In 29th Annual International Conference IEEE EMBS: 2007; Lyon, France; 2007. [9] K. J. Friston, A. P. Holmes, K. J. Worsley, J. B. Poline, C. Frith and R. S. J. Frackowiak, "Statistical Parametric Maps in Functional Imaging: A General Linear Approach," Human Brain Mapping, vol. 2, pp. 189-210, 1995. [10] F. Crivello, T. Schormann, N. Tzourio-Mazoyer, P. E. Roland, K. Zilles and B. M. Mazoyer, "Comparison of spatial normalization procedures and their impact on functional maps," Hum. Brain Mapp. (USA), vol. 16, pp. 228, 2002. [11] M. Ozcan, U. Baumgartner, G. Vucurevic, P. Stoeter and R. Treede, "Spatial resolution of fMRI in the human parasylvian cortex: Comparison of somatosensory and auditory activation," Neuroimage, vol. 25, pp. 877-887, 4/15. 2005. [12] F. L. Bookstein, ""Voxel-based morphometry" should not be used with imperfectly registered images," Neuroimage, vol. 14, pp. 1454-1462, 2001. [13] B. Ng, R. Abugharbieh, X. Huang and M. J. McKeown, "Characterizing fMRI activations within regions of interest (ROIs) using 3D moment invariants," in Computer vision and Pattern Recognition Workshop, pp 63, 2006, 167 [14] Martin J. McKeown, Ashish Uthama, Rafeef Abugharbieh, Samantha Palmer, Mechelle Lewis and Xuemei Huang, "Shape (But Not Volume) Changes in the Thalami in Parkinson Disease", Submitted to BioMed Central Neuroscience, 2007 [15] Uthama A, Abugharbieh R, Traboulsee A, McKeown M. J, "Invariant SPHARM Shape Descriptors for Complex Geometry in MR Region of Interest Analysis" , in preparation. [16] Uthama A, Abugharbieh R, Traboulsee A, McKeown MJ, "3D Regional fMRI Activation Analysis Using SPHARM-Based Invariant Spatial Features and Intersubject Anatomical Variability Minimization", in preparation [17] Uthama A, Abugharbieh R, Palmer S.J, Traboulsee A and McKeown M. J, "Characterizing fMRI Activations in ROIs while Minimizing Effects of Intersubject Anatomical Variability", in preparation. 168 8 APPENDIX All data obtained and analyzed during the course of this thesis was obtained from human subjects after due consent was obtained. The studies were approved by the concerned ethics review boards. This appendix lists the ethics approval certificates obtained from the University of British Columbia Clinical Research Ethics Board. 169 ]taps r1,4 obc e,t,r1., 1^0 SI ADN1331M81' 11RPNVKRCirYUl39 fromStrriv 11114 The University of British Columbia Office of Research Services Clinical Research Ethics Board — Room 210, 828 West 10th Avenue, Vancouver, BC V5Z 1L8 ETHICS CERTIFICATE OF EXPEDITED APPROVAL: RENEWAL WITH AMENDMENTS TO THE STUDY PRINCIPAL INVESTIGATOR: Martin J. McKeown DEPARTMENT: UBC/Medicine, Faculty of/ Medicine, Department of/ Neurology - Med UBC CREB NUMBER: H04-70177 INSTITUTION(S) WHERE RESEARCH WILL BE CARRIED OUT: Institution^1^Site^ I UBC^ Point Grey Site Other locations where the research will be conducted: nia CO-INVESTIGATOR(S): Joseph K.C. Tsui Chong S, Lee Donald B. Caine A Jon Stoessl SPONSORING AGENCIES: Canadian Institutes of Health Research (CIHR) - ''Making the connection: Methods to Infer Functional Connectivity in brain studies" Natural Sciences and Engineering Research Council of Canada (NSERC) - "Making the connection: Methods to Infer Functional Connectivity in brain studies" UBC Internal Grant - "Motor & Cognitive Functional Magnetic Resonance Imaging Studies in Parkinson's & Related Disorders" PROJECT TITLE: Making the connection: Methods to Infer Functional Connectivity in brain studies The current UBC CREB approval for this study expires: June 1, 2008 IAMENDMENT(S):  AMENDMENT APPROVAL 3ATE: idt^icbc c^0,SIADNI334N1814181'NNKRulTUV34,fronttintitg land (I of 2161 ,2007 &41.39 AM 1 70

Cite

Citation Scheme:

    

Usage Statistics

Country Views Downloads
China 24 13
United States 5 2
United Kingdom 3 0
Iran 1 1
City Views Downloads
Beijing 24 2
Unknown 4 18
Mountain View 3 1
Ashburn 2 0

{[{ mDataHeader[type] }]} {[{ month[type] }]} {[{ tData[type] }]}

Share

Share to:

Comment

Related Items