Research Article - Journal of Biomedical Imaging and Bioengineering (2018) Volume 2, Issue 2
Times series evaluation of fine structure in neurological images.
2Department of Neurosurgery, RUSH University Medical Center
- *Corresponding Author:
- Richard D. Penn
Department of Bioengineering University of Illinois, Chicago
E-mail: [email protected]
Accepted date: October 25, 2018
Citation: Penn RD, Neill W, Werner MC. Times series evaluation of fine structure in neurological images. J Biomed Imag Bioeng. 2018;2(2): 102-106.
Image Modeling and Classification, Regression Analysis, MRI, Dementia
ADNI: Alzheimer’s disease Neuroimaging Initiative; MMSE: Mini Mental State Exam; OASIS: Open Access Series of Imaging Studies; OLS: Ordinary Least Square regression; KST: Kolmogorov-Smirnov Test; ROC: Receiver Operating Characteristic
Background and Purpose
When viewing a medical image one typically makes an overall assessment of the image and surveys region by region looking for elements of normal structure and pathology. The eye cannot easily view the tens of thousands of small regions of the image and compare them to each other. The information in the fine structure of the image is lost and neglected. Whether such fine structure of closely spaced pixels will be useful for medical purposes is an open question since it has been difficult to test images at this fine level.
We have developed a new time series -like computational method based on ordinary least squares regression and show that, indeed, the fine structure that is present in an image can be used to classify patients into mutually exclusive groups. A detailed mathematical explanation of the method has been published . This paper will outline the principles underlying the method, apply it to MR images of Alzheimer patients, investigate classification dependence on image resolution, and provide an example of how the analysis can be extended to the three-dimensional information in a stack of brain scan images. Time series typically refer to discrete measurements of a variable over time such as the Dow Jones industrial average, average temperature at a given location, or the output of a factory. Mathematical analysis of the changes or trends in the series of numbers is used to construct a model of the ordered numbers. A time series model of such data can be used to make predictions of future trends or to characterize a particular process . Given the known data, the model makes predictions of what the Dow will be in the future, what the next days’ temperature may be, or the possible output of a factory. The analysis of the time series can also be used to dissect out the important factors that influence the series over time and how many factors are necessary to describe the time series to a specified degree of accuracy .
The linear streams of numbers in the times series seem far removed from a two-dimensional array of pixel values of a typical MR image. However, it is easy to see that if a stacking procedure is used the columns of the image matrix can be arranged one after another in a long linear output of pixel values. A 208 x 176 matrix image would be arranged into a 1- dimensional array of 36,608 (as 208x176= 36,608) pixels in length (figure 1).
Figure 1: This figure illustrates how stacking image columns of pixels produces a single column that can be used in a regression model of the image. The image shown is 208 by 176 pixels, thus the entire image would produce a single column of 208x176=36608 pixels. (A): MR scan image of a subject classified as demented. The 2 vertical lines define which pixel intensities will be joined as a column. (B): The column of 416 pixel intensities made up from the first vertical line joined by the second vertical line to produce a time series-like data set of intensities indexed by image row indices.
Looking at such a sequence of pixel values, it is not obvious how such a string could be analyzed in a useful way. The key steps in our modeling of an image are to create one dimensional array of images whose pixels are offset by a step in the column direction, a step in the row direction, and a step in both directions. These steps, or lags in space, are directly analogous to the temporal lags in the time series modeling of temporal signals. Each of these three different strings of numbers provides information about how the image varies if the shift is row-wise, column-wise or both.
We show in  this image modeling procedure is equivalent to an hypothesis the image in question is the solution of a first order, two dimensional, linear, partial difference equation and that the hypothesis can be tested by using linear regression to estimate an image’s B values (regression coefficients) .
The way the structural analysis works can be illustrated by considering each step in the process using an MR image (Figure 2). If the full image is shifted up, over or diagonally by one pixel, visually the images appear similar, 2a. However if a region of 4 x 4 pixels is enlarged and then shifted, 2b, the magnified regions look quite different. The long columns made from stacking these pixels are also clearly different .
Figure 2: Steps used to stack an image. A: The original image, left, and the image moved up, over, and up and over one pixel. B: An enlarged area of the scan of 3 by3 pixels showing the difference between each due to the shift up, over, and up and over. C: The stacking of columns for each image. Note that the one pixel shifts which are hard to see in the top row, A , have a large effect on the values of the constructed columns and this reflects the fine structure of the original image.
The next step is to create a model of the original image by combining the information from the three columns, up, over, and diagonal to best approximate the original column of numbers. Each of these three columns is weighted mathematically to optimize the match to the original using a least squares method. This is done using linear algebra programs in Matlab. The weight of each column is its B value and the group of three B values is used to categorize each scan. This process takes only milliseconds to compute. The next step is to see if the B values of the scans can be used to separate subgroups that may be different in fine structure. The B values, which by design are sensitive to fine structure of the image, are used to separate one group of images from another. In the example to be analyzed, MR images of normal and demented patients are used.
Sample training scan sets of subjects of known cognitive state are used to calculate each subject’s B values. These values produce a single number (called a z projection) for each subject using a Fisher linear discriminant that optimally separates the demented and non-demented z projections . Each subject’s z projection determines subject cognitive group membership. If there is significant class difference between z projections, the image groups will separate into distinguishable distributions. Then, any scan not from the training sets can be modeled and its B values projected to determine class membership.
Since image fine structure is not visually apparent, a quantitative means would be useful to distinguish classes. This is what we have done with high probability for scans of normal and Alzheimer patients available from the Internet.
Materials and Methods
Our small sample study employed OASIS (Open Access Series of Imaging Studies) axial MR image slices (www.oasis-brains.org). The normal subjects were screened as right-handed males of ages 65 to 75 years with MMSE scores of 29 to 30 and the database produced 14 such images. Screening Alzheimer patients who were right-handed males of ages 67 to 84 years with MMSE scores of 17 to 26 produced 12 images.
Representative images are seen in Figure 3a. OASIS provided standardization for brain landmark location.
Our large sample study employed the larger 416 subjects and patient base of ADNI (Alzheimer’s disease Neuroimaging Initiative) axial MR image slices. (www.adni.loni.usc.edu). The 50 normal subjects were screened as males, ages 65 to 84 years with MMSE scores of 29 to 30 and the 47 demented patients were of ages 65 to 84 years with MMSE scores of 17 to 24. Representative images are seen in Figure 3b.
Figure 3b: Representative images from the Alzheimer’s disease Neuroimaging Initiative (ADNI). The top row shows normal images, the bottom row shows demented images.
Axial images were chosen which were one slice above the slice that showed the anterior commissure. In both studies the presented images were cropped of the cranium, dura, and artifacts manually.
MR Imaging Protocols
The OASIS images for each subject were a single axial slice, “snapshot” from the cross-sectional data. All the images were taken on a 1.5 Tesla Siemens Vision scanner using a MP RAGE sequence. The scans were taken in a single imaging session and a movement correlated co-registration average image was created. The ADNI images were on 1.5 Tesla machines and also used a MP RAGE sequence. The different centers had different MRI scanners that which were made by General Electric Healthcare, Phillips Medical Systems and Siemens Medical Solutions. The distribution of normal and demented patients for each scanner type was close to that of the overall group. The scans used were chosen so that no one single type of machine was used predominately for demented or normal patients.
Image Models and Discrimination Technique
The OASIS and ADNI images were of sufficient complexity as to support an expanded parameter set of 8 B values resulting from one and two shifts in rows, columns, and combined row and columns of the original scan. Additionally a column of constants added an additional parameter to account for image intensity mean values. Ordinary least squares regression estimated these 9 parameters (also called a subject’s parameter vector) using the shifted image columns and rows such that the variance between image and the explanatory column variables was minimized. The explanatory shifted columns and rows multiplied by their estimated B values are called the image model. Since the OASIS data set consisted of 14 normal and 12 demented subjects and each subject produced a regression-estimated 9 element parameter vector, the normal class parameter matrix containing all of the parameters estimated for normal subjects was 14 rows by 9 columns. Similarly, the demented parameter matrix was 12 rows by 9 columns. Comparable data matrices for the ADNI MR scan models were 50 by 9 and 47 by 9.
The demented and normal class parameter matrices were inputs to a Fisher linear discriminator that optimally projects each subject’s parameter vector onto a scalar projection axis. In this way each subject is assigned a unique single number called the subject’s z-projection number. That number is an optimal representation of the subject’s MRI scan because the Fisher linear discriminator maximally separates the normal subject projection numbers from the demented subject numbers . We used the Fisher optimization over both classes to produce a single 9 element vector which projects both the normal and demented subject parameter vectors onto the scalar z-projection axis. True classified and misclassified subjects generate the “true positive’’ and “false positive” events of a zero threshold ROC analysis . Therefore, if a subject outside of the database is to be tested for class membership, this vector tests that hypothesis by projecting the subject’s parameter onto the z-projection axis.
An experienced neuroradiologist was asked to separate our OASIS MR images into normal and demented groups. The MR images were provided in randomized order for visual assessment into their respective groups. No specific patient data was provided and the full series of image slices for each patient was not available for the neuroradiologist.
Time Series Method for 3D Image Models
As noted above, we employ linear regression to estimate B values for each shifted version of an original image. Each B value is the coefficient of a column of pixels representing a shifted version of the original image. This format suggests the hypothesis that a model of an original image could also be dependent on brain scans immediately above and below the image being modeled since neurons are not restricted to a 1 mm brain slice. If the hypothesis is true then the B values estimated for those columns representing shifted images above and below the original image should be significantly different from zero.
To test this hypothesis we made 3D models of the 50 normal ADNI images. Since adjacent images in this data set are not physically related, the B parameters relating the original image and the images above and below the original should not be significantly different from zero. Then we 3D modeled 60 axial slice images from a single subject. In this case the B parameters relating the original and adjacent images should be significantly different from zero. 3D axial image models incorporate sagittal and coronal information and weight the 3D information equitably only if the row, column, and row-column shifts are equal in number. Therefore for the original image to be modeled, we chose 3 parameters representing a column shift, a row shift, and a combined row-column shift. For the 2 adjacent slices 4 parameters sufficed as the unshifted images themselves are related to the original image so each 3D model produced B vectors with 11 estimated parameters.
A comparison of the ability of 2D and 3D models to classify images was performed with ADNI images, 20 demented and 20 normal slices from single subjects. All were estimated as 3D models. The same images were then estimated as 2D models and class discrimination computed with the Fisher linear discriminant.
When asked to place scans into normal or demented groups, an experienced neuroradiologist correctly identified 6/12 demented and 8/14 non-demented. This is statistically no better than guessing.
The MR OASIS images were 208 by 176 pixels and the ADNI images were 256 by 170 pixels that gave regression data column sizes of 36608 pixels and 43520 pixels respectively. This works out to 4068 and 4836 sample pixels per estimated regression parameter and is one reason why all of the estimated models had R-Square statistics of 99% or better. Further, all z projection distributions were found to be robustly normally distributed by the classical Kolmogorov-Smirnov test (KST) . Consequently, 25/26 or 96% of the OASIS images and 92/97 =94.8% ADNI images were correctly classified, see Figure 4.
Figure 4: (A) OASIS image classification distributions. The probability distributions derived from the z projections passed a rigorous test for normality. One image out of 2 was misclassified. (B) These image classification probability distributions are derived as above for the ADNI images. As above, the distributions passed a rigorous test for normality. Five images out of 97 were misclassified.
Figure 5: illustrates a test with 50 normal ADNI subjects (48 useable) of the hypothesis that the 2D image modeling technique can be extended to 3D. The graph indicates that none of the parameters expressing correlations between original and adjacent images are significantly different from zero. This result is consistent with the fact that the images are of distinct, normal subjects. Figure 6 shows the 11 parameters estimated for 60 normal ADNI subjects (58 useable) where the original and adjacent images are from the same subject. That virtually all estimated parameters are significantly different from zero supports the 3D hypothesis.
We also performed a cross validation test for the ADNI images by randomly selecting 12 normal and 12 demented subjects to be test subjects for a training classifier based on the remaining 38 normal and 35 demented subjects . The test subjects were classified with an accuracy of 22/24=91.7%.
Figure 7 and 8 illustrate the how estimating images by 2D and 3D differ in classification accuracy. 20 each of ADNI demented and normal subject images were modeled and classified as 2D images, Figure 7, and as 3D images, Figure 8. 3D imaging reduced incorrect classification from 11 to 1.
Our new technique using the fine structure of MRI images for classification is able to separate scans of non-demented and demented patients with high accuracy. It is more accurate than that of an expert neuroradiologist using the same single two dimensional images. This is not to say that our expert would not be more accurate if given the patient’s history and access to the full complement of scans for the patient but it illustrates how important information about fine structure cannot be seen even by a trained observer. Of course, additional information could also improve our classification accuracy. The addition of information from the scans below and above shows that three-dimensional information may be useful as 3D scans depend on sagittal and coronal information as well axial information in estimating model parameters. This suggests that changes in the model parameters may allow the progression of lesions over time to be followed and quantified. We are currently investigating this possibility.
Our computational technique is unique in that it relies on a rigorous statistical basis of time series and partial difference equations thereby emphasizing relationships of small groups of pixels over the whole area of interest. The technique will only be useful if future studies show that such relationship information is characteristic of images and that it is medically relevant.
The usefulness of separating scans into non-demented and demented groups in clinical practice remains to be seen. However, we have a tool to make such an investigation possible. Presently, our analysis is done in MATLAB© which is readily incorporated in clinical applications. The basic concepts presented could be applied to most computational imaging systems. Fortunately, these concepts are not complex nor are the model estimations computationally time consuming; the modeling and class separation optimization done for the 97 subjects of the type shown in Figure 3b required 540ms. In a previous paper  we compare our method to 13 image modeling methods and find it 14.3 to 74 times faster and more accurate at classification accuracy than all but one.
The value of the fine structure analysis will depend on the imaging modality and pathology for region of interest studies. The method we have outlined should be informative if regional variations in pixel values underlie a certain pathological process. The type of two-dimensional image could be MRI, CT, SPECT or ultrasound as long as the individual pixels can be stacked appropriately for analysis. More information might be found by limiting the image analysis to specific areas of pathology, such as a tumor mass. If progression of a tumor is reflected by changes in pixel fine structure then it could be valuable in following tumor progression.
This work was supported by grant; DK104393 National Institutes of Health, and Research to Prevent Blindness Foundation awarded to Dr. Mahnaz Shahidi.
We are indebted to Mahnaz Shahidi, PhD for her advice and help in developing our fine structure algorithm and applying it to other medical fields.
- O’ Neill W, Penn R, Werner M, et al. A theory of fine structure image models with an application to detection and classification of dementia. Quant Imaging Med Surg. 2015;5:356-67.
- Box G, Jenkins GM, Reinsel GC. Time Series Analysis. 3rd Edition, Upper Saddle River NJ: Prentice- Hall.1994.
- Fuller W. Introduction to Statistical Time Series. New York: Wiley & Son. 1996.
- Johnston J. Econometric Methods. New York: McGraw-Hill Book Company.1984.
- Magnus J. Linear Structures. Oxford U. Press, New York. 1988.
- Wilks S. Mathematical Statistics. New York: Wiley. 1962.
- Hanley JA, McNeil B. The meaning and use of the area under a ROC curve. Radiology. 1982;43:29-36.
- Massey F. The Kolmogorov-Smirnov test for goodness of fit. J Amer Stat Assoc.1951;46:68-78.
- Rebecca PA. Use of the Jackknife Statistic to Evaluate Result Replicability. J Gen Psychol. 1998;218-228.
- Jain A. Fundamentals of Digital Image Processing. Englewood Cliffs NJ: Prentice-Hall. 1989.