# Segmentation of Striatal Brain Structures from High Resolution PET Images

Authors: Ricardo J. P. C. Farinha; Ulla Ruotsalainen; Jussi Hirvonen; Lauri Tuominen; Jarmo Hietala; José M. Fonseca; Jussi Tohka
Journal: International Journal of Biomedical Imaging
 (2009)

---

## Abstract

We propose and evaluate an automatic segmentation method for extracting striatal brain structures (caudate, putamen, and ventral striatum) from parametricC11-raclopride positron emission tomography (PET) brain images. We focus on the images acquired using a novel brain dedicated high-resolution (HRRT) PET scanner. The segmentation method first extracts the striatum using a deformable surface model and then divides the striatum into its substructures based on a graph partitioning algorithm. The weighted kernel k-means algorithm is used to partition the graph describing the voxel affinities within the striatum into the desired number of clusters. The method was experimentally validated with synthetic and real image data. The experiments showed that our method was able to automatically extract caudate, ventral striatum, and putamen from the images. Moreover, the putamen could be subdivided into anterior and posterior parts. An automatic method for the extraction of striatal structures from high-resolution PET images allows for inexpensive and reproducible extraction of the quantitative information from these images necessary in brain research and drug development.


---

## Body

## 1. Introduction

POSITRON emission tomography (PET) is a widely used functional imaging technology. PET imaging allows measuring physiological processes in human’s in vivo and computing three dimensional (3D) images quantifying the parameters of these processes. This functional information enhances the understanding of biochemical basis of normal and abnormal bodily functions [1]. The type of the functional information provided by the PET images depends on the applied radiopharmaceutical that is injected into the subject. To create quantitative 3D images of the physiological parameter of interest—so called parametric images—one needs to study the dynamics of the concentration of the applied radiopharmaceutical at each voxel and relate these to the input function. The input function is the concentration of the radiopharmaceutical in either arterial plasma or in a specific brain region known as the reference region [2]. In this study, raclopride labeled with 11C is used as the radiopharmaceutical. 11C-Raclopride is a selective antagonist and binds reversibly to dopamine D2 receptors. The parameter of interest in this case is the binding potential (BP) of D2 dopamine receptors which reflects the receptor density. The striatum contains the highest density of D2 receptors in the brain. In this work, we propose an automatic method to segment striatal substructures in 11C-raclopride PET images. Here, we consider BP as the ratio of specific to free plus nonspecific binding in tissue, and refer to it this point forward as BPND [2]. An example of a parametric image representing the BPND of 11C-raclopride is shown in Figure 1.An example of one11C-raclopride PET image of a human head: (a) transaxial, (b) sagittal, (c) coronal views.
(a)(b)(c)Most image segmentation in clinical environments is currently done by manual slice editing, meaning that a human operator manually delineates the region of interest (ROI) on each2D slice of the 3D image. This process has several disadvantages. For example, the results will be dependent on any subjective bias on part of the operator deducing the spatial shape of 3D brain structures by viewing the 2D slices. Thus, the final result might not be reproducible by another operator, hampering the comparability between individual raters. Second, the enormous amount of data that must be manually or semiautomatically extracted makes the manual segmentation process time-consuming and expensive.We have earlier presented an automatic method to extract striatal structures, namely left and right caudate and putamen, from11C-raclopride PET images [3]. This method has been shown to produce more reproducible segmentations than the traditional manual slice editing [3, 4]. So far, this method has not been applied to the images acquired using the newest generation of PET scanners, the High-Resolution Research Tomograph (HRRT), because a new, improved method is needed for the segmentation of these images.There are two primary reasons for this need. First, the segmentation method presented in [3] relies on the computation of the eigenvalues and eigenvectors of an image derived affinity matrix. However, when the image size increases, eigenvalue computation becomes computationally prohibitive especially regarding the memory requirements. The images acquired using a high-resolution PET scanner (ECAT HRRT, Siemens Medical Solutions, Knoxville, Tenn, USA) are of a size that exceeds the limit of making the eigenvalue computations impractical. The reconstructed dynamic images have approximately 800 Mbytes of data and the size of the affinity matrix is approximately 16000 × 16000. Second, the method used in [3] was designed to extract only caudate and putamen, but, with the resolution enhancement offered by HRRT, more substructures can be separated.Accurate segmentation of the striatum from high-resolution PET images is of great interest because the functional anatomy of this structure is far more complex than what can be captured by previous generation PET scanners with lower spatial resolution. Indeed, the striatum can be divided based on topographically organized anatomical connections with cerebellar cortex into limbic (nucleus accumbens, and ventral parts of caudate and precommissural putamen), associative (most of the head of caudate and ventral parts of precommissural putamen), and sensorimotor (postcommissural putamen) divisions [5, 6]. This anatomical framework may have implications for the regional distribution of D2 receptors within the striatum. Therefore, in the present work, we introduce an automatic technique that is able to extract multiple striatal structures from parametric 11C-raclopride BPND images generated with HRRT PET.The segmentation problem is divided in three parts, the extraction of the left and right striatum, the creation of a striatal image graph, and partitioning the striatal image graph into a predefined number of segments corresponding to striatal substructures. The extraction of the striatum is performed by a global optimization based deformable surface model. Then, a weighted striatal image graph, whose nodes correspond to striatal voxels, is constructed.  The weights of the edges in the graph are computed based on the spatial locations and intensity values of voxels. The image graph can be represented by a matrix of the edge weights called the affinity matrix. Once the affinity matrix is constructed, we apply a graph partitioning algorithm to find the optimal partitioning. Spectral methods, such as normalized cuts algorithm applied in [3], are a widely used choice for partitioning. However, spectral clustering methods rely on the computation of the eigenvalues and eigenvectors of the affinity matrix to partition the graph. In cases where the large size of affinity matrix makes the eigenvector computation prohibitive, spectral methods will fail. The HRRT data used in this study produced large affinity matrix, rendering the spectral methods impractical. Therefore, our implementation of a graph partitioning method has to dwell with this data size problem. In [7] it was shown that the weighted kernel k-means clustering criterion with particular selections of kernel and weights is equivalent to the multiway normalized cuts criterion. The weighted kernel k-means criterion is a generalization of the well known k-means criterion. It can be locally optimized using a simple and fast algorithm similar to k-means clustering algorithm. This equivalence provides a principled way to avoid the eigenvalue computations to solve normalized cuts type graph partitioning criteria.


## 2. Methods

Our method for the segmentation of striatal structures from HRRT PET images proceeds in a top-down manner.  The automatic method consists of a series of consecutive steps shown in the hierarchy graph in Figure2.Figure 2
Steps of the segmentation of caudate, putamen, and ventral striatum: (a) preprocessing steps, (b) affinity matrix construction method, (c) clustering process.The preprocessing step consists of the extraction of the striatum, from HRRT PETBPND images, see Figure 2(a). This extraction method is the same as the one used in the earlier method [3] except for slight tuning of the initial parameters.Once the striatum has been extracted an affinity matrix containing the similarity values between all pairs of striatal voxels is constructed using voxel coordinates and intensities as features, see Figure2(b).  We call this step feature extraction and it is described in Section 2.2.Finally, the clustering algorithm segments the striatum into components corresponding to caudate, ventral striatum and anterior and posterior putamen, see Figure2(c). The clustering process is treated as a graph partitioning problem, where the affinity matrix defines to the graph to be partitioned. The clustering algorithm is described in Section 2.3.


### 2.1. Preprocessing

Before we can proceed with the striatum segmentation, we need to segment left and right striatum from the background. This preprocessing can be subdivided in three steps: extraction of the brain surface using a deformable surface model [8], segmentation of the brain image in right and left brain hemispheres based on the extracted brain surface [8], and finally extraction of the striatum from each hemisphere [3]. The previous method [3] was found to be applicable for also HRRT-PET images with a minor parameter tuning and we refer to [3] for more detailed information on these methods.


### 2.2. Feature Extractor

The feature extractor aims to construct an affinity matrix using features extracted from the analyzed image.  As we have already mentioned, the affinity matrix can be understood as a graph and vice versa. This graph contains the similarity values between all pairs of voxels belonging to the striatum. The image features used in our method are coordinates of voxel centers and voxel intensity values, that is,BPND values. Let D denote the set of voxels within left or right striatum and let I(i) be the intensity value of the voxel i=(ix,iy,iz) belonging to D.  We define the similarity A(i,j)′ between voxels i and j as
(1)A(i,j)′={exp(-2|I(i)-I(j)|maxk∈DI(k))ifiis connected toj0otherwise.The voxeli is said to be connected to the voxel j if j is included in the neighborhood of i. Voxels i and j are neighbors if the Euclidean distance between their centers is less than given constant. This constant should be selected to reflect the noise level and spatial resolution of the images to be segmented: with noisy images the constant needs to be lower than with less noisy images.  Assuming that I(i)>0 for all i, A(i,j) takes values between exp(-2) and 1 if the voxels i and j are connected. We still enhance the boundaries between the different structures by thresholding the affinity matrix. The final thresholded A(i,j) affinity matrix is defined as
(2)A(i,j)={A(i,j)′A(i,j)′>PT0otherwise.In (2), A(i,j) corresponds to the final affinity value between voxels i and j, while A(i,j)′ is computed by (1).  This second thresholding removes connections between any two voxels with   A(i,j)′ inferior to a chosen threshold. This threshold is defined in (2) as Percentage Threshold (PT). The values of PT can vary from 0 to 1. In this work, we set PT=0.5 based on experimental data.


### 2.3. Clustering Process

We partition the graph represented by the affinity matrixA into subgraphs corresponding the caudate, anterior and posterior putamen, ventral striatum, and background using weighted kernel k-means algorithm.As already mentioned, the weighted kernelk-means algorithm is a generalization of the well known k-means algorithm. The weighted kernel k-means criterion function is [7]
(3)d({πi}i=1k)=∑i=1k∑a∈πiW(a)∥ϕ(a)-mi∥2.In (3), πi denotes ith cluster with the cluster centre mi=∑b∈πiW(b)ϕ(b)/∑b∈πiW(b) and W(b) is the weight assigned to the voxel b. These weight values are computed based on the affinity matrix (2) as we will show in what follows. The nonlinear function ϕ maps the original (input) space to a higher-dimensional feature space. This mapping allows extracting the clusters that are not linearly separable in the input space. This mapping is also computed based on the affinity matrix.LetW be a diagonal matrix of the weights W(a) for all voxels in D and Wj be a diagonal matrix of weights in the cluster πj. Further, collect all the inner products ϕ(a)·ϕ(b),a,b∈D into the kernel matrix K. The weighted kernel k-means criterion can be locally optimized by Algorithm 1Algorithm 1:Graph partition algorithm (Weighted Kernel k-means). This algorithm assigns the voxel ai to the cluster πc if the cluster πc minimizes the distance function from point 2 [7].
Function[{πc}c=1k]=KernelKmeans(K,k,W,tmax,{πc(0)}c=1k)Input→K:Kernel matrix;k:number of clusters;W:sum of weights for each point;tmax=maximumnumberofiterations;{πc(0)}c=1k:intial cluster.Output→{πc}c=1k:final clusters.1.  If no intial clustering is given, initialize thek clusters (i.e., randomly).Set the iteration countert=0.2. For each pointai and every cluster c, compute:d(ai,mc)=Kii-2∑aj∈πc(t)Wjkij∑aj∈πc(t)Wj+∑aj,al∈πc(t)WjWlKjl(∑aj∈πc(t)Wj)23.Findc*(ai)=argmincd(ai,mc),resolvingtiesarbitrarily.Compute the updated clusterπc(t+1)={a:c*(ai)=c}.4. If not converged ortmax>t,set t=t+1 and go to step 2;Otherwise, stop and output the final clusters{πc}c=1k={πc(t+1)}c=1kThe weighted form of the kernelk-means is mathematically equivalent to a general weighted graph clustering [7]. Such equivalence implies that the weighted kernel k-means algorithm can be used to locally optimize a number of graph clustering objectives such as normalized cuts. The equivalence of (3) to the normalized cuts criterion is obtained by selecting
(4)W(i,i)=∑jA(i,j),K=W-1AW-1,
where A is the affinity matrix defined in (2).Although the weighted kernelk-means algorithm allows for segmenting a large number of striatal voxels into the desired substructures of interest, it has some subtleties. Namely, the final result depends on the desired number of clusters and the initial clustering. These parameters have to be carefully tuned to provide the most accurate result. We examined several different choices for these parameters and found an initial configuration that was both reliable and accurate in this application. Particularly, the initial clustering was formed by dividing the striatal domain D into k=4 components of the equal size based on a lexicographical ordering. Voxels in the striatum were lexicographically ordered first scanning the domain from left to right direction, then from anterior to posterior direction, and finally from inferior to superior direction. Technically, for the image with dimensions X×Y×Z, each voxel i=(ix,iy,iz)  in D was assigned a score s(i)=ix+Xiy+XYiz, where ix,iy, and iz are coordinates of i in the left-right, anterior-posterior, and inferior-superior axes, respectively. The voxels were ordered based on these scores: the voxel i was initially assigned to the cluster πc if   (c-1)|D|/4<r(i)≤c|D|/4, where  r(i) is  the rank of s(i) among the voxels in D. No differences to the final results were observed if instead of scanning from left to right we would scan from lateral to medial in the case of the right striatum. Note that the initialization process is deterministic, that is, there is no random component in our algorithm. This is important for the current application because it ensures the reproducibility of the automatic segmentations.


## 2.1. Preprocessing

Before we can proceed with the striatum segmentation, we need to segment left and right striatum from the background. This preprocessing can be subdivided in three steps: extraction of the brain surface using a deformable surface model [8], segmentation of the brain image in right and left brain hemispheres based on the extracted brain surface [8], and finally extraction of the striatum from each hemisphere [3]. The previous method [3] was found to be applicable for also HRRT-PET images with a minor parameter tuning and we refer to [3] for more detailed information on these methods.


## 2.2. Feature Extractor

The feature extractor aims to construct an affinity matrix using features extracted from the analyzed image.  As we have already mentioned, the affinity matrix can be understood as a graph and vice versa. This graph contains the similarity values between all pairs of voxels belonging to the striatum. The image features used in our method are coordinates of voxel centers and voxel intensity values, that is,BPND values. Let D denote the set of voxels within left or right striatum and let I(i) be the intensity value of the voxel i=(ix,iy,iz) belonging to D.  We define the similarity A(i,j)′ between voxels i and j as
(1)A(i,j)′={exp(-2|I(i)-I(j)|maxk∈DI(k))ifiis connected toj0otherwise.The voxeli is said to be connected to the voxel j if j is included in the neighborhood of i. Voxels i and j are neighbors if the Euclidean distance between their centers is less than given constant. This constant should be selected to reflect the noise level and spatial resolution of the images to be segmented: with noisy images the constant needs to be lower than with less noisy images.  Assuming that I(i)>0 for all i, A(i,j) takes values between exp(-2) and 1 if the voxels i and j are connected. We still enhance the boundaries between the different structures by thresholding the affinity matrix. The final thresholded A(i,j) affinity matrix is defined as
(2)A(i,j)={A(i,j)′A(i,j)′>PT0otherwise.In (2), A(i,j) corresponds to the final affinity value between voxels i and j, while A(i,j)′ is computed by (1).  This second thresholding removes connections between any two voxels with   A(i,j)′ inferior to a chosen threshold. This threshold is defined in (2) as Percentage Threshold (PT). The values of PT can vary from 0 to 1. In this work, we set PT=0.5 based on experimental data.


## 2.3. Clustering Process

We partition the graph represented by the affinity matrixA into subgraphs corresponding the caudate, anterior and posterior putamen, ventral striatum, and background using weighted kernel k-means algorithm.As already mentioned, the weighted kernelk-means algorithm is a generalization of the well known k-means algorithm. The weighted kernel k-means criterion function is [7]
(3)d({πi}i=1k)=∑i=1k∑a∈πiW(a)∥ϕ(a)-mi∥2.In (3), πi denotes ith cluster with the cluster centre mi=∑b∈πiW(b)ϕ(b)/∑b∈πiW(b) and W(b) is the weight assigned to the voxel b. These weight values are computed based on the affinity matrix (2) as we will show in what follows. The nonlinear function ϕ maps the original (input) space to a higher-dimensional feature space. This mapping allows extracting the clusters that are not linearly separable in the input space. This mapping is also computed based on the affinity matrix.LetW be a diagonal matrix of the weights W(a) for all voxels in D and Wj be a diagonal matrix of weights in the cluster πj. Further, collect all the inner products ϕ(a)·ϕ(b),a,b∈D into the kernel matrix K. The weighted kernel k-means criterion can be locally optimized by Algorithm 1Algorithm 1:Graph partition algorithm (Weighted Kernel k-means). This algorithm assigns the voxel ai to the cluster πc if the cluster πc minimizes the distance function from point 2 [7].
Function[{πc}c=1k]=KernelKmeans(K,k,W,tmax,{πc(0)}c=1k)Input→K:Kernel matrix;k:number of clusters;W:sum of weights for each point;tmax=maximumnumberofiterations;{πc(0)}c=1k:intial cluster.Output→{πc}c=1k:final clusters.1.  If no intial clustering is given, initialize thek clusters (i.e., randomly).Set the iteration countert=0.2. For each pointai and every cluster c, compute:d(ai,mc)=Kii-2∑aj∈πc(t)Wjkij∑aj∈πc(t)Wj+∑aj,al∈πc(t)WjWlKjl(∑aj∈πc(t)Wj)23.Findc*(ai)=argmincd(ai,mc),resolvingtiesarbitrarily.Compute the updated clusterπc(t+1)={a:c*(ai)=c}.4. If not converged ortmax>t,set t=t+1 and go to step 2;Otherwise, stop and output the final clusters{πc}c=1k={πc(t+1)}c=1kThe weighted form of the kernelk-means is mathematically equivalent to a general weighted graph clustering [7]. Such equivalence implies that the weighted kernel k-means algorithm can be used to locally optimize a number of graph clustering objectives such as normalized cuts. The equivalence of (3) to the normalized cuts criterion is obtained by selecting
(4)W(i,i)=∑jA(i,j),K=W-1AW-1,
where A is the affinity matrix defined in (2).Although the weighted kernelk-means algorithm allows for segmenting a large number of striatal voxels into the desired substructures of interest, it has some subtleties. Namely, the final result depends on the desired number of clusters and the initial clustering. These parameters have to be carefully tuned to provide the most accurate result. We examined several different choices for these parameters and found an initial configuration that was both reliable and accurate in this application. Particularly, the initial clustering was formed by dividing the striatal domain D into k=4 components of the equal size based on a lexicographical ordering. Voxels in the striatum were lexicographically ordered first scanning the domain from left to right direction, then from anterior to posterior direction, and finally from inferior to superior direction. Technically, for the image with dimensions X×Y×Z, each voxel i=(ix,iy,iz)  in D was assigned a score s(i)=ix+Xiy+XYiz, where ix,iy, and iz are coordinates of i in the left-right, anterior-posterior, and inferior-superior axes, respectively. The voxels were ordered based on these scores: the voxel i was initially assigned to the cluster πc if   (c-1)|D|/4<r(i)≤c|D|/4, where  r(i) is  the rank of s(i) among the voxels in D. No differences to the final results were observed if instead of scanning from left to right we would scan from lateral to medial in the case of the right striatum. Note that the initialization process is deterministic, that is, there is no random component in our algorithm. This is important for the current application because it ensures the reproducibility of the automatic segmentations.


## 3. Experiments

### 3.1. Phantom Study

We performed experiments with simulated HRRT PET phantoms to evaluate the automatic segmentation method. For this, we generated dynamic phantoms corresponding11C-raclopride brain images. The construction of the phantoms required a model of anatomy (a labeled 3D brain volume) and a model of physiology (Time Activity Curves (TACs) for each brain region).The anatomical model was derived from a manually segmented T1-weighted brain MR image (the Jacob phantom by Montreal Neurological Institute [9]). We reduced the number of brain structures represented in the Jacob phantom to 9 distinct anatomical structures for our anatomical model. The structures in the anatomical model were anterior (dorsal) and posterior putamen, ventral striatum (including only nucleus accumbens in our model), caudate (dorsal), white and gray matter, cerebellum, sinuses, scalp, and skull. The putaminal subdivision was defined according to previously published criteria [10].The physiological model consisted of a set of TACs, each one describing the dynamic behavior of the radiopharmaceutical in a different anatomical structure. The TAC’s for the individual structures were taken from [11]. In [11], the model of physiology did not include separate TACs for posterior/ anterior putamen or ventral striatum. Therefore, an estimation of those TAC’s had to be done in order to test our method properly. The TAC for the anterior putamen was obtained by multiplying the TAC of the posterior putamen by a constant (1.05). Similarly, the TAC for the ventral striatum was obtained by multiplying the caudate TAC by a constant (0.9). These are reasonable approximations based on our experience with these data. We additionally generated images using a set of different constant values. Different constants led to different BPND values in the structures, and for us, increased or decreased contrast between the structures.  However, slightly altered contrasts caused only minor changes to the segmentation results and therefore we only report quantitative results with phantoms constructed using the physiological model specified above.The noiseless dynamic PET phantom was generated by assigning each voxel and each time frame to a count value based on the anatomical model used and the time activity curves (TACs) of the radiopharmaceutical. The dimensions for the phantom were256×256×200 with 26 time frames. The voxel size was 1.22mm×1.22mm×1.22mm.  The phantom was convolved with a 2.5 mm3 Gaussian kernel to model the resolution of the ECAT HRRT scanner [12]. After generating radioactive decay effects, noise was added by drawing voxelwise counts from the Poisson distribution with the mean as the count value in the noiseless convolved phantom. This ideal level of noise is much lower than the one encountered in reality due to the effects of scatter, attenuation, randoms, and detector dead time. In addition to the phantom with ideal noise level, we generated a phantom with a more realistic noise level. For this we multiplied the voxel values by an estimated survival probability of 0.04 before generating the noise. The survival probability models the probability that a photon pair resulting from positron-electron annihilation is detected by the scanner. Our rough estimate of the survival probability was based on the approximation that 80% of the true counts are lost due to attenuation and other 80% are lost due to detector dead-time. Other factors affecting the noise level were ignored in this estimate. The performance evaluations of HRRT PET scanners are available [12], but linking these performance characteristics to image noise estimates is challenging [13], and thus we settled for this rough noise level estimate.  The images were decay corrected, multiplied by the inverse of the survival probability and the units were converted to nCi/ml. Finally, the BPND image was computed using the simplified reference region model using the cerebellum as the reference region [14]. Examples of the phantom with ideal and realistic noise levels can be seen in Figure 3.HRRT phantom. (a), (b), and (c) represent  the transaxial, sagittal, and coronal views of the HRRT phantom with the ideal noise level, while (d), (e), and (f) represent  to the same phantom with the more realistic noise level. Figures3(g) to 3(l) are the zoomed (400%) to the striata in the images presented in (a) to (f).
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)With the phantom images, it is possible to quantitatively evaluate our method because we know the ground truthBPND values for each structure. Because the main interest was in evaluating the accuracy of the extracted regional BPND values, we limit the discussion on the evaluation of the accuracy of structure-wise BPND values computed based on automatic segmentation. The structure-wise BPND values were computed by averaging the voxel-wise BPND values in a parametric image for a given structure (left and right caudate, ventral striatum, and anterior and posterior putamen). The performance of our method was measured using the normalized absolute similarity (NAS), which is defined:
(5)NAS=1-(2|BPa-BPb|(BPa+BPb)).In (5), BPa and BPb represent, respectively, the average BPND values computed based on the automatic segmentation and the ground truth BPND value computed based on the physiological ground truth TAC. The NAS values range from -1 to 1 assuming that BPND values are positive. This and closely related performance measures are widely applied to assess test-retest variability within PET imaging (e.g., [3, 4, 15]) which simplifies its interpretation relative to other studies.We also compared the automatically segmented volumes to the anatomical ground-truth volume. The applied figure of merit was positive predictive value (PPV) introduced in our earlier work [3] to estimate the quality of PET image segmentations. PPV was computed for every structure separately. We define true positives (TPs)as voxels which truly belong to a certain brain structure and are also labeled as belonging to that same brain structure in an automatic segmentation. False positives (FPs) are voxels which do not truly belong to a brain structure in the anatomical ground-truth but are labeled as belonging to it in an automatic segmentation. PPV is now defined as
(6)|TP||TP|+|FP|,
where |TP| (|FP|) is the number of TPs (FPs). PPV is used as the figure of merit here, because, as argued in [3], it is important that an extracted brain structure contains only voxels truly belonging to that brain structure and it is not that important whether or not we detect all the voxels belonging to that brain structure.


### 3.2. PET Studies with Healthy Volunteers

We compared our automatic segmentation method to manual ROI segmentation. This comparison was done using11C-raclopride BPND high-resolution PET images. The PET scans were provided by Turku PET Center and were acquired from four right-handed, nonsmoking healthy male volunteers. This study was approved by the Joint Ethical Committee of the University of Turku and Turku University Central Hospital, and was conducted according to the Declaration of Helsinki. All subjects gave ethical committee-approved written informed consents. All subjects were free of any somatic or psychiatric illness, illicit drug abuse, and alcoholism, as confirmed by blood and urine tests, somatic assessment, and a structured clinical interview for DSM-IV Axis I disorders conducted by a psychiatrist. The age, height, and weight of the subjects were 27.5±7.74; (20.3 to 33.1) years, 175±1 (174 to 176) cm, and 76.1±5.16 (70.0 to 82.5) kg, respectively (mean, s.d.; range). All subjects underwent T1-weighted magnetic resonance imaging at 1.5T to rule out structural brain abnormalities.The radioligand11C-raclopride was prepared as previously described [15]. PET imaging was performed by a brain-dedicated ECAT HRRT, Siemens Medical Solutions, Knoxville, Tnee, USA) [12].  In the acquired measurements, a molded thermoplastic mask was used as a head fixation to reduce the noise due to head movements. The right antecubital vein was cannulated and 11C-raclopride was administered as an intravenous rapid bolus flushed with saline. Injected dose and the specific radioactivity were 443.25±104.69 MBq and 185.7±218.4 MBq/nmol, respectively (mean ± S.D.). The PET 11C-raclopride scans were acquired in list mode and histogrammed by axial compression of span 9. The 11C-raclopride uptake was measured for fifty one minutes after injection and the frame sequence consisted of three frames of sixty seconds, four frames of one hundred and eighty seconds, and six frames of three hundred and sixty seconds. The image size of each scan was 200×200×150 slices and they had an isotropic voxel dimension of 1.22×1.22×1.22 mm3. The raw data was reconstructed with a speed-optimized version of OP-OSEM-3D (Ordinary Poisson-OSEM in full 3D) [16], with sixteen subsets and eight full iterations [17]. The frame-images of each dynamic PET image were first realigned to correct for head movement during the scan.  Integral image (summed over time) was calculated from the realigned PET-image and the MRI was coregistered to this PET-image. All realignment and coregistration procedures were performed using Statistical Parametric Mapping software version 2 (SPM2, Wellcome Department of Cognitive Neurology, University College London, UK). The parametric BPND images were created from the dynamic images using the simplified reference tissue model with the cerebellar cortex as the reference region [14].Regions of interest (ROIs) were manually drawn to coregistered transaxial MRI using Imadeus software (version 1.2., Forima Inc., Turku, Finland). ROIs onto dorsal caudate, dorsal putamen, ventral striatum, thalamus, and cerebellum were defined as previously described [16]. Simplified reference tissue model with the cerebellar cortex as the reference region [14] was applied to TAC data to obtain the regional BPND values.Average NAS (ANAS) over the four subjects was used to measure the similarity of theBPND values based on the automatic segmentation and the manual segmentation. In this case the BPa and BPb in (5) represent, respectively, the BPND values computed based on the automatic segmentation and manual segmentation.  In this case, however, it cannot be claimed that the BPND values based on the manual segmentation would be correct although they can be used as benchmarks.


## 3.1. Phantom Study

We performed experiments with simulated HRRT PET phantoms to evaluate the automatic segmentation method. For this, we generated dynamic phantoms corresponding11C-raclopride brain images. The construction of the phantoms required a model of anatomy (a labeled 3D brain volume) and a model of physiology (Time Activity Curves (TACs) for each brain region).The anatomical model was derived from a manually segmented T1-weighted brain MR image (the Jacob phantom by Montreal Neurological Institute [9]). We reduced the number of brain structures represented in the Jacob phantom to 9 distinct anatomical structures for our anatomical model. The structures in the anatomical model were anterior (dorsal) and posterior putamen, ventral striatum (including only nucleus accumbens in our model), caudate (dorsal), white and gray matter, cerebellum, sinuses, scalp, and skull. The putaminal subdivision was defined according to previously published criteria [10].The physiological model consisted of a set of TACs, each one describing the dynamic behavior of the radiopharmaceutical in a different anatomical structure. The TAC’s for the individual structures were taken from [11]. In [11], the model of physiology did not include separate TACs for posterior/ anterior putamen or ventral striatum. Therefore, an estimation of those TAC’s had to be done in order to test our method properly. The TAC for the anterior putamen was obtained by multiplying the TAC of the posterior putamen by a constant (1.05). Similarly, the TAC for the ventral striatum was obtained by multiplying the caudate TAC by a constant (0.9). These are reasonable approximations based on our experience with these data. We additionally generated images using a set of different constant values. Different constants led to different BPND values in the structures, and for us, increased or decreased contrast between the structures.  However, slightly altered contrasts caused only minor changes to the segmentation results and therefore we only report quantitative results with phantoms constructed using the physiological model specified above.The noiseless dynamic PET phantom was generated by assigning each voxel and each time frame to a count value based on the anatomical model used and the time activity curves (TACs) of the radiopharmaceutical. The dimensions for the phantom were256×256×200 with 26 time frames. The voxel size was 1.22mm×1.22mm×1.22mm.  The phantom was convolved with a 2.5 mm3 Gaussian kernel to model the resolution of the ECAT HRRT scanner [12]. After generating radioactive decay effects, noise was added by drawing voxelwise counts from the Poisson distribution with the mean as the count value in the noiseless convolved phantom. This ideal level of noise is much lower than the one encountered in reality due to the effects of scatter, attenuation, randoms, and detector dead time. In addition to the phantom with ideal noise level, we generated a phantom with a more realistic noise level. For this we multiplied the voxel values by an estimated survival probability of 0.04 before generating the noise. The survival probability models the probability that a photon pair resulting from positron-electron annihilation is detected by the scanner. Our rough estimate of the survival probability was based on the approximation that 80% of the true counts are lost due to attenuation and other 80% are lost due to detector dead-time. Other factors affecting the noise level were ignored in this estimate. The performance evaluations of HRRT PET scanners are available [12], but linking these performance characteristics to image noise estimates is challenging [13], and thus we settled for this rough noise level estimate.  The images were decay corrected, multiplied by the inverse of the survival probability and the units were converted to nCi/ml. Finally, the BPND image was computed using the simplified reference region model using the cerebellum as the reference region [14]. Examples of the phantom with ideal and realistic noise levels can be seen in Figure 3.HRRT phantom. (a), (b), and (c) represent  the transaxial, sagittal, and coronal views of the HRRT phantom with the ideal noise level, while (d), (e), and (f) represent  to the same phantom with the more realistic noise level. Figures3(g) to 3(l) are the zoomed (400%) to the striata in the images presented in (a) to (f).
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)With the phantom images, it is possible to quantitatively evaluate our method because we know the ground truthBPND values for each structure. Because the main interest was in evaluating the accuracy of the extracted regional BPND values, we limit the discussion on the evaluation of the accuracy of structure-wise BPND values computed based on automatic segmentation. The structure-wise BPND values were computed by averaging the voxel-wise BPND values in a parametric image for a given structure (left and right caudate, ventral striatum, and anterior and posterior putamen). The performance of our method was measured using the normalized absolute similarity (NAS), which is defined:
(5)NAS=1-(2|BPa-BPb|(BPa+BPb)).In (5), BPa and BPb represent, respectively, the average BPND values computed based on the automatic segmentation and the ground truth BPND value computed based on the physiological ground truth TAC. The NAS values range from -1 to 1 assuming that BPND values are positive. This and closely related performance measures are widely applied to assess test-retest variability within PET imaging (e.g., [3, 4, 15]) which simplifies its interpretation relative to other studies.We also compared the automatically segmented volumes to the anatomical ground-truth volume. The applied figure of merit was positive predictive value (PPV) introduced in our earlier work [3] to estimate the quality of PET image segmentations. PPV was computed for every structure separately. We define true positives (TPs)as voxels which truly belong to a certain brain structure and are also labeled as belonging to that same brain structure in an automatic segmentation. False positives (FPs) are voxels which do not truly belong to a brain structure in the anatomical ground-truth but are labeled as belonging to it in an automatic segmentation. PPV is now defined as
(6)|TP||TP|+|FP|,
where |TP| (|FP|) is the number of TPs (FPs). PPV is used as the figure of merit here, because, as argued in [3], it is important that an extracted brain structure contains only voxels truly belonging to that brain structure and it is not that important whether or not we detect all the voxels belonging to that brain structure.


## 3.2. PET Studies with Healthy Volunteers

We compared our automatic segmentation method to manual ROI segmentation. This comparison was done using11C-raclopride BPND high-resolution PET images. The PET scans were provided by Turku PET Center and were acquired from four right-handed, nonsmoking healthy male volunteers. This study was approved by the Joint Ethical Committee of the University of Turku and Turku University Central Hospital, and was conducted according to the Declaration of Helsinki. All subjects gave ethical committee-approved written informed consents. All subjects were free of any somatic or psychiatric illness, illicit drug abuse, and alcoholism, as confirmed by blood and urine tests, somatic assessment, and a structured clinical interview for DSM-IV Axis I disorders conducted by a psychiatrist. The age, height, and weight of the subjects were 27.5±7.74; (20.3 to 33.1) years, 175±1 (174 to 176) cm, and 76.1±5.16 (70.0 to 82.5) kg, respectively (mean, s.d.; range). All subjects underwent T1-weighted magnetic resonance imaging at 1.5T to rule out structural brain abnormalities.The radioligand11C-raclopride was prepared as previously described [15]. PET imaging was performed by a brain-dedicated ECAT HRRT, Siemens Medical Solutions, Knoxville, Tnee, USA) [12].  In the acquired measurements, a molded thermoplastic mask was used as a head fixation to reduce the noise due to head movements. The right antecubital vein was cannulated and 11C-raclopride was administered as an intravenous rapid bolus flushed with saline. Injected dose and the specific radioactivity were 443.25±104.69 MBq and 185.7±218.4 MBq/nmol, respectively (mean ± S.D.). The PET 11C-raclopride scans were acquired in list mode and histogrammed by axial compression of span 9. The 11C-raclopride uptake was measured for fifty one minutes after injection and the frame sequence consisted of three frames of sixty seconds, four frames of one hundred and eighty seconds, and six frames of three hundred and sixty seconds. The image size of each scan was 200×200×150 slices and they had an isotropic voxel dimension of 1.22×1.22×1.22 mm3. The raw data was reconstructed with a speed-optimized version of OP-OSEM-3D (Ordinary Poisson-OSEM in full 3D) [16], with sixteen subsets and eight full iterations [17]. The frame-images of each dynamic PET image were first realigned to correct for head movement during the scan.  Integral image (summed over time) was calculated from the realigned PET-image and the MRI was coregistered to this PET-image. All realignment and coregistration procedures were performed using Statistical Parametric Mapping software version 2 (SPM2, Wellcome Department of Cognitive Neurology, University College London, UK). The parametric BPND images were created from the dynamic images using the simplified reference tissue model with the cerebellar cortex as the reference region [14].Regions of interest (ROIs) were manually drawn to coregistered transaxial MRI using Imadeus software (version 1.2., Forima Inc., Turku, Finland). ROIs onto dorsal caudate, dorsal putamen, ventral striatum, thalamus, and cerebellum were defined as previously described [16]. Simplified reference tissue model with the cerebellar cortex as the reference region [14] was applied to TAC data to obtain the regional BPND values.Average NAS (ANAS) over the four subjects was used to measure the similarity of theBPND values based on the automatic segmentation and the manual segmentation. In this case the BPa and BPb in (5) represent, respectively, the BPND values computed based on the automatic segmentation and manual segmentation.  In this case, however, it cannot be claimed that the BPND values based on the manual segmentation would be correct although they can be used as benchmarks.


## 4. Results

### 4.1. Phantom Study

Table1 presents the BPND and NAS values from automatic segmentations of caudate, ventral striatum, putamen, anterior and posterior putamen. These results were obtained when applying the automatic method to ideal and realistic noise phantoms. The voxel connectivity threshold in (1) was set to 10 with the phantom images.Table 1
BPND values, NAS values, the number of voxels in each brain structure, and PPV values  computed based on the ideal phantom and the realistic noise phantom.
MeanBPND valueNAS# voxelsPPVIdeal phantomCaudateAverageAutomatic1.6140.983436720.8295AverageGround truth1.6416190Ventral striatumAverageAutomatic1.3510.981723250.0714AverageGround truth1.376954PutamenAverageAutomatic2.6980.989748740.9984AverageGround truth2.7266289Posterior PutamenAverageAutomatic2.5240.923026390.4433Averagegroung truth2.7261929Anterior PutamenAverageAutomatic2.9110.999322350.9817AverageGround truth2.9134360Realistic noise phantomCaudateAverageAutomatic1.6550.991543400.7005AverageGround truth1.6416190Ventral striatumAverageAutomatic1.4030.980623990.0354AverageGround truth1.376954PutamenAverageAutomatic2.8540.954150630.9891AverageGround truth2.7266289Posterior PutamenAverageAutomatic2.5910.949222710.5042AverageGround truth2.7261929Anterior PutamenAverageAutomatic2.8040.961927920.9362AverageGround truth2.9134360Table1 shows that the ideal phantom NAS results were between 92.3% and 99.9%, while the phantom with a realistic noise level had NAS values between 94.9% and 99.9%. In a test-retest setting, where the same subject is studied twice with a short interval between the scans, the similarity of BPND estimations is typically in the range of 90%–95% [3, 18]. Therefore, the automatic segmentations provided accurate BPND values in the phantom experiments, in that the deviation from the ground truth was less than the typical test-retest variability. Moreover, the NAS values presented in Table 1 also indicate that the increase of noise level did not affect the NAS values, suggesting that the automatic method was resistant to noise.  The PPV values were high for caudate and anterior putamen indicating that the functionally driven segmentations were also in an excellent agreement with the anatomy with respect to these structures. The PPV values for posterior putamen were lower resulting from the posterior putamen in automatic segmentation containing many partial volume voxels between posterior and anterior putamen, and also, between anterior putamen and surrounding white matter tissue. This was expected since the BPND values of the partial volume voxels between anterior putamen and surrounding white matter were similar to BPND values of voxels in posterior putamen. The PPV values for ventral striatum were low reflecting the mixing of ventral striatum and caudate in automatically segmented volumes. This was due to small size of the structure, 954 voxels when both left and right hemispheres were considered, majority of which were partial volume voxels (922 voxels had a neighboring voxels belonging to another structure) in the anatomical ground truth. However, while the segmentation of ventral striatum was anatomically incorrect, the BPND values obtained based on it were accurate. This was explained by the fact that the majority of the voxels incorrectly labeled as ventral striatum were partial volume voxels near the boundary of caudate and ventral striatum. Thus, we do not consider the low PPV values of ventral striatum particularly worrisome. Note also that we have taken nucleus accumbens to represent ventral striatum while ventral striatum also encompasses some adjacent areas.


### 4.2. PET Studies with Healthy Volunteers

Table2 presents the average BPND values obtained based on the automatic and the manual segmentation of the HRRT images of four healthy subjects. The average BPND values attained when applying manual segmentation to the real test data were 3.54±0.53 for caudate, 4.26±0.61 for putamen, and 3.42±0.59 for ventral striatum. Note that the BPND values are average values across hemispheres and subjects. The average BPND values obtained with the automatic method were 2.97±0.59 for caudate, 3.85±0.64 for putamen, and 3.24±0.57 for ventral striatum.Table 2
BPND  and ANAS values obtained from the average of the HRRT PET data. These measurements allow the similarity analysis between the automatic and manual segmentations.
mean±SDANASCaudateLeftManual3.56±0.460.78LeftAutomatic2.86±0.41RightManual3.53±0.590.84RightAutomatic3.07±0.79AverageManual3.54±0.530.81AverageAutomatic2.97±0.59PutamenLeftManual4.29±0.630.89LeftAutomatic3.86±0.60RightManual4.24±0.610.90RightAutomatic3.84±0.69AverageManual4.26±0.610.90AverageAutomatic3.85±0.64Ventral striatumLeftManual3.41±0.570.91LeftAutomatic3.16±0.57RightManual3.43±0.610.86RightAutomatic3.35±0.71AverageManual3.42±0.590.90AverageAutomatic3.24±0.57As can be seen from Table2, the BPND values for the automatic and manual segmentations were similar for all structures and ANAS varied from 78% to 91%. These ANAS values indicate that the variations in regional BPND due to ROI extraction method were small compared to variations due to other choices in the image processing (e.g., the reconstruction and the parametric image generation method). The ANAS values were smallest in caudate. By visually inspecting the automatic segmentation results, it could be seen that the caudate cluster did not only include caudate, but also some white matter. This was probably the reason for the lower BPND, and thereby, lowers NAS values for the caudate.Figure4 shows the result of applying the automatic segmentation method to both left and right striatum from one of the HRRT PET scans. The results are from slices 115 (a), 112 (b), 109 (c), 106 (d), and 103 (e) from the coronal view. The segmentation results were appropriate, namely, the algorithm found four clusters corresponding to anterior putamen, posterior putamen, ventral striatum, and caudate. In Figure 4 the orange and white clusters correspond, respectively, to the posterior and the anterior putamen, the red cluster corresponds to the caudate and, the yellow cluster corresponds to the ventral striatum and a part of the caudate head.Segmentation of the right and left striatum using the automatic method. Coronal (a) slice 115, (b) slice 112, (c) slice 109, (d) slice 106, and (e) slice 103 magnified to show only the striatum. The orange and white clusters correspond, respectively, to the posterior and the anterior putamen, the red cluster corresponds to the caudate and, the yellow cluster corresponds to the ventral striatum. It can be seen that ventral striatum contains also a part of the caudate head.
(a)(b)(c)(d)(e)The results presented in Figure5(a) are from slices 110, 115, and 120, from coronal, sagittal, and transverse views, while in Figure 5(b) they are from slice 110, 145, and 120. The images presented in Figures 4 and 5 are taken from the same clustering result.Segmentation result using the automatic method: (a) left striatum and (b) right striatum. In both images (a), (b), the orange and white clusters correspond, respectively, to the posterior and the anterior putamen, the red cluster corresponds to the caudate and, the yellow cluster corresponds to the ventral striatum and a part of the caudate head.
(a)(b)Interestingly, the automatic method was able to divide the putamen into two substructures. A slight difference betweenBPND values within anterior and posterior putamen was observed in this study.  Table 3 shows the BPND values from anterior and posterior putamen from all the four HRRT PET images. Unfortunately, the manual segmentation of these putaminal subdivisions was not available and we cannot compare the automatic BPND values with the manual ones. However, by inspecting Table 2, it can be seen that the similarity between the putamen BPND values obtained with the help of automatic and manual segmentations was approximately 90%. This indicates that the whole putamen was segmented very similarly with automatic and manual methods what reassures the correctness of the automatic segmentation. Moreover, when the method was applied to the phantoms, where the BPND values from those subdivisions were known, the algorithm accurately segmented both subdivisions (Table 1). These two facts allow us to assume that this separation was correct.Table 3
BPND values from anterior and posterior putamen from all 4 HRRT PET scans.
Binding potential valuesAnterior putamenImage 1Left3.5810Right3.2744Image 2Left3.8194Right3.7010Image 3Left3.8194Right3.5784Image 4Left4.4652Right4.7607Posterior putamenImage 1Left2.5193Right1.9575Image 2Left4.3798Right4.6810Image 3Left4.3961Right4.2066Image 4Left4.5961Right4.3590In Table3, BPND values from Image1 are lower than for other images. The same was visible from the whole putamen BPND values based on the manual segmentation: left and right putamen average BPND was 3.41 for Image1 while it was greater or equal than 4.27 for other images. In the visual inspection, the BPND Image 1 appeared considerably noisier compared to the other images featured in this study, and especially, posterior putamen had also lower BPND values compared to the other images. The reasons for this might include the low injected dose (292 MBq), imaging artifacts (e.g., subject movements), and individual variability.The final segmentation results were found to be highly dependent on the initial clustering guesses when applying the kernel weightedk-means algorithm to subdivide the striatum. The best segmentation results were obtained when the initial clusters were separated into four clusters each with the same size. A random initialization resulted usually in a single cluster consisting of the entire striatum. The neighborhood parameter in (1) needed also to be tuned to achieve the best segmentation. A voxel connectivity threshold of 4 voxel-sides achieved the best results.  This value was different from the phantom experiments due to more complex noise patterns with real data. It was also noted that the method always resulted in four clusters even if the number of the initial clusters was greater than four. In other words, starting with more than four clusters resulted in empty clusters.


## 4.1. Phantom Study

Table1 presents the BPND and NAS values from automatic segmentations of caudate, ventral striatum, putamen, anterior and posterior putamen. These results were obtained when applying the automatic method to ideal and realistic noise phantoms. The voxel connectivity threshold in (1) was set to 10 with the phantom images.Table 1
BPND values, NAS values, the number of voxels in each brain structure, and PPV values  computed based on the ideal phantom and the realistic noise phantom.
MeanBPND valueNAS# voxelsPPVIdeal phantomCaudateAverageAutomatic1.6140.983436720.8295AverageGround truth1.6416190Ventral striatumAverageAutomatic1.3510.981723250.0714AverageGround truth1.376954PutamenAverageAutomatic2.6980.989748740.9984AverageGround truth2.7266289Posterior PutamenAverageAutomatic2.5240.923026390.4433Averagegroung truth2.7261929Anterior PutamenAverageAutomatic2.9110.999322350.9817AverageGround truth2.9134360Realistic noise phantomCaudateAverageAutomatic1.6550.991543400.7005AverageGround truth1.6416190Ventral striatumAverageAutomatic1.4030.980623990.0354AverageGround truth1.376954PutamenAverageAutomatic2.8540.954150630.9891AverageGround truth2.7266289Posterior PutamenAverageAutomatic2.5910.949222710.5042AverageGround truth2.7261929Anterior PutamenAverageAutomatic2.8040.961927920.9362AverageGround truth2.9134360Table1 shows that the ideal phantom NAS results were between 92.3% and 99.9%, while the phantom with a realistic noise level had NAS values between 94.9% and 99.9%. In a test-retest setting, where the same subject is studied twice with a short interval between the scans, the similarity of BPND estimations is typically in the range of 90%–95% [3, 18]. Therefore, the automatic segmentations provided accurate BPND values in the phantom experiments, in that the deviation from the ground truth was less than the typical test-retest variability. Moreover, the NAS values presented in Table 1 also indicate that the increase of noise level did not affect the NAS values, suggesting that the automatic method was resistant to noise.  The PPV values were high for caudate and anterior putamen indicating that the functionally driven segmentations were also in an excellent agreement with the anatomy with respect to these structures. The PPV values for posterior putamen were lower resulting from the posterior putamen in automatic segmentation containing many partial volume voxels between posterior and anterior putamen, and also, between anterior putamen and surrounding white matter tissue. This was expected since the BPND values of the partial volume voxels between anterior putamen and surrounding white matter were similar to BPND values of voxels in posterior putamen. The PPV values for ventral striatum were low reflecting the mixing of ventral striatum and caudate in automatically segmented volumes. This was due to small size of the structure, 954 voxels when both left and right hemispheres were considered, majority of which were partial volume voxels (922 voxels had a neighboring voxels belonging to another structure) in the anatomical ground truth. However, while the segmentation of ventral striatum was anatomically incorrect, the BPND values obtained based on it were accurate. This was explained by the fact that the majority of the voxels incorrectly labeled as ventral striatum were partial volume voxels near the boundary of caudate and ventral striatum. Thus, we do not consider the low PPV values of ventral striatum particularly worrisome. Note also that we have taken nucleus accumbens to represent ventral striatum while ventral striatum also encompasses some adjacent areas.


## 4.2. PET Studies with Healthy Volunteers

Table2 presents the average BPND values obtained based on the automatic and the manual segmentation of the HRRT images of four healthy subjects. The average BPND values attained when applying manual segmentation to the real test data were 3.54±0.53 for caudate, 4.26±0.61 for putamen, and 3.42±0.59 for ventral striatum. Note that the BPND values are average values across hemispheres and subjects. The average BPND values obtained with the automatic method were 2.97±0.59 for caudate, 3.85±0.64 for putamen, and 3.24±0.57 for ventral striatum.Table 2
BPND  and ANAS values obtained from the average of the HRRT PET data. These measurements allow the similarity analysis between the automatic and manual segmentations.
mean±SDANASCaudateLeftManual3.56±0.460.78LeftAutomatic2.86±0.41RightManual3.53±0.590.84RightAutomatic3.07±0.79AverageManual3.54±0.530.81AverageAutomatic2.97±0.59PutamenLeftManual4.29±0.630.89LeftAutomatic3.86±0.60RightManual4.24±0.610.90RightAutomatic3.84±0.69AverageManual4.26±0.610.90AverageAutomatic3.85±0.64Ventral striatumLeftManual3.41±0.570.91LeftAutomatic3.16±0.57RightManual3.43±0.610.86RightAutomatic3.35±0.71AverageManual3.42±0.590.90AverageAutomatic3.24±0.57As can be seen from Table2, the BPND values for the automatic and manual segmentations were similar for all structures and ANAS varied from 78% to 91%. These ANAS values indicate that the variations in regional BPND due to ROI extraction method were small compared to variations due to other choices in the image processing (e.g., the reconstruction and the parametric image generation method). The ANAS values were smallest in caudate. By visually inspecting the automatic segmentation results, it could be seen that the caudate cluster did not only include caudate, but also some white matter. This was probably the reason for the lower BPND, and thereby, lowers NAS values for the caudate.Figure4 shows the result of applying the automatic segmentation method to both left and right striatum from one of the HRRT PET scans. The results are from slices 115 (a), 112 (b), 109 (c), 106 (d), and 103 (e) from the coronal view. The segmentation results were appropriate, namely, the algorithm found four clusters corresponding to anterior putamen, posterior putamen, ventral striatum, and caudate. In Figure 4 the orange and white clusters correspond, respectively, to the posterior and the anterior putamen, the red cluster corresponds to the caudate and, the yellow cluster corresponds to the ventral striatum and a part of the caudate head.Segmentation of the right and left striatum using the automatic method. Coronal (a) slice 115, (b) slice 112, (c) slice 109, (d) slice 106, and (e) slice 103 magnified to show only the striatum. The orange and white clusters correspond, respectively, to the posterior and the anterior putamen, the red cluster corresponds to the caudate and, the yellow cluster corresponds to the ventral striatum. It can be seen that ventral striatum contains also a part of the caudate head.
(a)(b)(c)(d)(e)The results presented in Figure5(a) are from slices 110, 115, and 120, from coronal, sagittal, and transverse views, while in Figure 5(b) they are from slice 110, 145, and 120. The images presented in Figures 4 and 5 are taken from the same clustering result.Segmentation result using the automatic method: (a) left striatum and (b) right striatum. In both images (a), (b), the orange and white clusters correspond, respectively, to the posterior and the anterior putamen, the red cluster corresponds to the caudate and, the yellow cluster corresponds to the ventral striatum and a part of the caudate head.
(a)(b)Interestingly, the automatic method was able to divide the putamen into two substructures. A slight difference betweenBPND values within anterior and posterior putamen was observed in this study.  Table 3 shows the BPND values from anterior and posterior putamen from all the four HRRT PET images. Unfortunately, the manual segmentation of these putaminal subdivisions was not available and we cannot compare the automatic BPND values with the manual ones. However, by inspecting Table 2, it can be seen that the similarity between the putamen BPND values obtained with the help of automatic and manual segmentations was approximately 90%. This indicates that the whole putamen was segmented very similarly with automatic and manual methods what reassures the correctness of the automatic segmentation. Moreover, when the method was applied to the phantoms, where the BPND values from those subdivisions were known, the algorithm accurately segmented both subdivisions (Table 1). These two facts allow us to assume that this separation was correct.Table 3
BPND values from anterior and posterior putamen from all 4 HRRT PET scans.
Binding potential valuesAnterior putamenImage 1Left3.5810Right3.2744Image 2Left3.8194Right3.7010Image 3Left3.8194Right3.5784Image 4Left4.4652Right4.7607Posterior putamenImage 1Left2.5193Right1.9575Image 2Left4.3798Right4.6810Image 3Left4.3961Right4.2066Image 4Left4.5961Right4.3590In Table3, BPND values from Image1 are lower than for other images. The same was visible from the whole putamen BPND values based on the manual segmentation: left and right putamen average BPND was 3.41 for Image1 while it was greater or equal than 4.27 for other images. In the visual inspection, the BPND Image 1 appeared considerably noisier compared to the other images featured in this study, and especially, posterior putamen had also lower BPND values compared to the other images. The reasons for this might include the low injected dose (292 MBq), imaging artifacts (e.g., subject movements), and individual variability.The final segmentation results were found to be highly dependent on the initial clustering guesses when applying the kernel weightedk-means algorithm to subdivide the striatum. The best segmentation results were obtained when the initial clusters were separated into four clusters each with the same size. A random initialization resulted usually in a single cluster consisting of the entire striatum. The neighborhood parameter in (1) needed also to be tuned to achieve the best segmentation. A voxel connectivity threshold of 4 voxel-sides achieved the best results.  This value was different from the phantom experiments due to more complex noise patterns with real data. It was also noted that the method always resulted in four clusters even if the number of the initial clusters was greater than four. In other words, starting with more than four clusters resulted in empty clusters.


## 5. Discussion

A new automatic method was developed to extract striatal brain structures, namely ventral striatum, caudate, anterior putamen and posterior putamen, from parametric HRRT11C-raclopride PET images. This anatomical subdivision of the striatum has relevant functional implications. An automatic method for the extraction of these structures from high-resolution PET images allows for inexpensive and reproducible extraction of the quantitative information from these images necessary in brain research and drug development. The method developed in this study proved to be accurate and reproducible, and was found to provide more information about the distribution of D2 receptors within the striatum than manual segmentation. We evaluated it using both images from healthy volunteers and simulated images (phantoms).The method is based on the extraction of the striatum using deformable surfaces. The extracted striatum is then sub-divided into four clusters corresponding to ventral striatum, caudate and anterior and posterior putamen using the weighted kernelk-means algorithm.  In this clustering, we take into account both the intensity value, in this case representing binding potential, of the voxel as well as its spatial location with respect to the other striatal voxels. We chose to use the weighted kernel k-means criterion for the clustering because of its recently established equivalence to the normalized cuts criterion. The weighted kernel k-means algorithm can be seen as a fast algorithm to optimize the multiway normalized cuts criterion. Earlier methods to optimize this criterion require a solution of eigenvectors of a large matrix, which could be computationally difficult. Indeed, the new algorithm, programmed in MATLAB, required only few minutes to segment the striatum.In the algorithm, a few parameter values needed to be set (see Sections2.2 and 2.3). All the parameter values were set in the same way in all the experiments with a single exception of the connectivity threshold set differently between the phantom and the real data experiments. The reason for this was the obvious difference between the phantoms and real PET studies. The noise model in the phantom studies was highly simplified and it either ignored or simplified many aspects that affect the appearance of parametric BPND images in reality. Also, the spatial resolution of the simulated phantom images was higher than in the actual studies because the effects of the image reconstruction were ignored. Simulating more realistic computer generated PET images requires advanced computational methods (see, e.g., [19]) and these Monte Carlo simulators are not yet directly available for modeling HRRT-PET images. We do not think that the setting of the parameter values for the method would present a major hurdle to the automatic application of the method to a larger set of HRRT-PET images.In the experiments with the simulated images, the numbers of voxels in the substructures in the automatic segmentations differed from the numbers of voxels in the anatomical ground truth. In the anterior putamen and caudate, the numbers of voxels in the automatic segmentation were lower than in the ground truth. This is a desired feature of the segmentation algorithm because otherwise partial volume effects (PVEs)—here referring to the combined effect of the tissue fraction effect and the point spread effect in the terminology of [20]—could compromise the quantitative accuracy of the regional BPND values (see discussions on this issue in [3, 4]). While the resolution of PET images acquired by the HRRT is superior to the resolution of PET images acquired by earlier PET scanners, the image resolution and PVEs still continue to be a limiting factors of the accuracy of the image  segmentation. The numbers of voxels in the ventral striatum and posterior putamen were greater in the automatic segmentation than in the anatomical ground truth. We attribute also this to the PVEs although it was not a desired feature of the segmentation algorithm. In the case of ventral striatum, the majority of extra voxels were partial volume voxels near the boundary of caudate and ventral striatum, thus having lower BPND values than “pure” caudate voxels. Similarly, the majority of extra voxels assigned to posterior putamen were partial volume voxels between the anterior and posterior putamen. Hence, the effects of these segmentation errors on the accuracy of the regional BPND values were negligible as it visible in NAS values in Table 1.Although the accuracy or regionalBPND values—the computation of which was the reason behind designing the segmentation algorithm—were accurate it would be advantageous if the segmentations themselves were more accurate, especially with respect to ventral striatum. There are a number of research avenues one could follow while trying to improve the algorithm. The extraction of striatum could be improved by taking the PVEs better into account. This could be achieved by adopting a volume based segmentation strategy based on, for example, Markov Random Fields, instead of the present edge-based strategy. However, this is not straight-forward because the BPND values within the striatum are highly nonuniform due to BPND differences between striatal substructures, noise, and PVEs. Moreover, in the case of the HRRT scanner, the PVE is challenging to model due to gaps between adjacent flat panel detector heads.   The segmentation of striatum could also include modeling of the PVE and neuroanatomical a priori knowledge could be better utilized either in the kernel weighted k-means criterion or in the initialization for the kernel weighted k-means algorithm.The ability of the automated method to divide the putamen into anterior and posterior components was an originally unexpected result. This ability is truly interesting as the distribution ofD2 dopamine binding within the putamen has been rarely studied. This is partly due to the fact that with the earlier generation PET scanners differences between anterior and posterior putamen are less evident due to low image resolution. Also, there are no generally accepted guidelines for dividing the putamen between anterior and posterior components manually; for instance, we previously employed a crude division half-way along the anterior-posterior axis of the putamen [10]. For this reason, we did not attempt putaminal subdivisions based on the manual segmentation a priori. Instead, we chose to evaluate the automatic segmentation results based on visual inspection, which indicated that the division was indeed highly reasonable. Moreover, in the studies with simulated images, the quantitative values obtained for the putamen subdivisions (anterior and posterior putamen) based on automatic segmentation were highly similar to the ground truth. The division to anterior and posterior putamen is consistent with a decreasing rostrocaudal gradient in dopamine D1/D2 receptor density in the human striatum [10], and may be physiologically relevant. For instance, D2 receptors are inhibitory on striatal output neurons, and differential expression of D2 receptor in functional subdivisions of the putamen may underlie differential dopaminergic modulation of associative versus motor functions [5, 6]. The fact that the direction of the contrast was not uniform across the subjects in this small sample (i.e., some had decreasing and some had increasing anterior-posterior gradients) might suggest inter-individual differences in preferential dopaminergic modulation of associate versus motor functions, but larger samples would be required to confirm this.This method for automatic extraction of volumes could be further developed to segmentation of other types of PET images. Especially, the method is applicable to extract striatal substructures from PET images acquired using other radiopharmaceuticals than11C-raclopride if the images share the characteristics of the 11C-raclopride images, that is, the image intensity is the highest in striatum. An example of such radiopharmaceutical is 18F-SPA-RQ used to study brain neurokinin NK1 receptors. The automatic segmentation of low resolution 18F-SPA-RQ PET images was studied in [4].It would be conceptually straight-forward to modify the proposed method to work directly with 4-D PET images, where each voxel would be represented by a time activity curve—instead of 3D parametric images where each voxel is represented by aBPND value computed based on the time activity curve. Technically, this would entail replacing the term |I(i)-I(j)|/maxk∈DI(k) in (1) by an appropriate distance measure between two time activity curves. The distance measure could be, for example, the Euclidean distance between the vector representations of radioactivity measurements in different time intervals or a more sophisticated L2-norm between the time activity curves viewed as continuous valued functions in a suitable function space. However, the specification of the meaningful distance measure might prove troublesome in practice. In the current work, our interest was in the segmentation of striatum based on BPND values, and other information derivable from time activity curves was therefore not relevant for this application. We would not rule out the possibility to segment striatum without BPND computation as a preprocessing step, however, the development of such a method could be challenging in practice.



---

*Source: 156234-2009-11-04.xml*