Self-normalized Classification of Parkinson’s Disease DaTscan Images

A. PET/SPECT Normalization

PET/SPECT imaging can be categorized into dynamic and static imaging. Dynamic imaging acquires images at multiple time points after tracer injection and is mainly used in research. Static imaging acquires a single image over a fixed period hence is suitable for clinical applications. Analyzing clinical (static) PET/SPECT images typically requires a preprocessing step which normalizes its intensity [1], [2]. Normalization is required because the amount of radioligand reaching the brain depends on multiple factors, including age, sex, and medication. This results in an unknown scaling factor for each image. Standard normalization selects a normalization region to calculate the BP, thereby eliminating this scaling factor, as mentioned before. Note that in PD DaTscans, the BP is called the striatal binding ratio (SBR) since dopamine transporters mainly exist in the striatum.

In PD DaTscans, the occipital lobe is usually selected as the normalization region since it contains few dopamine transporters [9], [10], [11] and most of the binding of the radiotracer in this region is non-specific. Other choices for the normalization region include the cerebellum [12] and the whole brain except the striatum [10]. Different normalization regions alter the predictive power of BP to classify PD, up to a difference of 0.147 in terms of area under the curve [5].

For PET/SPECT imaging using other tracers, e.g. [18 F]fluorodeoxiglucose for visualizing glucose (FDG-PET), the impact of the normalization regions is also significant [7], [8].

Besides the BP normalization mentioned above, other normalization techniques have also been proposed. For DaTscans, these include analyzing the distribution of intensity values in the whole brain except the striatum [13], and minimizing the squared error within a selected region between a template and the linear transformed image [1].

B. DaTscan Classification

Classification of DaTscan images into PD and HC has been achieved using standard machine learning techniques [14], [15], [10], [11] as well as deep learning methods [16]. Early classification work tends to use a small proprietary datasets and perform their own registration [17], [14], while later work has shifted to using a public dataset: the Parkinson’s Progression Marker Initiative (PPMI) dataset which provides already registered DaTscans [15], [11], [18].

Support vector machine (SVM) and logistic regression are probably the most popular PD vs. HC classifiers, appearing in most of the non-deep-learning studies [14], [15], [10], [11], [18], [19]. Graph-based transductive learning was introduced for classifying multi-modality neurodegenerative image data in [20]. Recently, convolutional neural networks (CNN) have been applied to classifying DaTscan images for PD diagnosis [21], [16], [22]. Most of these methods calculate the SBR using a normalization region as a pre-processing step before classification [15], [11], [18].

Other studies use geometric image features including the length and volume of the segmented striatum [23], shape fitting coefficients [24], [25], isosurfaces [16], and intensity summary statistics [22].

A. Multiplicative Equivalence

Let Ω be the set of voxels in an image, with the total number of voxels being d. Any nonzero image I defined on Ω is an element of ${ℝ}^d$. SPECT/PET images have non-negative voxels, i.e. these images lie in the non-negative orthant of ${ℝ}^d$ with the origin removed. Geometrically speaking, the set of all images related by positive scalar multiples is a half-ray passing through the origin of ${ℝ}^d$ (see Fig. 2(a)). We denote the half ray of the image I by [I]. We also denote the set of all half-rays passing through the non-negative orthant of ${ℝ}^d$ (with the origin removed) as $\mathcal{I}$. Classifying images under multiplicative equivalence means partitioning $\mathcal{I}$ into disjoint subsets.

B. Normalized-Classification

A normalization region is a set of voxels N ⊂ Ω. All images restricted to N form a subspace of ${ℝ}^d$, which we denote as Π_N. The projection operator ${\pi }_N:{ℝ}^d\to {\Pi }_N$ projects an image onto this subspace by setting values outside N to zero.

Let $1={[1,1,\dots ,1]}^T\in {ℝ}^d$ denote an image with all 1’s. Then π_N(1) is the image containing all 1’s in the voxels N and 0’s everywhere else. For any image I, the mean value in the normalization region is $\frac{1}{|N|}{\left({\pi }_N(1)\right)}^{T}I$, where |N| is the number of voxels in N. Denoting $\frac{1}{|N|}{\pi }_N(1)$ by 1_N, the mean value in the normalization region is simply $1_N^{T}I$ and the normalized image is $\tilde{I}=I/1_N^{T}I$. Normalization can be interpreted geometrically after noting that the mean value of $\tilde{I}$ in the normalization region is always 1, since $1_N^T\tilde{I}=1_N^T\left(I/1_N^{T}I\right)=1$. Thus the geometry of normalization is explained as follows (see Fig. 2(b)): From the image I construct the half-ray [I]. Then the normalized image $\tilde{I}$ is the intersection of the half ray [I] with ${\mathcal{A}}_N$, the d−1 dimensional affine subspace defined by ${\mathcal{A}}_N=\left\{x\in {ℝ}^d:1_N^{T}x=1\right\}$.

Once the image is normalized, it is classified using voxels in another region C ⊂ Ω. We assume that N ∩C = ∅, i.e. the normalization and classification voxels are disjoint. Similar to N above, images restricted to C form the subspace Π_C of ${ℝ}^d$. Because N and C are disjoint, we have Π_C ⊥ Π_N.

A linear classifier chooses a unit norm vector w ∈ Π_C and a real number b to define a decision boundary $w^T\tilde{I}-b=0$. We refer to linear classification using the normalized image $\tilde{I}$ as normalized-classification, and its boundary as the normalized-classification boundary. This boundary has an alternate description. First consider the equation w^Tx − b = 0 for $x\in {ℝ}^d$. Because w ≠ 0, this equation describes a d − 1 dimensional affine subspace of ${ℝ}^d$. Call it $\mathcal{B}$ (see Fig. 2(b)). Then the normalized-classification boundary $w^T\tilde{I}-b=0$ is the intersection of $\mathcal{B}$ and ${\mathcal{A}}_N$ (see Fig. 2(b)), which is the set of all $x\in {ℝ}^d$ satisfying

Because w and 1_N are linearly independent, $\mathcal{B}\cap {\mathcal{A}}_N$ is a d−2 dimensional affine subspace of ${ℝ}^d$.

C. From Normalized-Classification to Subspace Classification

A slight shift in point-of-view shows that normalized-classification is really just a classification of half-rays by subspaces. The key idea here is to take $\mathcal{D}=\text{span}\left(\mathcal{B}\cap {\mathcal{A}}_N\right)$ (see Fig. 2(c)). Because $\mathcal{D}$ is a span of a set of points, $\mathcal{D}$ is a subspace of ${ℝ}^d$. It has properties that are described below. Proofs of all properties can be found in the Appendix.

D. Self-normalized Classification

The above analysis strongly suggests how the dependence on the normalization region can be eliminated altogether: Take the intersection of each half-ray and the unit sphere in ${ℝ}^d$ (see Fig. 2(d)) and then classify the intersection points with a subspace. There is no normalization region involved in this; the data normalizes itself — it is self-normalizing.

This idea can be pushed further in two ways: First, we need not classify the intersection points on the sphere with only a subspace. We can use any d−1 dimensional affine subspace. Second, we can apply the same idea to the image restricted to the classification voxels C. That is, we take only intensities of voxels in C, project them to the unit sphere, and classify them using an affine subspace. This too eliminates the normalization region.

E. Non-linear Classification

The above ideas extend easily to nonlinear classification. To see how, first note that the normalization step is the same as before. It corresponds to moving the image I to the intersection of the half-ray [I] and ${\mathcal{A}}_N$. However, the classification boundary, which was the affine subspace $\mathcal{B}\cap {\mathcal{A}}_N$ in Fig. 2, is no longer linear. Suppose it is a d−2 dimensional sub-manifold of ${\mathcal{A}}_N$, which we will call $\mathcal{D}$ (see Fig. 2(e)). Every $x\in \mathcal{D}$ gives a half-ray [x] and this half ray intersects the unit sphere at a single point. Because $\mathcal{D}$ is a d−2 dimensional submanifold, the set of all half rays of points in $\mathcal{D}$ is a d − 1 dimensional submanifold of ${ℝ}^d$, and this manifold intersects the unit sphere transversely to give a d − 2 dimensional submanifold of the unit sphere. Denote this submanifold as $\widehat{\mathcal{D}}$ (see Fig. 2(e)). It is straightforward to see that the partition (i.e. classification) of $\mathcal{I}$ by $\mathcal{D}$ is identical to the partition of $\mathcal{I}$ by $\widehat{\mathcal{D}}$. The converse is also straightforward: By simply reversing the argument it is easy to see that any classification boundary, which is a d − 2 dimensional submanifold of the part of the unit sphere in the non-negative orthant of ${ℝ}^d$, induces a d−2 dimensinal submanifold as a classification boundary in ${\mathcal{A}}_N$. Thus, we have

A. PPMI Data

The PPMI dataset contains 449 early-stage PD subjects and 210 HC subjects. The PD subjects have scans at baseline, and at approximately 1, 2, 4, and 5 years from baseline, with missing scans. Most of the HC subjects have only a single scan. The images have a size of 109 × 91 × 91 voxels, with 2 mm³ voxels. The images are already registered by PPMI to the Montreal Neurological Institute (MNI) atlas. However, following the procedure in [26], we found and removed some misregistered images. This left us with 365 PD subjects (ages: 62.6 ± 9.8 years, male/female: 237/128) and 208 HC subjects (ages: 60.6 ± 11.2 years, male/female: 135/73) to analyze. All PD subjects had baseline images, but only 136 PD subjects had images at year 4.

Because PD affects the two brain hemispheres asymmetrically, the DaTscan images were flipped around the mid-plane so that the more affected side was on the right. For normalized-classification, we used the occipital lobe as the normalization region (N) and the striatum as the classification region (C). The striatum mask was derived by applying Otsu’s threshold [27] on the mean HC image to remove the background and then applying it again on the remaining voxels to remove the nonspecific binding voxels. The occipital lobe mask was taken from [26]. All masks were restricted to the 29–55th slices.

B. Experimental Setup

The ultimate goal of our experiments is to verify the classification accuracy of the self-normalized classification of Section III-D. However, an important subgoal is to determine whether self-normalized classification provided useful information about voxels that are important to classification. This can be achieved by analyzing the weights for linear classifiers and the saliency maps [28] for non-linear classifiers.

We carried out two classifications: 1) HC vs. PD using only the baseline images, 2) Baseline PDs vs. PDs at year 4. The PD vs. HC classification has obvious clinical importance. The baseline vs year 4 classification is not clinically useful on its own, but because it identifies voxels where the disease progresses from baseline to year 4, it provides a simple disease progression footprint.

In each classification problem, we used three normalization strategies:

1. classical normalization using SBR with occipital lobe normalization (referred to as SBR from now on),

2. self-normalization via a projection of all voxels from the striatum and the occipital lobe onto the unit sphere (referred to as S + O),

3. self-normalization via a projection of voxels only from the striatum to the unit sphere (referred to as S).

To avoid numerical underflow problems, the radius of the sphere used in self-normalization was set to $\sqrt{d}$, where d is the dimension (number of voxels) of the region being projected. The dimension d for the above three strategies are d = 9948, 43077, and 9948 respectively.

Along with each normalization, we used three classifiers, two of which were linear: logistic regression (LR) and SVM, and one of which was nonlinear: convolutional neural network (CNN). The linear classifiers had a sparsity constraint as implemented in the fitclinear function in Matlab with a sparsa optimizer [29]. The nonlinear classifier was implemented in PyTorch with 2 convolutional/pooling layers (kernel size 5 × 5 × 5, 6 and 16 feature maps respectively) followed by 2 linear layers (120 hidden nodes). We used the ReLu activation function after the convolutional/linear layers. For the CNN, voxels were inserted into a 3D image cube padded by zeros. The combination of normalization method with the classifier led to 3 × 3 = 9 classification experiments, whose results are reported below.

For every classification task, we randomly split the data set 100 times into a training set (80%) and a test set (20%). The sparsity parameter for the linear classifiers was chosen using 10-fold cross validation on the 100 training sets. Fig. 3 shows the cross validation results of PD vs. HC. Similar curves (not shown) were obtained for PD baseline vs. PD year 4 classification. Because CNN training is computationally expensive, we did not use cross validation to set the hyperparameters for the nonlinear classifier. Instead, the learning rate was set manually to 0.01 and the number of epochs was set to 500. We used a scheduler which reduced the learning rate by half after half of the epochs. Furthermore, 10% of the training set was kept as the validation set and the remaining was used to train the neural network. We pick the network parameters with the highest accuracy (over the epochs) on the validation set for test.

C. Classification Accuracy

The classification accuracies of the 100 training/test sets for PD vs. HC and PD baseline vs. PD year 4 are shown in Fig. 4. For PD vs. HC, classical SBR (SBR), the self-normalized striatum plus occipital lobe voxels (S + O), and the self-normalized striatum voxels (S) lead to almost identical results for the test set (see Fig. 4(a) and Table I). The classification accuracies for the test set are quite high (the classification accuracies for the training set are similar). In fact, the classification accuracies in Table I are noticeably higher than those reported in the literature that use the same PPMI data (typical reported accuracies are in the range 95.1–97.9%) [15], [11], [19], [23], [16].

To compare the classification accuracies of standard SBR vs. self-normalized methods, we calculated p-values from a t-test comparing the mean classification accuracy of the self-normalizing methods to that of standard SBR. The p-values (see Table I) are all significantly above 0.05, showing that there is no significant performance difference between the methods.

Classification accuracies for PD baseline vs. PD year 4 show a similar pattern (see Fig. 4(b) and Table II). While the classification accuracies are not as high as PD vs. HC (this is discussed further in Section V), the differences in the performance of the classifiers are again not statistically significant, except for one case (LR on S + O).

D. Salient Voxels

Recall that part of our goal is to identify salient voxels (voxels which contribute significantly to classification). As mentioned before, classification weights of linear classifiers indicate these voxels. For nonlinear classifiers, the saliency map [28], which calculates the derivative of the logit with respect to the input, can be used for similar purposes. Because the saliency map is calculated for each test example, we took one train-test split and averaged the saliency maps over all the test examples.

The linear weights and the averaged saliency maps are shown in Fig. 5. Over the 9 combinations of normalizing strategies and classifiers used for each classification task, the results are surprisingly consistent. For HC vs. PD, the negative classifier weights and negative saliency (both rendered in blue) for PD vs. HC are mostly in the putamen on the right side (the more affected putamen). The positive weights and positive saliency is mostly in the left caudate (the least affected side). Thus reduced values in the right putamen relative to the values in the left caudate are significant in classifying PD vs. HC.

PD baseline vs. PD year 4 shows negative coefficients and negative saliency in the putamen on the left side, showing that decreasing values in this putamen corresponds to progression from baseline to year 4.

To make the above observations more quantitative, we evaluated the contribution of each region to the classification. Noting that all classifiers give almost identical salient voxels, we focused on the sparse logistic linear classifier. We restricted w to each of the four caudate and putamen regions from the MNI atlas, and evaluated $w^T\widehat{I}$ using the restricted w and self-normalized $\widehat{I}$ for PD vs. HC and PD baseline vs. PD year 4. For PD vs. HC, the mean $w^T\widehat{I}$ has the largest difference when w is restricted to the right putamen, indicating that this is the region where PD has the largest effect (see Fig. 6). Similarly, for PD baseline vs. PD year 4, we evaluated the mean for PD baseline and year 4. The largest change in the mean is in the left putamen (see Fig. 6).