Figure 1.1: Spin of nucleus and magnetic moment along spin axis [18].
Different nuclei have stronger or weaker magnet moments, the stronger moments providing easier MR detection. The hydrogen nucleus, has the strongest magnetic moment and therefore is the most easily detected. This is useful considering that water accounts for approximately fifty-five percent of human body weight, and that the water content rises to ninety percent for some soft tissues. This makes water the greatest contributor to any MR signal, except in tissues where lipids are highly concentrated. Bone contains little water and therefore is almost invisible to MRI.
Under an external magnetic field all the nuclei move, from their random orientations, to align themselves with the magnetic field. For nuclei with a magnetic spin quantum number equal to a half there exist only two orientations of alignment, as shown in Fig 1.2, other nuclei have more than two orientations. These orientations, or spin states, have different energies involved with them. A nucleus with an orientation where the magnetic moment points in the same direction as the external field has parallel alignment, and is the lower energy state. Whereas a nucleus with magnetic moment pointing opposite to the external field has anti-parallel alignment, and is the higher energy state.
Figure 1.2: (a) Nuclei with no preferred orientation in the absence of a magnetic field. (b) Nuclei aligned with external magnetic field, B, in parallel and anti-parallel orientations.
Detection is accomplished through the use of radio frequency (RF) pulses, and their ability to cause resonance, the alternate absorption and dissipation of energy. RF has both electric and magnetic components. In the presence of a continuos magnetic field, the magnetic component of a specific frequency RF pulse can be used to rotate the average direction of the nuclear magnetic moment by 90û from the longitudinal to the transverse plane. This specific frequency is called the Larmor frequency. When the RF is removed, the magnetic moment of the nuclei rotates in the transverse plane about the vector prescribing the constant field, Fig 1.3. This rotating magnetisation, which occurs at the Larmor frequency and has a magnitude which decreases over time, can induce an analogue current in a receiver coil and therefore be examined. The voltage of the induced current is a damped sinusoid decreasing exponentially with time and is termed a free induction decay (FID). The decay rate is characterised by a time constant T2*.
Figure 1.3: Under the influence of an RF pulse, M1 rotates to M2. In the transverse plane M2 rotates around B inducing a current in the receiver coil.
The decay rate, which is due to relaxation of the magnetic moments to the original position, varies depending on the nuclei. Two different mechanisms cause the relaxation to occur, providing two different relaxation times.
T1, the longitudinal relaxation time, describes the transfer of energy between the spin state and the lattice. This is caused by the transfer of energy between spins, in order to restore population equilibrium, giving some energy to the surrounding environment or lattice at the same time. The second of these, T2, the transverse relaxation time, is the loss of magnetisation from the x-y plane of the nuclei. It occurs particularly due to the inhomogenous nature of the constant magnetic field, yielding a range of possible Larmor frequencies and causing an average cancellation of magnetic moments.
The FIDs obtained using a receiver coil are all measures of amplitude variations with time. In order to construct an image, a signal of the amplitude variations with frequency is required. As an FID consists of a damped sinusoid, this would suggest that such a problem is easy. However, the receiver coil actual gets the FIDs from all locations simultaneously, and therefore the received signal actually consists of a number of overlapped FIDs each with different frequency components and phases. In order to create the relevant signal required to construct an image, the 2D Fourier transform is used, along with complex signal analysis.
The resulting images exist in one of three main types, or a mixture of the three, these are T1, T2 and proton density. As has been explained T1 and T2 refer to the longitudinal and transverse relaxation times of nuclei in the tissue. Proton density simply refers to the magnitude of the magnetic moments.
Until the arrival of MRI only axial images of the body had been available, using CT scanning. However MRI provided a method to obtain both sagittal and coronal images, in addition to oblique imaging. This along with the lack of image degradation that plagues CT scans due to the presence of bone, allows an improved method of imaging, especially for soft tissue regions surrounded by bone such as the brain.
While these properties offer great advantages for imaging, it is particularly the ability to provide contrast between different soft tissues, specifically grey and white matter, that makes MRI so important for brain imaging. Certain diseases which were only with difficulty diagnosed without MRI, such as multiple sclerosis and Alzheimer's disease, are more easily diagnosed with it. CT scans are however more useful where there is calcification of the brain. While these show up as areas of decreased signal in MRI, these areas provide no features which are of any diagnostic use.
1.2.1 Alzheimer's Disease
Alzheimer's disease is the commonest cause of dementia in the western world, effecting more than 500,000 people in the United Kingdom alone, some as young as 40 years old. It is an interesting disorder in that is very difficult to diagnose, and at present has no known cause. Patients are only diagnosed as probably having Alzheimer's disease, with confirmation being obtained only after death by an autopsy. The histological signs of the disease consist of grey matter thinning, the presence of amyloid plaques, starch like fibrous cores surrounded by deteriorating neurones, and neurofibrillary tangles.
Two specific forms of Alzheimer's disease exist, familial and sporadic. Approximately 15 percent of all cases are familial, where the cause of the disease is known to be genetic. This Familial Alzheimer's disease was discovered by noticing that people affected by Down's Syndrome acquired both a dementia and amaloid plaques similar to those found in Alzheimer's disease patients. Down's Syndrome patients are known to have an extra chromosome 21, and when examined it was found that some Familial Alzheimer's disease patients had defects on chromosome 21, 14 or 19, work is still in progress to identify other mutations. Familial Alzheimer's is of great interest since it provides a number of possible patients who will, with all probability, develop the disease. This means that MR images taken of these patients in their early 40s will provide images indicating earliest onset of the disease. Also sufferers of Familial Alzheimer's disease are known to develop the disease at an early age, avoiding such things as strokes, hypertension, vascular disease etc, which can all cause dementia.
It is hoped that MRI will provide the possibility of being able to exactly diagnose the presence of Alzheimer's disease before death. At present no early signs of the dementia are visible using any method of medical imaging. Instead, the probable diagnosis is made when physical symptoms, such as short-term memory loss, are severe enough to be recognised by the kin of the patient, which unfortunately occurs some time after the first physical signs manifest themselves. MRI of patients with the later stages of Alzheimer's disease shows a decrease in brain size, especially around the temporal lobes and the area of the hippocampus, as can be seen in Fig 1.4, in addition to sulcal broadening, enlarged ventricles and the reduction of the gyral size.
Figure 1.4: Cortical MR Images of (right) a patient with well developed Alzheimer's disease, note the atrophy that has occurred, and (left) a control subject.
Two problems therefore exist. Firstly, to diagnose Alzheimer's disease before the manifestation of any physical symptoms, and secondly, to discover the cause of the degeneration of the brain. Detection of the earliest changes in the brain which cause Alzheimer's disease will allow any drugs, experimental or not, available to stop the atrophy of the brain, to be administered at the earliest time. This is necessary, since once destroyed, neurones in the brain cannot be replaced. While the deteriorating condition of a person in the latter stages of the disease might be stabilised, using drugs, they will live for an extra twenty or so years with the disabling effects of the disease, rather than dying a few years later from the disease itself. The solution is thought to be image analysis of the MR images, obtained from patients with different stages of Alzheimer's disease.
Different methods of image analysis and calculation for predicting the onset of Alzheimer's disease have been tried. Volumetric measurements of the hippocampi, the cerebrospinal fluid and some of the ventricles have been used. J. Seab et al [23] first showed that hippocampal volume decreased significantly in patients with Alzheimer's disease. However volumetric methods have suffered from the problem that actual brain size and characteristics vary quite significantly between different people of the same age and overlaps occur in the hippocampal volume between Alzheimer's disease suffers and controls [15]. Texture and morphological analysis of MR images may offer the best solution. Both these methods are well suited to dealing with the initially tiny changes that occur in the brain.
1.2.2 MRI Analysis
Usually MR images are examined by the clinician by eye. While this is often sufficient due to the excellent ability of the eye at detecting small changes in image intensity, shape and texture, in cases where a tumour is just beginning to form, or where the exact size and nature of a lesion is required, this is not enough. In situations such as these, MR images need to be analysed and enhanced, so as to enhance the region which is of interest to the clinician, whether it be a tumour or some haematoma.
In the area of tissue characterisation and diagnosis, MRI is much more successful than CT, for example MRI is much better at separating the oedema from the tumour, allowing the surgeon, extracting the tumour, to remove the minimum amount of brain tissue. Variations between similar lesions in patients of different ages or in patients with lesions at different stages of development, provide different characteristics. This problem is further enhanced by the lack of a standard in MRI. At present similar images of a subject, obtained using different MRI machines, have different characteristics and therefore two images obtained under different situations cannot be compared.
Texture analysis examines and quantifies the small regular variations in intensity which make up a texture. Rather than simply examining intensity, the relation of the intensity of a point is examined with respect to its surroundings. While white and grey matter in MR images may vary in intensity dependent on the machine, person and magnetic field used, their texture will always remain the same. Further more, similar tumours would also have a similar texture defined by the pathology of the tissue causing the tumour.
The effect of noise in MRI is of particular importance when dealing with image analysis, specifically texture analysis. The signal to noise ratio (SNR) of any image depends inversely upon resolution. The equation 1.1 is a simplified formula of the SNR for a pixel [10].
(1.1)
While increasing the size of the image pixels, dx and dy, increases the SNR for each pixel, it also has the effect of increasing the size of the smallest object that can be perceived, and therefore decreasing the resolution. Increasing the slice thickness, dsl, has a similar effect. It improves the SNR, but by increasing the slice thickness, blurring of the image occurs, because the tissue in the pixel is Ôaveraged' to create the overall intensity, these are called partial volume effects. Increasing Nx and Ny, the number of samples in the x and the y direction in the image, improves the SNR partly by the averaging that occurs, and partly by reducing the bandwidth per pixel.
The dominant source of noise in MRI is Johnson, or thermal noise, caused by resistance in the process. An increase in the bandwidth, BW, of the signal causes a decrease in the SNR, equivalent to the root of the bandwidth, due to the effects of Johnson noise. In practice the magnitude of the bandwidth is constrained by the sampling time. The longer the sampling time the smaller the bandwidth required, but the more susceptible the process is to movement of both the subject being scanned and the internal movement of both blood and organs. The SNR increases with increasing operating frequency, w, where K and C are constants corresponding to resistance in the process. This occurs because when under the influence of the external magnetic field, more nuclei move into the anti-parallel alignment at higher frequencies than at lower frequencies. There is a maximum frequency at which surface conduction starts to occur on the skin, effectively screening the inside of the body from RF.
From a view of texture analysis, especially texture mapping, any method used must deal with a small enough segment of the image so that resolution is not lost. However if there is a large amount of noise in the image, it can affect the texture mapping when the segment size is very small. In such a case the texture analysis could be examining the texture of the noise rather than that of the image, causing errors.
There is no exact or even well defined description of texture. When one observes sand on the beach, or the bark of a tree, the homogenous visual pattern one observes is texture. It not only corresponds to variations in intensity, but also comes from the variations of the surface characteristics of the item, which can themselves be felt. Brodatz [5], was the first person to provide photographic examples of natural texture, ranging from grass to piles of coins. His images are often used as test images for texture analysis, because of their advantage over computer generated textures. A computer texture generated by an algorithm has the disadvantage in that what ever method is used to create it, is also the texture method that has the highest probability of identifying it.
Much of the interest in texture stems from the ease with which the eye can perceive it. Many textural models based on the human visual system [24], [25], have been implemented, with varying degrees of success. Neural networks have also been used [6], in the hope that texture analysis methods can be improved by simulating the method by which the brain, through the use of the eyes, detects texture.
Texture analysis has been used extensively, from practical uses, in robot vision on an assembly line [20], to classification of terrain from aerial images by Haralick et al [14]. The two main uses of texture analysis seem to be, firstly the characterisation of unknown textures corresponding to different objects, examples being different tissue types and various types of woodland. Secondly, the identification of different areas of texture for enhancement, an example being the enhancement of tumours for clinical study.
Texture can be divided into two very distinct types, namely statistical and structural (Fig 2.1). For statistical definitions, texture can be defined as the arrangement or spatial distribution of intensity variations in an image according to some underlying probabilistic model. Structural texture is the placement or spatial distribution of a set of primitives in an images based on some predefined placement rules [26].
Figure 2.1: Examples of computer generated texture.
Examples of (left) statistical texture and, (right) structural texture
2.1 Statistical Texture Analysis Methods
The importance of human visual perception means that much effort has been made to identify the characteristics of texture which are recognisable by it. Julesz [16], [17], was the first to attempt to determine this. He found that human observers could distinguish between textures with different first order statistics, and also between textures with the same first order statistic, but with different second order statistics (Fig 2.2). However observers were unable to perceive any differences, where there was a difference in the third or higher statistics. First order statistics deal with the intensity of a single pixel, while second order statistics compare the differences between the intensities of two pixels at different points and so on.
Figure 2.2: Textures with different (left) second and (right) third order statistics [16].
Much criticism has been made of Julesz conclusions [26]. However a great amount of the research into statistical methods continues to concentrate on the first and second order statistics, because of Julesz's conjecture that the majority of textural features is contained therein.
2.1.1 Histogram Based Methods
The histogram of an image, a first order statistical method, provides the simplest method for describing texture. The histogram of an image is p(zi), for i=0, 1, 2, ..., n, where n is the total number of discrete grey levels (Fig 2.3) and zi is a random variable denoting discrete image intensity. Each point in the histogram is the occurrence probability of that intensity level i.
Figure 2.3: A simple grey scale image and its histogram.
Different parameters can be worked out from the histogram, three examples being the mean, equation 2.1, and the moments, equation 2.2, and the entropy, equation 2.3. The second, third and fourth moments are know as the standard deviation, skew and kurtosis respectively.
(2.1)
(2.2)
(2.3)
The mean takes the average level of intensity of the image or texture being examined, and therefore can be used to blur images by taking the mean in small windows around each pixel. The standard deviation provides useful information about the variation of intensity through the image. The skew is zero if the histogram is symmetrical about the mean, and is otherwise either positive or negative depending on whether it has been skewed right or left of the mean. Kurtosis is a measure of the flatness of the histogram, and the entropy, a measure of its uniformity.
Histogram techniques suffer from their reliance on the intensity of a single pixel. If the images contain the same original texture but the mean and standard deviation have been shifted for each, due to a different scanning method for example, histogram based texture methods will obviously identify this one texture as many different textures. For this reason images are often normalised to have the same mean and standard deviation.
2.1.2 Co-Occurrence Matrix
Texture parameters computed using only the histogram of an image are susceptible to the limitation that they provide no information regarding the relative position of pixels, both of the same and different intensities. The co-occurrence matrix is a Ôtwo dimensional histogram' which overcomes this limitation, and is therefore a second order statistical method.
The co-occurrence matrix is an estimate of the second order joint probability, p(i,j), of two pixels, a distance d apart along a given direction q, having the same value (co-occurring). An example is given in equation 2.4. Two forms of co-occurrence matrix exist, one symmetric, where pairs separated by d and -d are counted, and other not symmetric where only separation by a distance d is counted.
(2.4)
The matrix is square, with the width and height being set by the number of discrete grey levels, or intensities, in the image to be examined. Due to the intensive nature of the computation often only the angles of 0û, 45û, 90û and 135û are measured with a distance, d, of 1 or 2 pixels as suggested by Haralick et al [14]. The value of d used is of moderate importance. The classification of fine texture usually requires small values of d, whereas coarse textures require large values of d. As many texture features are independent of any rotation, this can often be reduced to examine a single angle. Further increases in speed can be made by quantizing the image to a fewer levels of intensity, with some loss of textural information.
Haralick et al [14] defined a number of descriptors based on the co-occurrence matrix, some of those are given in equations 2.5 to 2.8, where mx, my and sx, sy are the mean and standard deviations of the row and column sums of the matrix.
(2.5)
(2.6)
(2.7)
(2.8)
Good results have been obtained using these methods [3], [4], on LANDSAT and artificial images. Additional uses of the co-occurrence matrix include edge detection and maximisation.
2.1.3 Grey Level Difference Method
The grey level difference method (GLDM) is a subset of the co-occurrence matrix. The resultant array, of size equal to the number of discrete grey levels in the image, indicates the probability, p(k), of there being a set of two points, a distance d apart in a direction q, with a difference between their intensity equal to k, the position in the grey level array. The array can be obtained from the co-occurrence matrix using equation 2.9.
(2.9)
There is a significant spread in p(k) where the texture is fine, if a coarse texture is being examined, the probabilities are much greater near to k=0. Four features based on the GLDM are shown in equations 2.10 to 2.13.
(2.10)
(2.11)
(2.12)
(2.13)
In order to allow rotationally dependent textures to be characterised with ease, at the cost of its detection of directionality of textures, a variation of the GLDM exists [1], called the neighbourhood grey-tone difference matrix. This method operates by comparing the grey level of a pixel with the mean of the grey levels of a square neighbourhood around that centre pixel.
2.1.4 Run Length Matrix
The run length matrix, p(i,j), is a method of statistical analysis, that represents the frequency that j points with a grey level i continue in the direction q [26]. The i dimension of the matrix corresponds to the grey level and has a length equal to the maximum grey level, n, while the j corresponds to the run length and has length equal to the maximum run length, l. As with the co-occurrence matrix, q = 0û, 45û, 90û and 135û offer the greatest interest, and once again there is rotation dependence. Five features can be calculated from the run length matrix (equations 2.14 to 2.18). A is the area of the image.
(2.14)
(2.15)
(2.16)
(2.17)
(2.18)
The image is often quantized to fewer levels to cut down on the dimension length and therefore the processing time. This however does reduce texture quality. The run lengths are very large for macrotextures, especially structural textures, but can be quite small for fine textures. The nonuniformity features are small if the grey levels or the run lengths are similar throughout the matrix, while the long run length is large if there is high intensity clustering in the texture.
2.2 Spectral-Based Statistical Methods
2.2.1 Fourier Power Spectrum
Texture analysis in the frequency domain can, by its quasi-periodic nature, provide definite results. The Fourier transform can be used to provide images in the frequency domain, although the Fast Fourier Transform (FFT) is quicker where images have a size which is the power of two. In order to characterise some unknown texture, the ring or wedge sum of the power spectrum (Fig 2.4 and equation 2.20) can be calculated.
(2.20)
Ring sums are used to calculate the magnitude of the image at a particular range of frequencies, whereas wedge sums are useful in detecting directionality, as texture which is rotationally dependent gives lines of interest in specific directions.
Figure 2.4: Fourier power spectrum: (a) Ring filter; (b) Wedge filter [26].
If the texture being analysed is vaguely sinusoidal in nature, the magnitude of frequency corresponding to that sinusoid provides a direct measure of its area. This stems from the Fourier transform representing the image as a two dimensional array of varying sinusoids. Generally fine texture is scattered throughout the spectrum, whereas coarse texture is concentrated at the lower frequencies.
2.2.2 Autocorrelation Function
The autocorrelation function is similar to the power spectrum. It is obtained by doing the reverse Fourier Transform on the power spectrum of an image. Until the arrival of the FFT it was generally easier to use the equation 2.21, but this is very computational intensive.
(2.21)
The effect of the autocorrelation function is to determine the coarseness of a texture. The coarser the texture, the larger the texture primitives and so the autocorrelation function drops to zero more quickly. A fine texture will give a smaller gradient. If the there is any rise in the function, this indicates that the texture involved is periodic. The value of 1/e is usually used to characterise the texture [10], [13].
2.3 Model-Based Statistical Methods
Model-based analysis methods exist because of their use in creating artificial texture. This means, however, that when tested an analysis method based on one model, works best with textures created using that same model. As well as being computationally intensive, their general structure is usually complicated and for that reason they are not included here.
2.4.1 Edge and Spot Detection
Edges occur and are defined as occurring where the intensity of an image varies abruptly. The edges in an image can be obtained either through convolution with an edge detection matrix (Fig 2.5), or by using the strength and direction equations given in equation 2.22.
(2.22)
where 
The edge detection method above allows the derivation of texture properties. The coarseness of the texture is measured by the density of the edge elements, a coarse texture having a low edge density. Texture contrast can be measured as the mean of the edge strengths. A histogram of both the edge strength and direction can be calculated, allowing textural features, similar to first order histogram features, to be obtained. Randomness of the texture is the edge strength histogram entropy, while a rough measure of texture directivity is the entropy of the edge directivity histogram.
Figure 2.5: A texture image (left) with its edges detected (right).
The co-occurrence matrix of the edge strength and direction can be extracted allowing the angular second moment, correlation, etc to be calculated. Edge analysis is particularly successful on images containing large homogenous areas of texture.
In a similar way spot detection measures texture as a function of the amount of spots of different sizes in the image. In practice a minimum and a maximum spot size are specified, and any spots outside these ranges are ignored. Spot analysis works well where the texture has a repetitive pattern of high intensity blobs, for example sandpaper.
2.4.2 The Peak and Valley (Extrema) Method
The local maxima and minima of certain sizes can be used to analyse texture. The method proposed by Mitchell et al [26] first used smoothing to remove any minor extrema and small variations caused by noise. The value at which the smoothing cuts off acts to define how many of the extrema are maintained, and therefore is of prime importance. This method proceeds to work along the lines of the image detecting changes of the gradient sign, when pixel intensity peaks or troughs. The number of extrema of different sizes characterise the texture, and the cut-off value can be varied to provide different features.
2.4.3 Random Walks
The random walk is based on the statistical movement of a particle. The random walk pointer can move in one of the four horizontal and vertical directions, and the direction of walk is determined by the pixels in these directions. If the pixels are all of the same intensity the probability of moving to any one of them is 0.25, however if they are all different the probabilities vary favouring a move to the pixel with smallest difference in intensity. If a small window is move around the image, the probable number of steps before hitting the window boundary can be calculated and used as a textural feature.
2.4.4 Moments
The moments of an image provide an analytical texture method that can be used to characterise both statistical and structural texture. For a 2D digital image f(x,y), the moment of order (p+q) , as defined in equation 2.23, exists so long as f(x,y) is continuous and non zero only in a finite part of the image.
(2.23)
The use of this simple moment is limited since its size depends greatly on the size of the image, making it sensitive to variations in scale. Normalised central moments, equation 2.25, a variation of central moments, equation 2.24, which are calculated from the moments about the mean, have the advantage of being less susceptible to variations in size.
(2.24)
(2.25)
Neither of these moment measures are invariant to rotation or a large degree of scaling. For this reason a set of seven invariant moments exist (equations 2.26) [12], which are, to a large degree, invariant to rotation, translation and scale change.
(2.26)
This texture analysis method is of particular use where the orientation and size of different images varies, one example is the texture characterisation of meteorological images.
While structural approaches to texture analysis are important, their uses are limited in the analysis of MR images. Biomedical textures tend to be statistical in nature because of the resolution of most imaging methods, rather than because they actually are. Cells in the human body are arranged in a regular way, demonstrating a very strong structural pattern. Statistical analysis of MR images is more likely to succeed since firstly, there are many more statistical analysis methods around in general, and more are suited to this particular application than structural methods. Secondly, the brain as it appears in MR images does not appear to have any structural texture contained within.
Structural texture analysis describes texture in terms of primitives and placement rules. These texture primitives can be used to create more complex patterns through repetition and variations in orientation and size. Structural approaches usually attempt to recognise the primitives in an unknown texture, and then to work out the placement rules for the texture and so characterise it. Training images are often used, so that an analysis method only needs to compare primitives and placement rules to characterise the texture.
Texture mapping differs from texture analysis in that the output is an image, preferably of the same size and with the same spatial characteristics. One method of doing this is to divide the image to be analysed into a number of square windows and to set the value of each pixel in the window to that of a calculated textural feature. This method, however, gives a resultant image which has a very low resolution. An improved method takes a window, or neighbourhood, of m by n pixels, where m and n are usually equal and odd, and calculates a textural feature for each point in the image and the defined neighbourhood around it (Fig 2.6). This method allows the result for each pixel to be summed into a texture mapped image (Fig 2.7) and for the purposes of this report is called the neighbourhood texture mapping method (NTMM).
Figure 2.6: Neighbourhood texture mapping method
The size of the neighbourhood is of prime importance. If a small neighbourhood is used the resolution of the mapped image will be high. However, if the texture in the image is too coarse, then the measurement of a textural feature will be incorrect for that texture. In addition if the neighbourhood is too small then variations in the image due to noise will be texture mapped rather than the texture variations of the underlying image itself. This is of particular concern in MRI, where the majority of the noise is Johnson noise (Section 1.3), as this provides a noise texture which exists in the whole image. If the neighbourhood is large, the resolution of the mapped image is much lower and there is a greater chance of trying to define more than one texture in the image as a single texture in the mapped image.
(a) (b) (c)
Figure 2.7: Texture mapped standard deviation of (a) with neighbourhood sizes of: (b) 5; (c) 17.
The optimum size of the neighbourhood is equal to the size of the Ôstatistical primitive' making up the image. This is easy to determine if the texture is recognised, however in MR images of the brain the number of textures and the sizes of their primitives are unknown and therefore experimentation with different neighbourhood sizes is required.
The Medical Image Display and Analysis System (MIDAS) is a interface that allows the user to display and manipulate images, specifically medical images such as MR images. MIDAS runs in X-Windows on a UNIX based computer and provides the facilities for loading any number of images, displaying up to six full size 256 by 256 images on the screen at any one time, and processing and saving images. While these are usually grey scale images, there is a colour palette function which allows the colour map to be varied from grey scale to the full colour spectrum, a function that is useful where the grey scale images are of a particularly low intensity.
MIDAS includes a number of tool boxes. These tool boxes contain image processing functions for manipulating the images in MIDAS and then displaying the results (Fig 3.1). Developers can write their own tool boxes, in the C programming language, which use a few basic commands to obtain images from the MIDAS front end, display them and to even plot graphs [11]. Tool box development in MIDAS enables the user to ignore the problems of loading, saving and displaying images and concentrate on the actual development of the tools.
Figure 3.1: The MIDAS window showing the tool box menu on the left.
3.2 Proposed Texture Mapping Routines
All texture mapping routines were programmed in C as tool box additions for MIDAS. These routines are all available in Appendix A, listed in the order of textural features.
To confirm that the computer method was correct, each mapping routine was tested on a 256 by 256 image, containing a white box of size 128 by 128 (Fig 2.7a). The expected texture value for different points in the images were compared with a value calculated on paper. The five points chosen were, one on the vertical line of the white box, the horizontal, top left hand corner, top right hand corner, and in the background outside the box
3.2.1 Requirements for Mapping Methods of Brain MR Images
In order to evaluate texture maps using different routines it was first necessary to decide the textural features that the brain would most likely exhibit, and secondly what features may be able to delineate between a normal brain and a brain with Alzheimer's' disease. It was decided that the neighbourhood texture mapping method (NTMM) would provide the best solution as to a mapping method. The list of requirements for any texture feature are as follows:
1. Must be a low order statistic, because that is where it is suggested the majority of statistical textural information is [16].
2. Needs to be fast and therefore not computationally intense, since the time available for running the algorithms is firstly limited and secondly could be used in a more productive fashion.
3. Must be invariant or have the possibility to be made invariant to changes in the spread of grey level, caused by the use of difference MRI machines and different patients.
4. Requantisation to fewer levels must not to be necessary.
5. Results should be directly comparable to the image in spatial features.
The first requirement suggests that all the first and second order statistical methods should be implemented, including spectral methods such the power spectrum. However the limitations of the FFT, which would have to be used, requires a window size of the power of two. This makes where to place the central pixel a problem, as well as the feasibility of taking the FFT of a neighbourhood containing a minimum of 4 to 16 pixels.
The time limitations of this project mean that high speed methods are required which work in a matter of minutes and not hours. The very nature of the NTMM, the calculation of a texture feature in the neighbourhood of a pixel for all pixels in the image, is computationally intensive. The internal processing required per neighbourhood makes model-based analysis methods, power spectrum routines, as well as the edge detection methods, unfeasible. The size of the co-occurrence matrix would suggest a similar problem. However as well as Argenti et al [3], [4], offering a solution for fast calculation of the co-occurrence matrix, the author has some ideas which proceed to work.
Moments and the autocorrelation function are invariant to rotation, scaling and translation while the histogram and co-occurrence methods, including difference statistics, are partly invariant, considering the intensity distribution of the brain.
Requantization to fewer grey levels is often necessary for the co-occurrence matrix and histogram. However a method is proposed which makes the co-occurrence matrix and, to a lesser extent, the histogram independent of the number of grey levels.
Finally, all the methods map to an image which can be directly compared, by overlay for example, with the original image.
3.2.2 Histogram Routines
Histogram routines based on the NTMM included the mean, standard deviation, skew and kurtosis (flatness). While all the programs had similar layout, they were programmed as separate files to improve the ease in which they could be debugged. At the centre of each program were two loops which moved a pointer through the image, and four variables, which specified the neighbourhood limits as a function of the image geometry.
It is suggested that Requantization of the MR images to fewer levels would speed up the image processing [14], however with the relatively small variations in the grey and white matter this would almost certainly remove important textural features. In order to obtain the texture value for a neighbourhood, its histogram had to be calculated. The histogram was an array of floating point variables of length equal to the number of grey levels. This was done by looping through the neighbourhood, incrementing the values of the histogram corresponding to the grey levels of the pixels, divided by the area of the neighbourhood, resulting in a probability function. To obtain the texture value, for the centre pixel in the neighbourhood, a second loop went through the histogram as suggested by the equations 2.1 and 2.2.
However it was found that for a 256 by 256 image with 256 grey levels and a neighbourhood size of 11, at least 256x256x(256+11x11) = 24.7 million loops were needed. This was too slow, and, therefore a method was divised where the loop number was independent of the number of grey levels. The grey level values of the pixels in the neighbourhood in the image were used as pointers to the corresponding grey level value in the histogram (Fig 3.2). Once the pdf corresponding to a grey level value had been read from the histogram, it was necessary to reset it to zero, so as to avoid it being read again.
Figure 3.2: Image and histogram showing how the neighbourhood could be used to reference the histogram.
This method has the effect of making the number of loops completely dependent upon the window size, but since this was partly the case anyway it is the optimum method available. For the previous example this was found to reduce the number of loops from 24.7 million to 16 million, 256x256x2x11x11, however this will vary with the neighbourhood size.
MIDAS was found to contain a built in histogram function, however with MR images of the brain the percentage of low grey levels was such that the pdf of the region of interest, the brain itself, was difficult to examine. For this reason a second histogram function was written that in addition to allowing the user to specify a bottom cut-off and so remove any black values from the histogram, allowed the histogram values to be saved to a file.
The five histogram NTMM routines, mean, standard deviation, skew, kurtosis and nth moment require the user to input an odd neighbourhood (window) size. They were tested using the white on black box image in fig 2.7a, and found to be working as their algorithms suggested (Also Fig 2.7).
3.2.3 Co-occurrence Routines
The implementation of the co-occurrence routines, angular second moment, element difference moment, inverse element difference moment, and correlation, was similar to the implementation for the histograms.
The co-occurrence matrix was initialised as a two dimensional floating point array where the size of both dimensions was equal to the number of grey levels. Inside the main loop a second loop worked through the neighbourhood calculating its co-occurrence matrix, after which a third loop calculated the required textural feature from the matrix. The neighbourhood limits had to be checked and changed for each pixel. This occurred because the neighbourhood size was not its specified size, for example 7 by 7, it was dependent upon the distance and angle variables and was always less than the specified size. In order to avoid having to requantize the image to fewer levels, it was necessary to make the program avoid doing a total of 4300 million, 256x256x(11x11+256x256), loops per image. As this relied on working through the whole co-occurrence matrix for each pixel in the image, a solution similar to that for the histogram was devised.
The grey levels of the pixels in the neighbourhood were used to set and obtain the probability values from the co-occurrence matrix. In order to optimise the basic routine further, it was necessary to remove the in loop calculation of the pixels in the direction and at the angle required from the neighbourhood pixels. To this end a second image matrix was used that set the value of a pixel in the image to the value of the Ôco-occurring' pixel at a distance and angle set by the user-specified parameters for the co-occurrence matrix. These improvements made the routines invariant to changes in the number of grey levels and only dependent on the neighbourhood size, the result being an approximate speed increase of 300 times.
The four routines were tested on the square test image, fig 2.7a, and the results, fig 3.3, were verified by hand. The user was required to enter the distance d and two offsets, dx and dy, corresponding to the angle q, as well as the neighbourhood size. For the angle q, +dx corresponds to right across the image and +dy to down the image, meaning that for this project an angle of 45û is 45û downwards from the horizontal.
(a) (b)
(c) (d)
Figure 3.3: Co-occurrence test results of fig 2.7a with d=1, q=45û and a neighbourhood size of 7: (a) Angular second moment; (b) 1st element difference moment (e.d.m.); (b) 1st inverse e.d.m.; (c) Correlation with q=0û.
Due to the large area of black surrounding an MR Image, the co-occurrence matrix has an equally large spike at its zero, zero point. This spike is undesirable since some of the routines set the centre pixel as a value based only on the contents of the co-occurrence matrix and not the grey level associated with it. In order to reduce the magnitudes involved a second set of routines was written that allowed a bottom grey level cut-off to be specified and so remove the problem of the spike.
3.2.4 Difference Statistics
These routines, angular second moment, contrast and mean, are all based on the co-occurrence routine. In addition to the co-occurrence matrix, a difference statistics array was necessary with length equal to the total number of grey levels.
Once the co-occurrence matrix has been obtained the difference statistics array is calculated using equation 2.9. This can then be used to calculated the different texture features. As with the co-occurrence routines optimisation was carried out using a second image matrix and by using the neighbourhood values as pointers to the positions in both the co-occurrence and difference arrays. In addition a second set of routines were written to allow the existence of a bottom cut-off.
The routines were tested and the results found to coincide with those calculated on paper (Fig 3.4).
(a) (b) (c)
Figure 3.4: Difference Statistics results of fig 2.7a with d=1, q=45û and a neighbourhood size of 7: (a) Angular second moment; (b) Contrast; (c) Mean.
3.2.5 Moments
The moment routines, original, central, normalized central, and invariant, proved to be simpler than those for the histogram or the co-occurrence matrix. Equations 2.23 to 2.26 are as optimised as necessary for the purposes of computer implementation and in addition to this there are no extra arrays required. The input variables required are the neighbourhood size, and the movement orders p and q.
The central moments routine incorporated the original moments routine, while the normalized central moments routine incorporated the central moments routine. In order to create a suitable method for calculating the seven invariant moments a function was added which returned the value of a normalized central moment for a specific p and q, and for a specific pixel and its neighbourhood. Due to the old version of C being used it turned out to be quite difficult to create this function, however in the end it was incorporated as part of the first invariant moments routine. Each of these routines consist of simple mathematical calculations carried out upon the results of the normalized central moment of each pixel.
It was found to be very difficult to calculate the invariant moments on paper, in order to get round this, simple mathematic routines were written that allowed the image calculation of the invariant moments using the images obtained from the correct normalized central moments. The other routines were tested on fig 2.7a and the results, fig 3.5, were confirmed on paper.
(a) (b) (c)
Figure 3.5: Moment results of fig 2.7a with p=q=1and a neighbourhood size of 7:
(a) Original moment; (b) Central moment; (c) Normalized central moment.
3.2.6 Autocorrelation Function
The autocorrelation function routine is not a NTMM, while it would be desirable to make it so, it would be very computationally intensive. Originally a simple routine was written that would calculate the autocorrelation function of an image using the algorithm shown in equation 2.21. However it turned out that this algorithm took in the region of an hour to run, and therefore an improved method was necessary.
MIDAS includes a forward and reverse FFT built in. It was found that the end result of taking the FFT of an image, the power spectrum, and then the reverse FFT, was the autocorrelation function, a procedure that took one minute. The routine was tested, but it proved to be difficult to ensure that the result was correct, however after reducing the test image size the routine was confirmed (Fig 3.6).
Figure 3.6: Autocorrelation of fig 2.7a with a neighbourhood size of 7.
3.3 Tools Developed and Assumptions
3.3.1 MR Images to be Used
For purposes of evaluating the texture mapping routines a number of MR images had to be selected. Each patient had a volume containing 124 images, where each image was a coronal brain slice. The images were scanned using a 1.5T GE Signa scanner with a TR=35ms, TE=5ms, the number of excitations=1, and a flip angle of 35û. The volume images were obtained using the spoiled gradient echo technique (SPGR). Images were contiguous, 1.5mm thick, and were initially 256 by 128 in size, interpolated to 256 by 256. As a maximum of one or two images per volume could be selected, it was necessary to decide what would be the most effective image to map. After experimentation an image slice approximately half way through, in line with the lobes of ears, was decided upon. This slice showed the largest area of temporal lobe and the hippocampi also appeared prominently. This is desirable since in Alzheimer's disease atrophy first starts to occur in these areas.
Since the purpose of this project is to compare the ability of different textural features in delineating between the brains of patients with and without Alzheimer's disease, images of patients in the different stages of Alzheimer's disease were required, including some controls. In the end a main group of six patients were selected (Fig 3.7), consisting of, two age-matched controls (C1&C2), one at risk (R1), one mild (M1) and two severe (S1&S2). The coronal slice corresponding as near as possible to the ear lobes was used for each patient. The additional images of six patients (C3, C4, R2, R3, M2, and S3) were used where the texture mapping method looked hopeful. The at risk patients are cases of familial Alzheimer's disease where the patient is roughly at the age at which symptoms are expected to appear.
Figure 3.7: MR evaluation images.
3.3.2 Tools for Image Manipulation
The C listings of the software tools for image manipulation can be seen in Appendix B.
The MR evaluation images cannot be compared with one another, since they are of different patients with a different range of grey levels corresponding to the same soft tissues. It was therefore necessary to find some normalisation method that would allow the textured mapped results of the images to be compared with one another.
One method often used is to normalise the images so that they have the same mean and standard deviation, because these two statistics are affected by the input variations. However MR images do not have a set number of grey scales and so a normalisation routine of this type would not work. Instead it was decided to normalise the images using points corresponding to two materials which were known to remain the same whoever the patient and whatever the MRI method.
The background content of an image remains the same for all images, it usually has a grey level value varying between 0 and 5. The content of cerebrospinal fluid (CSF) is known to stay constant in people, except those who have undertaken very rigorous exercise, for example boxers. The CSF was found to have a minimum grey level that varied between 30 and 50 for different images. In the normalized image, the CSF and background were set to two arbitrary values, 40 and 0 respectively, and the image was scaled to these values. The normalisation routine worked by taking and averaging 9 pixels, around a user specified point, to obtain the image values of CSF and background. These values were then used to scale the image so that the grey level of the CSF became 40 and the grey level of the background became 0.
The MR images of the brain contain a great deal of useless data. It would be desirable to have images only containing the brain and not the fat surrounding the skull or the material under the brain. Three routines were written that each removed some of the undesirable features and minimise the head scan to a brain scan. The first routine was based on two horizontal cut-off, specified by the user, which removed the top of skull fat and the material beneath the brain, and a function which attempted to locate the CSF points in each horizontal line using intensity thresh-holding with values specified by the user (Fig 3.8). While this routine worked well, it often left bits of skull fat on the corners of the brain and sometimes a little bit at the side.
(a) (b)
Figure 3.8: (a) Original image with horizontal cut-offs indicated in white; (b) Minimised image, note the remaining skull fat at the corners of the brain.
Two routines, a diagonal and a vertical cut-off, enabled the user to specify four points corresponding to a 45û diagonal cut-off and two vertical lines for left and right cut-offs (Fig 3.9).
Figure 3.9: The final minimised MR evaluation images from fig 3.7, after diagonals and horizontals have been removed, with almost no excess tissue. The images in fig 3.7 were also CSF normalized to give these images.
Once the texture mapped image is obtained, it would be good if the image could be directly compared with the original using an overlay routine. This can simply be done by adding two images, with a similar grey scale range, together. In order to scale an image two routines were written that multiplied an image by a user specified number, and added a specified number to an image. To overlay two images a multiply images routine and an add images routine were written.
In order to keep the continuity of the text, the resulting images, figures 4.1 to 4.18, have all been put at the end of this chapter.
The simple histograms, graph 4.1, were calculated for the evaluation images in fig 3.9. The peak value, at an approximate grey level of 180, corresponded to white matter, and a small maximum at a grey level of about 130 corresponded to grey matter. It was noted that the spread between grey and white matter peaks seemed to decrease in the Alzheimer's patients S1 and S2 as compared to the controls C1 and C2. However these results cannot be trusted since are cross-section of the brain and not the histogram of the brain as a whole, where the predominance of grey matter at the ends of the brain would give much higher grey matter levels. It has been shown [7], [8], that there is a definite variation between the grey and white matter peaks in the histograms of the entire brain.
Graph 4.1: Histograms of C1, C2, S1 and S2, note the larger spread for the controls.
The mean NTMM was not examined since it simply blurs the image by setting every pixel to the neighbourhood average.
The standard deviation texture map was calculated on C1 for neighbourhood sizes from 3 to 47 (Fig 4.1) in an effort to obtain the optimum neighbourhood size for the map. From the images it is clear that a neighbourhood of 5 is too small while 25 is too large. A neighbourhood size of 17 was decided upon and the standard deviation was mapped for the evaluation images (Fig 4.2). It can be seen from the results that there is a definite difference between the controls and the severe cases. The high intensity regions at the edges of the brain correspond to large grey level variations. This suggests that the interface between grey and white matter and the thin layer of external CSF, at the edge of the brain, varies less in intensity, or has a smaller gradient, for Alzheimer's patients than for the controls. While both M1 and R1 show a similar decrease in intensity variation, the difference between M1, R1 and C2 is not large enough to drawn any conclusions. Further mapping of M2 and R2 showed similar slight variations, but quantitative analysis is necessary to make sure. There is a noticeable high intensity region around the temporal lobe area of all the images, this is because in the original images these areas had the highest grey level values.
A neighbourhood value of 17 was also found to be optimum for the histogram skew (Fig 4.3). The resulting texture maps (Fig 4.4) for the evaluation images show that the interface between the grey and white matter, the inner of the two bands of high intensity, is being enhanced as well as the CSF-matter interface. The average intensity and size of this inner band can be seen to be much larger for the controls, especially C1, than for all the other images. This is a much better indication of changes between controls and the affected patients, since the intensity variation of the two bands are more clear cut, and more importantly the variation across the grey matter, white matter interface is being shown. The same pattern was found in the additional images, and if this is an indication of the physical effects of Alzheimer's disease on the brain, then R1 and possibly R2 are in the early stages of Alzheimer's disease.
When different neighbourhoods were tried using the kurtosis (flatness) map, the results for neighbourhoods greater than 15 seemed to low in resolution to be useful (Fig 4.5). The single band of intensity around the brain CSF interface is similar to that found for the standard deviation except with a wider band. The only different feature of interest (Fig 4.6) was the increase in average grey level intensity inside the band for S1 and S2 as compared to the other images. However this can be put down to variations cause by an increased interface due to widening of the sulci and the narrowing of the gyri.
4.2 Co-occurrence Matrix and Difference Statistics
For each of the basic routines, the difference between the original routine and the cut-off routine was examined using different images with different neighbourhood sizes. In addition to this, different distances and angles were investigated, so as to see if any improvements might result from adjusting these variables. The angle of 45û is of particular interest since visual observation of the images (Fig 3.9) shows vague line-like textural features coming from the centre of the brain outwards, mostly, at the diagonals corresponding to this angle. It was also noted that for some of the images the variations in intensity (grey level) could be enhanced by the use of a full colour spectrum, where red corresponded to high intensity and blue, low intensity.
Both the co-occurrence matrix and the difference statistic angular second moment (a.s.m.) mapping routines were found to not be of much use. The angular second moment acts to enhance the brain boundary variations to such an extent that the internal variations are not visible, the intensity across the interface seemed to vary according to the neighbourhood size and not according to the image differences. A similar result for the co-occurrence correlation map images, both with and without cut-off, was found.
The element difference moment was examined with different order moments and at different angles and neighbourhoods (Fig 4.7) on C2. From these images it was decided to only examine the first moment, since the effect of the higher order moments was to multiply the difference between the co-occurrence pixels by the order, and therefore to simply shift the grey levels up by this amount. The differences between the different co-occurrence angles and directions was so difficult to identify that an angle of 45û and a distance of 1 were chosen, because of Haralick et al [14]. Neighbourhood values between 3 and 15 appeared to be of interest, since there appeared to be low resolution with larger neighbourhoods. The e.d.m. of the evaluation images was obtained for neighbourhoods of 3 (Fig 4.8) and 7. No deductions can be made from these images alone, however further texture mapping of these mapped images might reveal something of interest, but it is unlikely considering the nature of the e.d.m. which gives low values for neighbourhoods with pixels of a similar grey level.
Evaluation of the texture map for the inverse element difference moment, with different variables (Fig 4.9), revealed the same optimum value of the variables as for the e.d.m., with the exception of the neighbourhood, where a value of between 7 and 20 was perceived to be useful. The evaluation images were mapped using a neighbourhood size of 7 and 15 (Fig 4.10), and while further analysis and mapping seemed to be necessary on the size 7 images, the size 15 images seemed to be of more use. Rather than the simple variation in intensity across an interface being of interest, as was found for the histogram routines, the variation across the entire area is so. The temporal lobe and hippocampi regions of S1, S2 and M1 can be seen to be much reduced, when compared with C1 and C2, and the sulcal area of these three images seems to be more visible. The i.e.d.m. enhances the small intensity variations in the matter of the brain, and provides some differentiation between grey and white matter. The area variations do seem to be of interest, especially since they are extreme in the areas, the temporal lobes etc, where Alzheimer's disease has its initial effects.
C2 was texture mapped using the difference statistics contrast with different neighbourhood sizes and co-occurrence angles (Fig 4.11). This resulted in an angle of 45û and a neighbourhood of 15, being chosen, and the evaluation images being calculated for these values (Fig 4.12). While there are slight variations in the grey levels of the resulting images, these variations are two small to indicate much. The intensity does seem to be larger in the white matter regions around the temporal lobes for S1 and S2, but, although similar results were obtained for S3, it is not enough to allow any firm deductions.
The texture map of the difference statistics mean, with the same variables, proved to be of greater interest (Fig 4.13), probably because the lack of the square in the k term (Equation 2.13 as compared to 2.11) provides a desirable lower contrast. The connection between the high intensity around the temporal lobes and the low intensity at the upper interface between the brain matter and the CSF, for S1 and S2, compared to the lower lobe and higher upper interface intensities, for C1 and C2, provide a region of interest. High intensities correspond to a large number of very different grey levels in each neighbourhood. The cause of these variations is not known, but it can be theorised, by noting that the features of M1 and R1 are intermediate, that changes in the temporal lobe regions gives rise to a much higher water content around that area as compared to other areas in the brain.
The original moment routine was ignored since it varied only slightly from the central moment routine, and these variations were in favour of the central moment.
4.3.1 Central and Normalized Central Moments
The central moment of order (p+q) was mapped for C2 with a number of different neighbourhood sizes and various values of p and q (Fig 4.14). A value of p=q=1 was decided upon because it offered the clearest images and it also reduced time taken for the routine. The evaluation images were mapped for the optimum neighbourhood size of 11 (Fig 4.15). In the resulting images, pixels with a similar grey level had a similar neighbourhood in the original images. High grey level regions correspond to original neighbourhoods with large grey level values and large variations in grey level. The interconnectivity of the different areas of high intensity appear interesting. There are more medium intensity interconnections, for S1 and S2, with smaller high intensity spots, than for C1 and C2. This suggests that the difference in grey level and therefore water content between the grey and white matter is decreasing for Alzheimer's patients, which has be shown to be the case [15] for the complete brain histograms as well as demonstrated for the standard deviation and skew. Mapping C3, C4 and S3 confirmed these findings, however quantitative analysis as well as further texture mapping is required to substantiate these findings.
Images for the normalized central moment of order (1+1) were found to appear exactly the same as those for the central moment of order (1+1), with a neighbourhood of 11. However with a neighbourhood of 25, which decreases image resolution by a large amount, the two routines gave different results (Fig 4.16). The resulting images show bands of intensity, similar to those found for the histogram standard deviation, at the interface between the CSF and brain. This occurs because with a large neighbourhood the specific features are averaged out by the routine and it basically determines the grey level variation across the interface.
4.3.2 Invariant Moments
Each of the invariant moment mapping routines were examined using different neighbourhood sizes and on the different evaluation routines. Nothing could be determined by eye from any of the invariant moments except for moments 2 and 4.
The evaluation images were mapped using invariant moment 2, with an optimum neighbourhood size of 15, and the results were examined (Fig 4.17). In the resulting images regions of the same intensity corresponded to neighbourhoods which, when rotated or scaled, were the same. Once again more areas of high intensity can be seen, at the interface between CSF and brain, for the controls than for the Alzheimer's patients. However the intensity difference across the brain can be seen to be smaller in the Alzheimer's patients, especially S1, suggesting that there is a decrease in contrast between the grey and white matter. This is particularly obvious around the temporal lobes where there has been a high intensity region in the other mapped images.
The invariant moment 4, requiring a neighbourhood size of 15, was used to map the evaluation images (Fig 4.18). It can be seen that there are considerable differences between the controls and patients. The moment enhances the sulci and the hippocampi, by recognising the similar textures of their boundaries and, because of the large variation in intensity at the boundaries, gives these regions a high intensity. The intensity of the interfaces around S1 and S2 are much greater than those in the other images. The hippocampi in S1 and S2 also appear to be smaller and less visible than for the controls. Again, the quantitative analysis of many more images, is necessary to ensure these are in fact correct deductions.
The autocorrelation function of the each of the evaluation images was obtained (Fig 4.19). Readings of horizontal intensity showed that the autocorrelation of the two controls had a larger gradient than those of the Alzheimer's patients. If the images were single textures this would show that C1 and C2 are coarser textures than S1 and S2, a result that corresponds to the effects of Alzheimer's disease on grey and white matter variations. A fine texture would suggest a small variation between the water content of the predominant tissues in the brain, grey and white matter, whereas a coarse texture would suggest a more normal variation. While this is a very inexact hypothesis, it does tie in with the results so far.
THESE HAVE BEEN LEFT OUT SINCE NO MEANINGFUL RESULTS WERE OBTAINED... IF YOU WOULD LIKE TO SEE THEM EMAIL ME.
5.1 Further Possibilities of the Present Methods
The amount of possibilities available with the methods used in this project have made it impossible to cover even a fraction of one percent of them. To improve on the results obtained, the different routines should be used to evaluate as many different images from as many different patients, with as many different values of the user-definable variables, as possible. Some of the texture mapping methods showed remarkably good edge detection, these could be used for unsupervised minimisation of the head MR images to the brain, which would allow overnight multiple image minimisation.
Edge detection would also be useful in confining the neighbourhood within the different texture edges, thus providing an irregular neighbourhood, and also removing the texture overlap problem with the NTMM. The edge detection could be done by intensity or by texture, but it has been shown [10] that intensity edge detection methods work much better than their textural counterparts.
This project only looked at one stage texture mapping. Many of the texture maps could provide interesting results if they were mapped again using other methods, or if texture analysis was carried out on them. This could be expanded to provided a multistage mapping routine that would be combined with the original images to provide regions of interest and probable areas where the disease is starting, thus allowing high resolution MRI of these areas.
Alzheimer's disease is known to have an early effect on the temporal lobes. High resolution images of the temporal lobes could provide textural information that was not visible before, due to the lack of resolution. Resolution is the big problem with texture mapping these MRI images, the texture may be so small, cell size for example, that at the present resolution, no texture features are being perceived, just the orderly construction of the brain.
5.2 Other Methods and Variations
Using the two ideas of Julesz [16], that people cannot differentiate between texture with differences only in the third order statistic, and that while the majority of textural information is in the first two orders, higher order statistics do contain some information, it would seem wise to investigate higher order statistics, such as the run length matrix. The reasoning behind this is that it is very difficult to visually detect any texture variations in the MR images, and as this suggests that there is not much first or second order texture, the assumption could be that the majority of the textural information for these MR images is contained in the high order statistics. While these statistics are computationally intensive, they could provided interesting texture mapped images where a lot of time can be allowed for image processing.
The Fourier transform of an MR image contains quite a lot of high frequency components which almost certainly correspond to noise in the image. If a ring cut-off filter is used to remove the high frequency components and the reverse FFT is used to transform back from the frequency domain, while it is likely that this will remove much of the noise in the image, there is a possibility that this will degrade the image thus removing some of the resolution and texture.
All the methods in this project deal with a 2D image. If texture mapping and analysis are carried out on the 3D volume, the results would be much improved [21], [22]. A single image contains only 2D of the 3D the actual texture has and if this image is taken from the wrong place, the results could be misleading. While texture methods in this project can be extended to three dimensions with some difficulty, usually texture methods involving 3D masks work more successfully and involve simpler and faster computations.
5.3 Presentation of Images to the Clinician
The most important aspect of any imaging method is how the information can be presented to the clinician, from which a diagnosis or the treatment decision can be made. While clinicians in a technical field are usually good at dealing with technology and its results as whole, those in a non-technical environment are sometimes technophobic and adverse to making decisions based on the results of a machine. For this reason, any material from MRI, which has been manipulated into a new format by image manipulation, including texture mapping, needs to be given to the clinician in a consistent format. A result consisting of a pseudo-polydimensional enhanced image which can be compared with or incorporates the original, highlighting regions of interest, and generally providing qualitative information.
Until analysis methods are proven to be 100% reliable, clinicians will have to evaluate an enhancement by eye, because of the human element of confirmation that is necessary.
The different texture methods have been evaluated and the results examined. These results are not definite, and must therefore be viewed with some scepticism, due to the number of variables involved. But they do provide an indication of any interesting directions a future research project might want to go. Each of the five methods, histograms, co-occurrence matrix, difference matrix and autocorrelation function have some interesting features, however none of these seem to suggest what the root cause of Alzheimer's disease might be. They do, however, all seem to suggest variations in the intensity difference between grey and white matter. Suggesting that for Alzheimer's patients a relative increase in the grey matter water content, or decrease in the white matter water content, occurs.
The main limitation of this project has been time. This acted to limit the amount of image processing that could be undertaken, and more importantly has limited the total number of images samples that could be texture mapped, resulting in inconclusive results.
This project has covered a tiny area of a very large field of research, and although the results are inconclusive, they do provide a guide to the specific regions of interest. This provides any future researcher with ideas and results that can be used as a basis for new methods.
BIBLIOGRAPHY
[1] M. Amadasun and R. King. Textural features corresponding to textural properties. IEEE Transactions on Systems, Man and Cybernetics, 19:1264-74, September 1989.
[2] R. Antonczyk. MRI Texture Analysis. Imperial College MSc in Engineering and Physical Science in Medicine, 1993.
[3] F.Argenti et al. Fast calculation of co-occurrence matrix parameters for image segmentation. Electronics Letters, 26:23-4, 4 January 1990.
[4] F.Argenti et al. Fast algorithms for texture analysis using co-occurrence matrices. IEE Proceedings for Radar and Signal Processing, 137:443-8, December 1990.
[5] P.Brodatz.Texture: A Photographic Album for Artists and Designers. Dover, 1956.
[6] F.A. DeCosta and M.F. Chouika. Neural network recognition of texture images using third order cumulants as functional links. IEEE International Conference on Acoustics, Speech and Signal Processing, 3:77-80, March 1992.
[7] R. Finn. Quantitative Assessment of Neurological Image Characteristics in Alzheimer's Patients. MSc Thesis, Imperial College Centre for Biological and Medical Systems, 1993.
[8] R. Finn. Results of present work on the grey and white matter variations in the brain using histograms. Biomedical Systems Group, Imperial Colleg Electrical and Electronic Engineering Department.
[9] M.A. Foster and J.M.S.Hutchinson. Practical NMR Imaging. IRL Press, 1987.
[10] P.A. Freeborough. The Development of a Software Kernel for Medical Image Analysis and an Investigation into the Application of Texture Analysis to Clinical MRI. MPhil Transfer Report, Imperial College, May 1994.
[11] P.A. Freeborough. MIDAS 1.0 Application Programmer's Manual, January 1993.
[12] R.C. Gonzalez and R.E. Woods. Digital Image Processing: Addison-Wesley, 1992.
[13] R.M. Haralick. Statistical and structural approaches to texture. IEEE Proceedings, 67:768-804, 1979.
[14] R.M. Haralick et al. Textural features for image classification. IEEE Transactions on Systems, Man and Cybernetics, 3:610-621, November 1973.
[15] C.R. Jack et al. MR-based hippocampal volumetry in the diagnosis of Alzheimer's Disease. Neurology, 42:183-188, 1992.
[16] B. Julesz. Experiments in the visual perception of texture. Scientific American, 232:2-11, 1975.
[17] B. Julesz. Visual pattern discrimination. IRE Transactions on Information Theory, 8:84-92, 1962.
[18] P.J. Keller. Basic Principles of Magnetic Resonance Imaging. General Electric Medical Systems, 1988.
[19] K. Straughan. His orginal project description which provided an introduction that could not be improved upon.
[20] T.A. Mitchell and M. Sarhadi. A machine vision system using fast texture analysis for automated visual inspection. IEEE International Conference on Systems Engineering, 399-402, September 1992.
[21] M.C. de Oliveira. Texture Analysis and Pseudo-Colour Enhancement of MR Images for 3D Visualization. MPhil Transfer Report, Imperial College Electrical and Electronic Engineering Department, June 1991.
[22] M.C. de Oliveria and R.I. Kitney. Texture analysis for descrimination of tissues in MRI data. Proceedings Computers in Cardiology, 481-4, September 1991.
[23] J.P. Seab et al. Quantitive NMR measurement of hippocampal atrophy in Alzheimer's Disease. Magnetic Resonance in Medicine, 8:200-208, 1988.
[24] T.N. Tan and A.G.Constantinides. Texture analysis based on a human visual model.IEEE International Conference on Acoustics, Speech and Signal Processing, 4:2137-40, April 1990.
[25] T.N. Tan. Texture feature extraction via visual cortical channel modelling. IAPR International Conference on Pattern Recognition, 3:607-610, March 1992.
[26] F. Tomita and S. Tsuji. Computer Analysis of Visual Textures: Kluwer Academic Publishers, 1990.