Patent Publication Number: US-2012027275-A1

Title: Disease determination

Description:
The present invention relates to methods and apparatus suitable for use in the determination of the presence or absence of disease. More particularly, but not exclusively, the invention relates to methods for analysing retinal images to determine an indication of likelihood of disease. 
     Screening of large populations for early detection of indications of disease is common. The retina of the eye can be used to determine indications of disease, in particular diabetic retinopathy and macular degeneration. Screening for diabetic retinopathy is recognised as a cost-effective means of reducing the incidence of blindness in people with diabetes, and screening for macular degeneration is recognised as an effective way of reducing the incidence of blindness in the population more generally. 
     Diabetic retinopathy occurs as a result of vascular changes in the retina which cause swellings of capillaries known as microaneurysms and leakages of blood into the retina known as blot haemorrhages. Microaneurysms may eventually become a source of leakage of plasma causing thickening of the retina, known as oedema. If such thickening occurs in the macular region, this can cause loss of high quality vision. Retinal thickening is not easily visible in fundus photographs. Fat deposits known as exudates are associated with retinal thickening, and the presence of exudates may therefore be taken to be an indication of retinal thickening. Exudates are reflective and are therefore visible in retinal photographs. 
     A currently recommended examination technique for diabetic retinal screening uses digital fundus photography of the eye. Fundus images are examined by trained specialists to detect indicators of disease such as exudates, blot haemorrhages and microaneurysms as described above. The trained specialists determine whether a patient should be referred to an ophthalmologist based upon detected indicators of disease. This is time consuming and expensive. 
     A known manual assessment scheme is the Scottish Diabetic Retinopathy Grading Scheme (SDRGS). The SDRGS has a number of categories into which a patient is placed based upon the detection of lesions on the retina of the eye. For example if four or more blot haemorrhages are detected in one hemi-field of the eye, a patient is categorised as having referable background retinopathy and is referred to an ophthalmologist. Other similar manual assessment schemes exist in other geographic regions. 
     Automated image analysis may be used to reduce manual workloads in determining properties of images. Image analysis is now used in a variety of different fields. In particular, a variety of image analysis techniques are used to process medical images so as to provide data indicating whether an image includes features indicative of disease. Image analysis techniques for the processing of medical images in this way must be reliable both from the point of view of reliably detecting all features which are indicative of disease and from the point of view of not incorrectly detecting features which are not indicative of disease. 
     An image of the retina of the eye has a large number of features including blood vessels, the fovea, and the optic disc. An automated system that is able to distinguish between indicators of disease and normal features of the eye needs to take into account characteristics of the retina so as to properly distinguish features of a healthy eye from features which are indicative of disease. 
     Known manual grading schemes for determining whether a patient should be referred to an ophthalmologist such as the SDRGS are designed to use data generated from manual inspection of images by trained specialists who are highly skilled at accurately recognising indicators of disease. Known automated systems are often partially successful in identifying features in retinal images which are indicative of disease. However these known systems often fail to detect all retinal features of interest. In particular, known automated systems often only provide an indication of a confidence that a detected feature is of a particular type. As such, an automatic grading system which is to rely on computer generated lesion data should take this uncertainty into account. More generally, an automatic grading system that is able to effectively use automatically detected lesion data is desirable. 
     It is an object of some embodiments of the present invention to obviate or mitigate at least some of the problems set out above. 
     According to a first aspect of the invention there is provided a method of generating output data providing an indication of the presence or absence of disease from at least one retinal image of a patient. The method comprises receiving first data associated with detection of a first lesion type in the at least one image, receiving second data associated with detection of a second lesion type in the at least one image and receiving third data associated with detection of a third lesion type in the at least one image. The first data, the second data and the third data are combined to generate the output data providing an indication of the presence or absence of disease. 
     Combining a plurality of data associated with a plurality of lesion types in this way provides an improved indication of the presence or absence of disease. The data associated with each lesion type may give an indication of the presence of that lesion type in the image. 
     Each of the first data, second data and third data may comprise a value selected from a range of values on a continuous scale. 
     At least one of the first, second and third data may be a quantitative indication of associated with detection of a respective lesion type in the image. For example, the data associated with each lesion type may comprise a number of occurrences of that lesion type in the image or data indicating a confidence that a predetermined number of occurrences of that lesion type are included in the image. The confidence may be a value on a continuous scale of values. The first, second and third data associated may each be a number or a confidence, or some of the first, second and third data may be a number whilst others of the first, second and third data may be a confidence. 
     The first, second and third data may be automatically generated and as such may contain a degree of error. The combining of the data may comprise arithmetically combining the first, second and third data and the uncertainty associated with each of the first, second and third data may be mitigated by the arithmetic combination of the data associated with each lesion type. The arithmetic combination may be, for example, addition. 
     Said first data, said second data and said third data may have associated weights and said first data, said second data and said third data may be arithmetically combined in accordance with the respective associated weights. For example, each of said first, second and third data may be multiplied by an associated weight and the results of the multiplications may be added together to generate the output data. 
     Associating a weight with each of the first, second and third data allows the contribution of the data associated with different lesion types to a disease or no disease determination to be controlled. 
     The output data may be a value on a continuous scale of values. The continuous scale may be a continuous scale of integers or real numbers, or any other suitable continuous scale. 
     The method may further comprise processing said output data with reference to a threshold to generate Boolean data indicating the presence or absence of disease. 
     The at least one retinal image may comprise a first image and a second image, and the first image may be a retinal image of the left eye of said patient and the second image may be a retinal image of the right eye of said patient. The first data may be generated by selecting one of data associated with said first lesion type in said first image and data associated with said first lesion type in said second image or by combining data associated with said first lesion type in said first image and data associated with said first lesion type in said second image. The second data may be generated by selecting one of data associated with said second lesion type in said first image and data associated with said second lesion type in said second image or by combining data associated with said second lesion type in said first image and data associated with said second lesion type in said second image. The third data may be generated by selecting one of data associated with said third lesion type in said first image and data associated with said third lesion type in said second image or by combining data associated with said third lesion type in said first image and data associated with said third lesion type in said second image. 
     The first and second images may each be sets of images and the first, second and/or third data may be generated by selecting data associated with the respective lesion type in one of the sets of images or by combining data associated with the respective lesion type in one or more of the images in the first set of images and one or more of the images in the second set of images. For example, a maximum value may be selected across both sets of images or a maximum value in said first set of images may be combined with a maximum value in said second set of images. 
     The first lesion type may be microaneurysm. The first data may indicate a number of microaneurysms detected in said at least one retinal image. 
     The second lesion type may be blot haemorrhage. The second data may be generated from only a part of said at least one retinal image. The part of said at least one retinal image may be selected based upon a location of the centre of the fovea in the at least one retinal image. The part of the at least one retinal image may be a connected region of the retinal image. The part of the at least one retinal image may represent a convex region of the retinal image. The part of said retinal image may have a size determined based upon a size of an optic disc. The part of said at least one retinal image may be a substantially circular portion having a radius substantially equal to the diameter of said optic disc. The part of said at least one retinal image may be centred on the location of the centre of the fovea in the at least one retinal image. 
     The second data may be a sum of a plurality of confidence values, each confidence value being associated with a respective area of said at least one retinal image determined to be a possible blot haemorrhage, and indicating a confidence that said respective area represents a blot haemorrhage. The plurality of confidence values may be a predetermined number of confidence values. The second data may be a sum of the three largest confidence values associated with respective areas of said at least one retinal image determined to be a possible blot haemorrhage. 
     The third lesion type may be exudate. The third data may be generated from only a part of said at least one retinal image. The part of said at least one retinal image may be selected based upon a position of the centre of the fovea in the at least one retinal image. The part of said at least one retinal image may have a size determined based upon a size of an optic disc. The part of said at least one retinal image may be a substantially circular portion having a radius substantially equal to the diameter of said optic disc. The part of said at least one retinal image may be centred on the position of the centre of the fovea in the at least one retinal image. 
     The third data may be a sum of a plurality of confidence values, each confidence value being associated with a respective area of said at least one retinal image determined to be a possible exudate, and indicating a confidence that said respective area represents an exudate. The third data may indicate a sum of the three largest confidence values associated with respective areas of said at least one retinal image determined to be a possible exudate. 
     The method may further comprise generating fourth data indicating the presence of said third lesion type wherein said third data and said fourth data are generated from different parts of said retinal image. The fourth data may be generated from a part of said retinal image larger than the part of said retinal image used to generate said third data. The larger part of said retinal image may wholly enclose said part of said retinal image used to generate said third data. The larger part of said retinal image may be a substantially circular portion having a radius substantially equal to twice the diameter of an optic disc. 
     The size of an optic disc may be a standardised disc diameter obtained by taking the mean of measurements of the diameter of the optic disc in a plurality of images, each image having been obtained from a respective one of a plurality of subjects. Alternatively, the size of an optic disc may be calculated based upon the retinal image being processed. 
     The weights may be generated by receiving a plurality of data items, each data item comprising a plurality of data values and each data item being based upon a respective subject. For each data item classification data indicating the presence or absence of disease in the respective subject may be received. The plurality of data items may be processed so as to generate a weight for each data value, the weights being such that when applied to said data values of said data items an output is generated for each data item indicating the presence or absence of disease and the weights being generated such that the correspondence of said outputs with said classification data is maximised. 
     The method may further comprise generating at least one of said first data associated with said first lesion type, said second data associated with said second lesion type and said third data associated with said third lesion type. 
     The at least one image may be a plurality of images. The plurality of images may comprise a plurality of images taken from the same eye or may comprise one or more images taken from the right eye and one or more images taken from the left eye. 
     According to a second aspect of the invention there is provided a method of generating output data providing an indication of the presence or absence of disease from at least one retinal image of a patient. The method comprises receiving first data associated with detection of a first lesion type in the at least one image receiving second data associated with detection of a second lesion type in the at least one image, and receiving third data associated with detection of the second lesion type in the at least one image, wherein the second data and the third data are generated from different parts of the retinal image; and arithmetically combining the first data, the second data and the third data to generate the output data providing an indication of the presence or absence of disease. 
     The method may further comprise receiving fourth data associated with detection of a third lesion type in the at least one image, wherein the arithmetic combining combines the first data, the second data, the third data and the fourth data. 
     The second lesion type may be exudate. The third data may be generated from a part of the retinal image larger than the part of the retinal image used to generate the second data. 
     The larger part of the retinal image may wholly enclose the part of the retinal image used to generate the second data. The larger part of the retinal image may be a substantially circular portion having a radius substantially equal to twice the diameter of an optic disc and/or the part of the retinal image used to generate the second data may be a substantially circular portion having a radius substantially equal to the diameter of an optic disc. The diameter of an optic disc may be a standardised disc diameter obtained by taking the mean of measurements of the diameter of the optic disc in a plurality of images, each image having been obtained from a respective one of a plurality of subjects. 
     A further aspect of the invention provides a method of generating output data providing an indication of the presence or absence of disease from first and second retinal images. The method comprises generating first data by selecting one of data associated with a first lesion type in said first image and data associated with said first lesion type in said second image, generating second data by selecting one of data associated with a second lesion type in said first image and data associated with said second lesion type in said second image and generating third data by selecting one of data associated with a third lesion type in said first image and data associated with said third lesion type in said second image. Said first data, said second data and said third data are processed to generate said output data providing an indication of the presence or absence of disease. 
     Generating an indication of the presence or absence of disease based upon a combination of data associated with lesions in first and second retinal images in this way allows indicators of disease in different images from a patient to be considered. This allows lesion detection processes to be carried out on different images from a patient and so for example if a lesion can only be detected in one of the first and second retinal images, the determination of the presence or absence of disease should still be achieved. 
     Each of the first, second and third data may be generated by selecting a maximum of data associated with said first image and data associated with said second image. The first image may be a retinal image of the right eye of a patient and the second image may be a retinal image of the left eye of said patient. This allows indicators of disease in both of a patient&#39;s eyes to be considered. 
     The first image may comprise a plurality of first images and the second image may comprise a plurality of second images. This provides a method of selecting an image from a number of images from a patient that provide the greatest indication of disease. 
     A still further aspect of the invention provides a method of generating data providing an indication of the presence or absence of disease. The method comprises generating data indicating the presence of blot haemorrhages in an eye from only a part of a retinal image, wherein said part of said retinal image is a connected region of said retinal image and is selected based upon the location of an anatomical feature, and processing said data indicating the presence of blot haemorrhages in an eye to generate data providing an indication of disease. 
     Only searching a restricted area of the image can decrease the processing time required to generate the indication of the presence or absence of disease. This allows more subjects to be screened with the same resources. 
     The part of the retinal image may be generated centred on the anatomical feature. The part of the retinal image may include the anatomical feature. The anatomical feature may be the fovea and the part of said retinal image may be selected based upon a position of the centre of the fovea in the retinal image. The part of said retinal image may have a size determined based upon a size of an optic disc. The part of said retinal image may be a substantially circular portion having a radius substantially equal to the diameter of said optic disc. The part of said retinal image may be centred on the position of the centre of the fovea in the retinal image. The part of said retinal image may be a square portion of said retinal image. 
     A further aspect of the invention provides a method of generating a set of weights for use in generating data providing an indication of the presence or absence of disease from a retinal image. The method comprises receiving a plurality of data items, each data item comprising a plurality of data values, wherein a first data value of said plurality of data values is associated with detection of a first lesion in an image, a second data value of said plurality of data values is associated with detection of a second lesion type and a third data value of said plurality of data values is associated with detection of a third lesion type in an image, wherein at least one of said first data value, second data value and third data value is a quantitative indication associated with detection of a respective lesion type in the image, and each data item being based upon a respective subject. For each data item classification data indicating the presence or absence of disease in the respective subject is received. Said plurality of data items is processed so as to generate a weight for each data value, the weights being such that when applied to said data values of said data items an output is generated for each data item indicating the presence or absence of disease and the weights being generated such that the correspondence of said outputs with said classification data is maximised. 
     Generating weights in this way allows data providing an indication of the presence or absence of disease to be more effectively generated. The set of weights is preferably generated based upon a training set derived using the same techniques as data which is to be processed. This in turn allows an indication of disease to be determined that takes into account factors such as the particular apparatus used, the nature of the computer processing and the photographic protocol and generates a more accurate indication of disease. 
     At least some of said data values may comprise data indicating the presence or absence of a respective lesion type. At least some of said data values may comprise data indicating a confidence of the presence of a respective lesion type. At least some of said data values may comprise data indicating a number of occurrences of a respective lesion type. 
     The or each lesion type may be selected from the group consisting of microaneurysm, exudate and blot haemorrhage. Each data item may indicate characteristics of a retinal image taken from said subject. Said classification data may comprise a Boolean value indicating the presence or absence of disease. Each output may comprise a value on a continuous scale. 
     Aspects of the invention provide methods for processing a retinal image to determine whether the retinal image includes indicators of disease. In particular, it is known that the occurrence of microaneurysms, blot haemorrhages and exudates can be indicative of various disease conditions, and as such methods are provided in which the identification of microaneurysms, exudates and blot haemorrhages is applied to generate data indicating whether a processed retinal image includes indicators of disease. The processing of retinal images in this way can determine whether the retinal image includes indicators of any relevant disease. In particular, the methods can be used to detect indicators of diabetic retinopathy, age-related macular degeneration, cardio-vascular disease, and neurological disorders (for example Alzheimer&#39;s disease and stroke) although those skilled in the art will realise that the methods described herein can be used to detect indicators of any disease which are present in retinal images. 
     Aspects of the invention can be implemented in any convenient form. For example computer programs may be provided to carry out the methods described herein. Such computer programs may be carried on appropriate computer readable media which term includes appropriate tangible storage devices (e.g. discs). Aspects of the invention can also be implemented by way of appropriately programmed computers. 
    
    
     
       Embodiments of the present invention will now be described, by way of example only, with reference to the accompanying drawings, in which: 
         FIG. 1  is a schematic illustration of a system for analysis of retinal images according to an embodiment of the present invention; 
         FIG. 1A  is a schematic illustration showing a computer of the system of  FIG. 1  in further detail; 
         FIG. 2  is an example of a retinal image suitable for processing using the system of  FIG. 1 ; 
         FIG. 3  is a further example of a retinal image, showing the location of important anatomical features; 
         FIG. 4  is a flowchart showing processing carried out to identify features of an eye; 
         FIG. 5  is a flowchart showing a process for vessel enhancement used in identification of temporal arcades in a retinal image; 
         FIG. 6  is a flowchart of processing carried out to fit semi-ellipses to the temporal arcades; 
         FIG. 7  is a schematic illustration of an eye showing areas to be searched to locate the optic disc; 
         FIG. 8  is a flowchart showing processing carried out to locate the optic disc in a retinal image; 
         FIG. 9  is a schematic illustration of an eye showing location of the fovea relative to the optic disc; 
         FIG. 10  is a flowchart showing processing carried out to locate the fovea in a retinal image; 
         FIG. 11  is a flowchart showing processing carried out to identify blot haemorrhages in a retinal image; 
         FIG. 12  is a flowchart showing normalisation processing carried out in the processing of  FIG. 11 ; 
         FIG. 13  is a flow chart showing part of the processing of  FIG. 11  intended to identify candidate blot haemorrhages in further detail; 
         FIG. 14  is a series of retinal images showing application of the processing of  FIG. 13 ; 
         FIG. 15  is a flowchart showing a region growing process carried out as part of the processing of  FIG. 11 ; 
         FIG. 16  is a flowchart showing a watershed region growing process carried out as part of the processing of  FIG. 11 ; 
         FIG. 17  is a flowchart showing a vessel detection process carried out as part of the processing of  FIG. 11 ; 
         FIGS. 18A and 18B  are each a series of images showing application of the processing of  FIG. 17 ; 
         FIG. 19  is a flowchart showing a process for classification of a candidate region; 
         FIG. 20  is a flowchart showing processing carried out to identify exudates in a processed image; 
         FIG. 21  is a flowchart showing part of the processing of  FIG. 20  intended to identify candidate exudates in further detail; 
         FIG. 22  is a flowchart showing processing carried out to classify regions as exudates; 
         FIG. 23  is a graph showing a plurality of Receiver Operator Characteristic (ROC) curves obtained from results of application of the method of  FIG. 11 ; 
         FIG. 24  is a graph showing a plurality of ROC curves obtained from results of application of the method of  FIG. 20 ; 
         FIG. 25  is a schematic illustration of an arrangement in which the methods described herein can be employed; 
         FIG. 26  is a schematic illustration of a further arrangement in which the methods described herein can be employed; 
         FIG. 27  is a flowchart showing processing carried out in the arrangement of  FIG. 26 ; 
         FIG. 28  is a flowchart showing part of the processing of  FIG. 27  in further detail; 
         FIG. 29A  is a graph showing variance of the area under a ROC curve in dependence upon a plurality of weights at a first scale; 
         FIG. 29B  is a graph showing variance of the area under a ROC curve in dependence upon a plurality of weights at a second scale; 
         FIG. 30  is a graph showing a plurality of ROC curves obtained from results of application of the method of  FIG. 27 ; and 
         FIG. 31  is a table showing results obtained from application of the method of  FIG. 27  to a set of images. 
     
    
    
     Referring now to  FIG. 1 , a camera  1  is arranged to capture a digital image  2  of an eye  3 . The digital image  2  is a retinal image showing features of the retina of the eye  3 . The image  2  is stored in a database  4  for processing by a computer  5 . Images such as the image  2  of  FIG. 1  may be collected from a population for screening for a disease such as, for example, diabetic retinopathy. The camera  1  may be a fundus camera such as a Canon CR5-45NM from Canon Inc. Medical Equipment Business Group, Kanagawa, Japan, or any camera suitable for capturing an image of an eye. 
       FIG. 1A  shows the computer  5  in further detail. It can be seen that the computer comprises a CPU  5   a  which is configured to read and execute instructions stored in a volatile memory  5   b  which takes the form of a random access memory. The volatile memory  5   b  stores instructions for execution by the CPU  5   a  and data used by those instructions. For example, in use, the image  2  may be stored in the volatile memory  5   b.    
     The Computer  5  further comprises non-volatile storage in the form of a hard disc drive  5   c . The image  2  may be stored on the hard disc rive  5   c . The computer  5  further comprises an I/O interface  5   d  to which are connected peripheral devices used in connection with the computer  5 . More particularly, a display  5   e  is configured so as to display output from the computer  5 . The display  5   e  may, for example, display a representation of the image  2 . Additionally, the display  5   e  may display images generated by processing of the image  2 . Input devices are also connected to the I/O interface  5   d . Such input devices include a keyboard  5   f  and a mouse  5   g  which allow user interaction with the computer  5 . A network interface  5   h  allows the computer  5  to be connected to an appropriate computer network so as to receive and transmit data from and to other computing devices. The CPU  5   a , volatile memory  5   b , hard disc drive  5   c , I/O interface  5   d , and network interface  5   h , are connected together by a bus  5   i.    
     Referring now to  FIG. 2 , a retinal image  6  suitable for processing by the computer  5  of  FIG. 1  is shown. The image  6  shows a retina  7  upon which can be seen an optic disc  8  and blood vessels  9 . Further areas  10  can be seen and these further areas can be classified by human inspection. Some of these further areas  10  are indicative of disease and detection and identification of such areas is therefore desirable. Each further area  10  may be, amongst other things, a lesion such as a microaneurysm, a blot haemorrhage, an exudate, drusen, or an anatomical feature such as the optic disc, the macula or the fovea. 
       FIG. 3  shows a further image of an eye.  FIG. 3  shows the green plane of a colour image, the green plane having been selected because it allows lesions and anatomical features of interest to be seen most clearly. The optic disc  8  can again be seen. The optic disc is the entry point into the eye of the optic nerve and of retinal blood vessels  7 . It can be seen that the appearance of the optic disc is quite different from the appearance of other parts of the retina. Retinal blood vessels  7  enter the eye through the optic disc  8  and begin branching. It can be seen that the major blood vessels form generally semi-elliptical paths within the retina, and these paths are known as temporal arcades denoted  11 . The fovea  12  is enclosed by the temporal arcades, and is the region of the retina providing highest visual acuity due to the absence of blood vessels and the high density of cone photoreceptors. The fovea appears as a dark region on the surface of the retina, although its location can be masked by the presence of inter-retinal deposits known as drusen, as well as by exudates or cataract. The region surrounding the fovea  12  indicated  13  in  FIG. 3  is known as the macula. 
     The methods described below benefit from accurate location of the optic disc  8  and the fovea  12 . This is because areas of an image representing the optic disc  8 , the fovea  12  and the macula  13  need to be processed in particular ways. More specifically, artefacts which would normally be considered as indicators of disease are not so considered when they form part of the optic disc. It is therefore important to identify part of a processed image representing the optic disc so as to allow appropriate processing to be carried out. Additionally, it is known that the presence of lesions within the macula  13  has a particular prognostic significance. Furthermore the fovea could be falsely detected as a lesion if it is not identified separately. It is therefore also important to identify part of a processed image representing the fovea  12  and the surrounding macula  13 . 
     Methods for locating the optic disc  8  and fovea  12  in an input image are now described.  FIG. 4  shows the processing at a high level. First, at step S 1  an input image is processed to enhance the detectability of blood vessels. Then, at step S 2 , semi-ellipses are fitted to the blood vessels so as to locate the temporal arcades within the image. At step S 3  the image is processed to locate the optic disc  8 , the processing being limited to an area defined with reference to the temporal arcades  11 . At step S 4  the image is processed to locate the fovea  12 , the processing being limited to an area defined with reference to the temporal arcades  11  and the location of the optic disc  8 . 
     As indicated above, at step S 1  an input image is processed so as to enhance the visibility of blood vessels. This aids the location of the temporal arcades at step S 2 . If the original image is a colour image then the processing to enhance the visibility of blood vessels is carried out using the green colour plane. The process of vessel enhancement is described with reference to a flowchart shown in  FIG. 5 . 
     The processing of  FIG. 5 , as is described in further detail below, is arranged to enhance vessels on the basis of their linear structure. Vessels are detected at a plurality of different angles which are selected such that substantially all vessels can be properly enhanced. Vessels will generally satisfy the following criteria, which are used in the processing of  FIG. 5  as is described below:
         (i) an intensity gradient will exist at all pixels along each vessel wall;   (ii) intensity gradients across opposite vessel walls will be in approximately opposite directions; and   (iii) vessels are expected to have a range of widths, for example from 5 to 15 pixels depending on the scale of the image.       

     For improved efficiency, the optic disc and fovea can be detected in images which have been sub-sampled. For example, vessel enhancement does not require an image greater than about 500 pixels per dimension for a 45° retinal image. Different parts of the analysis can be carried out on images which have been subjected to sub-sampling. For this reason, in the following description, dimensions are expressed in terms of the expected optic disc diameter (DD) whose value should be taken to be relevant to the current possibly sub-sampled image. The value 1 DD is a standardised disc diameter obtained by taking the mean of, possibly manual, measurements of the diameter of the optic disc in several images. 
     Referring to  FIG. 5 , at step S 5  the input image is appropriately sub-sampled. An appropriate ratio for sub-sampling may be determined based upon the size of the input image. A counter n is initialised to a value of 0 at step S 6 . A variable θ is set according to equation (1) at step S 7 : 
     
       
         
           
             
               
                 
                   θ 
                   = 
                   
                     
                       n 
                        
                       
                           
                       
                        
                       π 
                     
                     18 
                   
                 
               
               
                 
                   ( 
                   1 
                   ) 
                 
               
             
           
         
       
     
     Subsequent processing is arranged to enhance vessels extending at the angle θ. θ′ is an angle perpendicular to the angle θ. That is: 
     
       
         
           
             
               
                 
                   
                     θ 
                     ′ 
                   
                   = 
                   
                     θ 
                     - 
                     
                       π 
                       2 
                     
                   
                 
               
               
                 
                   ( 
                   2 
                   ) 
                 
               
             
           
         
       
     
     A filter kernel L(θ′) is defined by a pixel approximation to a line such that the gradient in direction θ′ can be evaluated, using convolution of the image with this kernel. An example of L(θ′) is: 
         L (θ′)=[−3,−2,−1,0,1,2,3]  (3)
 
     The appropriately sub-sampled green plane of the input image is convolved with the linear kernel L(θ 1 ) at step S 8 , as indicated by equation (4): 
         e   θ ( x,y )= I ( x,y )* L (θ′)  (4)
 
     where * denotes convolution. 
     Given that the linear kernel L(θ′) is arranged to detect edges in a direction θ′, the image e θ  indicates the location of edges in the direction θ′ and consequently likely positions of vessel walls extending in the direction θ. As explained above, opposite walls will be indicated by gradients of opposite sign. That is, one wall will appear as a ridge of positive values while the other wall will appear as a ridge of negative values in the image output from equation (4). This is indicated by criterion (ii) above. 
     An image having pixel values greater than 0 at all pixels which are located centrally between two vessel walls satisfying criterion (ii) is generated at step S 9  according to equation (5): 
         g   θ,w ( x,y )=min( e   θ ( x+u   θ,w   ,y+v   θ,w ),− e   θ ( x−u   θ,w   ,y−v   θ,w ))  (5)
 
     The vector (u θ,w ,v θ,w ) is of length w/2 and extends in a direction perpendicular to the angle θ. w is selected, as discussed further below to indicate expected vessel width. 
     It can be seen that a value for a particular pixel (x,y) in the output image is determined by taking the minimum of two values of pixels in the image e θ . A first pixel in the image e θ  is selected to be positioned relative to the pixel (x,y) by the vector (u θ,w ,v θ,w ) while a second pixel in the image (the value of which is inverted) is positioned relative to the pixel (x,y) by the vector −(u θ,w ,v θ,w ). Equation (5) therefore means that a pixel (x,y) in the output image g has a positive value only if the pixel at (x+u θ,w ,y+v θ,w ) has a positive value and the pixel at (x−u θ,w ,y−v θ,w ) has a negative value. Thus, equation (5) generates a positive value for pixels which are located between two edges, one indicated by positive values and one indicated by negative values, the edges being separated by the value w. 
     It can be appreciated that the value of w should be selected to be properly indicative of vessel width. No single value of w was found to enhance all vessels of interest. Therefore, applying processing with value of w of 9 and 13 has been found to provide acceptable results. 
     The preceding processing is generally arranged to identify vessels. However both noise and vessel segments extending at an angle θ will produce positive values in the output image g θ . Noise removal is performed by applying morphological erosion with a linear structuring element s(θ,λ), approximating a straight line of length λ extending at an angle θ, to the output image g θ . After morphological erosion a pixel retains its positive value only if all pixels in a line of length λ extending at the angle θ centered on that pixel also have positive values. 
     A greater value of λ increases noise removal but reduces the proportion of vessels that are properly enhanced. A value of λ=21 for a 45° image having dimensions of about 500 pixels (or 0.18 DD more generally) has been found to give good results in experiments. 
     Referring again to  FIG. 5  it will be recalled that at step S 9  an output image g θ,w  was formed. At step S 10 , an output image V θ  is created in which each pixel has a value given by the maximum of the corresponding pixels in two images created with different values of w (9 and 13) when eroded with the described structuring element s(θ,λ). This is expressed by equation (6): 
         V   θ =max[ε s(θ,21 ) g   θ,9 ( x,y ),ε s(θ,21)   g   θ,13 ( x,y )]  (6)
 
     At step S 11  a check is carried out to determine whether the value of n is equal to 17, if this is not the case, processing passes to step S 12  where the value of n is incremented before processing returns to step S 7  and is repeated in the manner described above. In this way, it can be seen that eighteen images V θ  are created for different values of θ. 
     When it is determined at step S 11  that processing has been carried out for all values of n which are of interest, processing continues at step S 13  where the maximum value of each pixel in all eighteen images V θ , is found so as to provide a value for that pixel in an output image V. At step S 14  the angle producing the maximum value at each pixel is determined to produce an output image Φ. That is, the output image Φ indicates the angle θ which resulted in each pixel of the image V having its value. 
     The processing described with reference to  FIG. 5  is arranged to produce an image in which vessels are enhanced. It will be recalled that it is desired to locate the semi-elliptical temporal arcades, as indicated by step S 2  of  FIG. 4 . This is achieved by applying a generalized Hough transform (GHT) to the images V and Φ. Use of the generalized Hough transform is explained in Ballard, D. H.: “Generalizing the Hough transform to detect arbitrary shapes”, Pattern Recognition, 13, 111-122, the contents of which are incorporated herein by reference. 
     The application of the GHT is shown, at a high level, in  FIG. 6 . 
     At step S 15  an image V +  is formed from the image V according to equation (7): 
     
       
         
           
             
               
                 
                   
                     
                       V 
                       + 
                     
                      
                     
                       ( 
                       
                         x 
                         , 
                         y 
                       
                       ) 
                     
                   
                   = 
                   
                     { 
                     
                       
                         
                           
                             
                               0 
                                
                               
                                   
                               
                                
                               if 
                                
                               
                                   
                               
                                
                               
                                 V 
                                  
                                 
                                   ( 
                                   
                                     x 
                                     , 
                                     y 
                                   
                                   ) 
                                 
                               
                             
                             ≤ 
                             0 
                           
                         
                       
                       
                         
                           
                             
                               1 
                                
                               
                                   
                               
                                
                               if 
                                
                               
                                   
                               
                                
                               
                                 V 
                                  
                                 
                                   ( 
                                   
                                     x 
                                     , 
                                     y 
                                   
                                   ) 
                                 
                               
                             
                             &gt; 
                             0 
                           
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   7 
                   ) 
                 
               
             
           
         
       
     
     The image V +  is then skeletonised at step S 16  to form an image U. That is: 
         U=SKEL ( V   + )  (8)
 
     To achieve acceptable execution times of the GHT, images V and Φ may need to be greatly sub-sampled. Tests have shown that the GHT performs satisfactorily after U and Φ have been sub-sampled to have each dimension being approximately 50 pixels. At step S 17  the image U is Gaussian filtered and at steps S 18  and S 19  the images U and Φ are appropriately sub-sampled. 
     At step S 20  the GHT is applied to the images U and Φ to locate vessels following semi-elliptical paths. 
     To enable acceptable execution time and memory usage Hough space is discretized, for example as five dimensions, as follows:
         p takes an integer value between 1 and 45 and is an index indicating a combination of ellipse aspect ratio and inclination;   q takes an integer value between 1 and 7 and is an index for a set of horizontal axis lengths linearly spaced from 23.5 to 55 sub-sampled pixels, at the sub-sampled resolution of U′;   h takes an integer value of 1 or 2 and indicates whether the semi-ellipse is the left or right hand part of a full ellipse; and   (a,b) is the location within the image of the centre of the ellipse.       

     Only some combinations of p and q are useful, given known features of retinal anatomy. For example, combinations of p and q giving rise to an ellipse whose nearest to vertical axis is longer than the anatomical reality of the temporal arcades are discarded. 
     The use of the GHT to locate the temporal arcades as described above can be made more efficient by the use of templates, as is described in Fleming, A, D,: “Automatic detection of retinal anatomy to assist diabetic retinopathy screening”, Physics in Medicine and Biology, 52 (2007), which is herein incorporated by reference in its entirety. Indeed, others of the techniques described herein for locating anatomical features of interest are also described in this aforementioned publication. 
       FIG. 7  is a schematic illustration of an eye, showing blood vessels  7   a  making up the temporal arcades. Two of the semi-ellipses  14  fitted using the processing described above are also shown. The semi-ellipses are used to restrict the search carried out at step S 3  of  FIG. 4  to locate the optic disc. 
     Experiments have shown that the optic disc is likely to lie near the right or left most point of the semi-ellipses. Experiments using training images also found that at least one point of vertical tangent of the three semi-ellipses defined in Hough space by (p n , q n , h n , a n , b n ) where n=1, 2, 3 was close to the position optic disc. The centre of the optic disc usually lies within an ellipse having a vertical height of 2.4 DD and a horizontal width of 2.0 DD centred on one of these points. Therefore, the union of the ellipses centred the point of vertical tangent of the three ellipses indicated above was used as a search region. 
     Referring again to  FIG. 7 , it can be seen that a point  15   a  on a semi-ellipse  14   a  has a vertical tangent, as does a point  15   b  on a semi-ellipse  14   b . An ellipse  16   a  having the experimentally determined dimensions centred on the point  15   a  is also shown, as is an ellipse  16   b  centred on the point  15   b . The union of the two ellipses (together with a third ellipse not shown in  FIG. 7  for the sake of simplicity) defines the area which is to be searched for the location of the optic disc. 
     A weighting function, W OD  is defined to appropriately limit the search area, such that all pixels outside the region of interest defined with reference to the union of ellipses described above have a zero weighting. 
     Within the search area, the optic disc is located using a circular form of the Hough transform, as is now described with reference to  FIG. 8 . Processing efficiency can be improved by sub-sampling the image. First, at step S 25  an anti-aliasing filter is applied to the input image. The optic disc is usually most easily detected in the green colour plane of a colour image. However in some cases, detection is easier in the red colour plane, and as such, at step S 26  both the green and red colour planes are sub-sampled to give image dimensions of about 250 pixels for a 45° fundus image, so as to improve processing efficiency. Gradient images are then formed by applying a Sobel convolution operator to each of the sub-sampled red and green planes at step S 27 . In order to remove the influence of vessels in the gradient images, a morphological closing is applied with a circular structuring element having a diameter larger than the width of the largest blood vessels but much smaller than the expected optic disc size at step S 28 . This morphological closing removes vessels but has little effect on the optic disc because it is usually an isolated bright object. Each gradient image after morphological closing is convolved with a Gaussian low pass filter with σ=1. 
     At step S 30 , the filtered gradient images produced at step S 29  from each of the red and green colour planes are combined, such that the value of any pixel in the combined image is the maximum value of that pixel in either the two filtered gradient images generated by processing the red and green image planes. 
     At step S 31  a threshold is applied to the image created at step S 30  so as to select the upper quintile (20%) of pixels with the greatest gradient magnitude. This threshold removes noise while maintaining pixels at the edge of the optic disc. 
     A circular Hough transform is applied to the image generated at step S 31  so as to locate the optic disc. The variety of radii for the optic disc observed in training images mean that the Hough transform is applied for a variety of radii. More specifically, nine radii arranged in a linear sequence between 0.7 DD and 1.25 DD were used. Experiments have shown that such radii represent 99% of actual disc diameters experienced. Using local gradient x and y components, the position of the optic disc centre can be estimated for each supposed pixel on the boundary of the optic disc and for each radius value. This means that, for each pixel, only a single Hough space accumulator need be incremented per radius value. Uncertainty in the location and inclination of the optic disc boundary is handled by applying a point spread function to the Hough space, which can be achieved by convolution with a disc of about ⅓ DD in diameter. 
     The optic disc location is generated at step S 33  as the maximum in Hough space from the preceding processing, bearing in mind the limitation of the search area as described above. 
     Referring back to  FIG. 4 , it was explained that at step S 4  the image is processed so as to locate the fovea. This is now described in further detail. The process involves locating a point in an input image which is most likely to represent the location of the centre of the fovea based upon a model of expected fovea appearance. The search is limited to a particular part of an input image, as is now described with reference to  FIG. 9 . 
       FIG. 9  is a schematic illustration of an eye showing a semi-ellipse  15  fitted using the GHT as described above. The optic disc  8  is also shown, together with its centre (x O , y O ) as located using processing described with reference to  FIG. 8 . The centre of a region to be used as a basis for location of the fovea is indicated (x F     —     EST , Y F     —     EST ). This point is positioned on a line extending from the centre of the optic disc (x O , y O ) to the centre of the semi-ellipse having centre (a 1 , b 1 ) as identified using the GHT. The centre of the region to be used as a basis for search is located 2.4 DD from the optic disc centre. A circular region  16  expected to contain the fovea has a diameter of 1.6 DD. The size of the region expected to contain the fovea, and its location relative to the optic disc were determined empirically using training images. 
     Processing carried out to locate the fovea is now described with reference to  FIG. 10 . This processing uses the green plane of an image of interest, and the green plane is sub-sampled at an appropriate ratio (down to a dimension of about 250 pixels) at step S 35  to produce an image I so as to improve processing efficiency. The sub-sampled image is then bandpass filtered at step S 36 . The attenuation of low frequencies improves detection by reducing intensity variations caused by uneven illumination and pigmentation. The removal of high frequencies removes small scale intensity variations and noise, which can be detrimental to fovea detection. The filtering is as set out in equation (9): 
         I   bpf   =I   lpf   −I   lpf   *gauss (0.7 DD )  (9)
 
     where: 
     I bpf  is the output bandpass filtered image; 
     gauss(σ) is a two-dimensional Gaussian function with variance σ 2 ; 
     I lpf =I*gauss(0.15 DD); and 
     I is the sub-sampled green plane of the input image. 
     At step S 37  all local minima in the bandpass filtered image are identified, and intensity based region growing is applied to each minima at step S 38 . The region generated by the region growing process is the largest possible connected region such that it includes the minimum of interest, and such that all pixels contained in it have an intensity which is less than or equal to a certain threshold. This threshold can be determined for example by taking the mean intensity in a circular region with a diameter of about 0.6 DD surrounding the minimum of interest. 
     Regions having an area of more than about 2.3 times the area of a standard optic disc are discarded from further processing on the basis that such areas are too large to be the fovea. Regions which include further identified minima are also discarded. 
     At step S 39  regions which do not intersect the circular region  16  expected to contain the fovea (as described above with reference to page 9) are discarded from further processing. At step S 40  a check is carried out to determine whether there are any regions remaining after the discarding of step S 39 . If this is not the case, the approximated position of the expected position of the fovea relative to the optic disc (x F     —     EST , y F     —     EST ) is used as the approximate location of the fovea at step S 41 . Otherwise, processing passes from step S 40  to step S 42  where regions intersecting the area in which the fovea is expected are compared with a predetermined model of the fovea which approximates the intensity profile of the fovea in good quality training images. The model has a radius R of 0.6 DD and is defined as: 
         M ( x,y )= B ( A −√{square root over (x 2   +y   2 )})  (10)
 
     where: 
     (x, y)εdisc(R); 
     disc(R) is the set of pixels within a circle of radius R centered on the origin; and 
     A and B are chosen so that the mean and standard deviation of M over disc(R) are 0 and 1 respectively. 
     The comparison of step S 42  is based upon a correlation represented by equation (11): 
     
       
         
           
             
               
                 
                   
                     C 
                      
                     
                       ( 
                       
                         x 
                         , 
                         y 
                       
                       ) 
                     
                   
                   = 
                   
                     
                       
                         ∑ 
                         
                           
                             ( 
                             
                               a 
                               , 
                               b 
                             
                             ) 
                           
                           ∈ 
                           
                             disc 
                              
                             
                               ( 
                               R 
                               ) 
                             
                           
                         
                       
                        
                       
                         
                           
                             I 
                             bpf 
                           
                            
                           
                             ( 
                             
                               
                                 x 
                                 + 
                                 a 
                               
                               , 
                               
                                 y 
                                 + 
                                 b 
                               
                             
                             ) 
                           
                         
                          
                         
                           M 
                            
                           
                             ( 
                             
                               a 
                               , 
                               b 
                             
                             ) 
                           
                         
                       
                     
                     
                       
                         
                           
                             
                               
                                 
                                   ∑ 
                                   
                                     
                                       ( 
                                       
                                         a 
                                         , 
                                         b 
                                       
                                       ) 
                                     
                                     ∈ 
                                     
                                       disc 
                                        
                                       
                                         ( 
                                         R 
                                         ) 
                                       
                                     
                                   
                                 
                                  
                                 
                                   
                                     
                                       I 
                                       bpf 
                                     
                                      
                                     
                                       ( 
                                       
                                         
                                           x 
                                           + 
                                           a 
                                         
                                         , 
                                         
                                           y 
                                           + 
                                           b 
                                         
                                       
                                       ) 
                                     
                                   
                                   2 
                                 
                               
                               - 
                             
                           
                         
                         
                           
                             
                               
                                 
                                   1 
                                   N 
                                 
                                  
                                 
                                   [ 
                                   
                                     
                                       ∑ 
                                       
                                         
                                           ( 
                                           
                                             a 
                                             , 
                                             b 
                                           
                                           ) 
                                         
                                         ∈ 
                                         
                                           disk 
                                            
                                           
                                             ( 
                                             R 
                                             ) 
                                           
                                         
                                       
                                     
                                      
                                     
                                       
                                         I 
                                         bpf 
                                       
                                        
                                       
                                         ( 
                                         
                                           
                                             x 
                                             + 
                                             a 
                                           
                                           , 
                                           
                                             y 
                                             + 
                                             b 
                                           
                                         
                                         ) 
                                       
                                     
                                   
                                   ] 
                                 
                               
                               2 
                             
                           
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   11 
                   ) 
                 
               
             
           
         
       
     
     Where N is the number of pixels in disc(R) and the mean of C is calculated for all pixels in a particular region. 
     Having determined a value indicative of the correlation of each region with the model at step S 42 , processing passes to step S 43 , where the candidate having the largest calculated value is considered to be the region containing the fovea, and the centroid of that region is used as the centre of the fovea in future analysis. 
     The preceding description has been concerned with processing images to identify anatomical features. As described above, the identification of such anatomical features can be useful in the processing of images to identify lesions which are indicative of the presence of disease. One such lesion which can be usefully identified is a blot haemorrhage. 
     Referring now to  FIG. 11 , processing to identify blot haemorrhages in an image is shown. At step S 51  an image A corresponding to image  2  of  FIG. 1  is input to the computer  5  of  FIG. 1  for processing. At step S 52  the image A is normalised as described in further detail below with reference to  FIG. 12  and at step S 53  the normalised image is processed to detect points which are to be treated as candidate blot haemorrhages as described in further detail below with reference to  FIG. 13 . Candidate blot haemorrhages are returned as a single pixel location in the original image A. At step S 54  the candidate blot haemorrhage pixels identified at step S 53  are subjected to region growing to determine the region of the image A that is a possible blot haemorrhage region corresponding to the identified candidate pixel, as described in further detail below with reference to  FIG. 15 . 
     At step S 55  a region surrounding the region grown at step S 54  is grown (using a technique called “watershed retinal region growing”) such that it can be used in determining properties of the background of the area which is considered to be a candidate blot haemorrhage, as described in further detail below with reference to  FIG. 16 . 
     At step S 56  a region surrounding each identified candidate region is processed to locate structures which may be blood vessels as described in further detail below with reference to  FIG. 17 . Areas where vessels, such as blood vessels  9  of  FIG. 2 , cross can appear as dark regions similar to the dark regions associated with a blot haemorrhage. It is possible to identify areas where vessels cross and this information can be useful in differentiating candidate regions which are blot haemorrhages from other dark regions caused by vessel intersection. 
     At step S 57  each identified candidate blot haemorrhage is processed to generate a feature vector. Features that are evaluated to generate the feature vector include properties of the candidate region together with features determined from the vessel detection of step S 56  and the watershed region growing of step S 55 . 
     At step S 58  each candidate blot haemorrhage is processed with reference to the data of step S 57  to determine a likelihood that a candidate is a blot haemorrhage. The determination is based upon the feature vector determined at step S 57  together with additional information with regard to the location of the fovea which can be obtained using the processing described above. The processing of steps S 57  and S 58  are described in further detail below with reference to  FIG. 19 . 
     Either one or zero candidates within 100 pixels of the fovea is classified as the fovea and removed from the set of candidate blot haemorrhages. All other candidates are then classified according to a two-class classification that produces a likelihood that each candidate is a blot haemorrhage or background. The two-class classification uses a support vector machine (SVM) trained on a set of hand-classified images. 
     Referring now to  FIG. 12 , processing carried out to normalise an image A at step S 52  of  FIG. 11  is described. At step S 60  the original image A is scaled so that the vertical dimension of the visible fundus is approximately 1400 pixels for a 45 degree fundus image. At step S 61  the scaled image is filtered to remove noise. The filtering to remove noise comprises first applying a 3×3 median filter, which removes non-linear noise from the input image, and second convolving the median-filtered image with a Gaussian filter with a value of σ=2. An image I is output from the processing of step S 61 . 
     At step S 62  an image of the background intensity K is estimated by applying a 121*121 median filter to the image A. Applying a median filter of a large size in this way has the effect of smoothing the whole image to form an estimate of the background intensity. 
     At step S 63  a shade-corrected image is generated by pixel-wise dividing the pixels of the noise-reduced image generated at step S 61  by the image K generated at step S 62  and pixel-wise subtracting 1. That is: 
     
       
         
           
             
               
                 
                   
                     
                       J 
                       ′ 
                     
                      
                     
                       ( 
                       
                         x 
                         , 
                         y 
                       
                       ) 
                     
                   
                   = 
                   
                     
                       
                         I 
                          
                         
                           ( 
                           
                             x 
                             , 
                             y 
                           
                           ) 
                         
                       
                       
                         K 
                          
                         
                           ( 
                           
                             x 
                             , 
                             y 
                           
                           ) 
                         
                       
                     
                     - 
                     1 
                   
                 
               
               
                 
                   ( 
                   12 
                   ) 
                 
               
             
           
         
       
     
     Where I and K are as defined above, and J′ is the output shade-corrected image. Subtracting the value 1 makes the background intensity of the image equal to zero objects darker than the background have negative values and objects brighter than the background have positive values which provides an intuitive representation but is not necessary in terms of the image processing and can be omitted in some embodiments. 
     At step S 64  the resulting image is normalised for global image contrast by dividing the shade-corrected image pixel-wise by the standard deviation of the pixels in the image. That is: 
     
       
         
           
             
               
                 
                   
                     J 
                      
                     
                       ( 
                       
                         x 
                         , 
                         y 
                       
                       ) 
                     
                   
                   = 
                   
                     
                       
                         J 
                         ′ 
                       
                        
                       
                         ( 
                         
                           x 
                           , 
                           y 
                         
                         ) 
                       
                     
                     
                       sd 
                        
                       
                         ( 
                         
                           J 
                           ′ 
                         
                         ) 
                       
                     
                   
                 
               
               
                 
                   ( 
                   13 
                   ) 
                 
               
             
           
         
       
     
       FIG. 13  shows the detection of candidate blot haemorrhages at step S 53  of  FIG. 11 . At step S 65  the normalised image J output from the processing of  FIG. 12  is smoothed by applying an anti-aliasing filter and at step S 66  a counter value s is set to 0. The image J is processed at successively increasing scales meaning the size of objects detected at each iteration tends to increase. The scaling can be carried out by reducing the image size by a constant factor such as √{square root over (2)} at each iteration and then applying the same detection procedure at each iteration. The counter value s counts through the scales at which the image J is processed as described below. At each scale, candidate regions are identified which tend to represent larger lesions as the image size is reduced. 
     At step S 67  an image J 0  representing the un-scaled image is assigned to the input image J. At step S 68  a counter variable n is assigned to the value 0 and at step S 69  a linear structuring element L n  is determined according to equation (14) below: 
         L   n =Λ( p,nπ/ 8)  (14)
 
     where p is the number of pixels in the linear structuring element and Λ is a function that takes a number of pixels p and an angle and returns a linear structure comprising p pixels which extends at the specified angle. It has been found that a value of p=15 is effective in the processing described here. 
     At step S 70  an image M n  is determined where M n  is the morphological opening of the inverted image J s  with the structuring element L n . The morphological opening calculated at step S 70  is defined according to equation (15) below, 
         M   n   =−J   s   ∘L   n   (15)
 
     where −J s  is the inversion of the image at scale s, L n  is the linear structuring element defined in equation (14) and ∘ represents morphological opening. 
     In the image M n , areas that are possible candidate blot haemorrhages, at the current scale, are removed and areas that correspond to vessels and other linear structures extending approximately at an angle nπ/8 are retained because the morphological opening operator removes structures which are not wholly enclosed by the structuring element. Since a linear structuring element is used, this means structures in the image that are not linear are removed, thus resulting in the removal of areas which are dark in J excluding vessel structures approximately at angle nir/8 but including the removal of candidate blot haemorrhages. 
     At step S 71  it is determined if n is equal to 7. If n is not equal to 7 then at step S 72  n is incremented and processing continues at step S 69 . If it is determined at step S 71  that n is equal to 7 then processing continues at step S 73  as described below. 
     The processing of steps S 69  to S 72  creates eight structuring elements which are arranged at eight equally spaced orientations. Applying these eight structuring elements to the image −J s  creates eight morphologically opened images, M n , each image only including vessels extending at a particular orientation, the orientation being dependent upon the value of n. Therefore, the pixel-wise maximum of M n =0 . . . 7 includes vessels at all orientations. 
     At step S 73  an image D s  is generated by subtracting pixel-wise the maximum corresponding pixel across the set of images M n , for n in the range 0 to 7, from the inverted image −J s . Given that each of the images M n  contains only linear structures extending in a direction close to one of the eight orientations nπ/8, it can be seen that the subtraction results in the removal from the image of all linear structures extending close to one of the eight orientations which is generally equivalent to removing linear structures at any orientation. This means that the image D s  is an enhancement of dark dots, at the current scale s, present in the original image with vessels removed and candidate blot haemorrhages retained. 
     As indicated above, an input image is processed at a variety of different scales. Eight scales are used in the described embodiment. The counter s counts through the different scales. At step S 74  it is determined if s is equal to 7. If it is determined that s is not equal to 7 then there are further scales of the image to be processed and at step S 75  the counter s is incremented. 
     At step S 76  an image J S  is determined by morphologically closing the image J s-1  with a 3×3 structuring element B and resizing this image using a scaling factor, √{square root over (2)}. The structuring element may be a square or approximately circular element and applying the element in this way eliminates dark areas which have at least one dimension with small extent. In particular, closing by the structuring element B removes or reduces the contrast of vessels in the image whose width is narrow compared to the size of the structuring element. Reducing the contrast of vessels can reduce the number of erroneously detected candidate blot haemorrhages. Closing by structuring element B, at each iteration, is particularly important when the morphological processing, at step S 73 , which distinguishes blot like objects from linear vessels, is applied at multiple scales. This is because when processing is carried out to identify large lesions, and the image is much reduced in size, the linear structuring element no longer fits within the scaled vessels, and as such large lesions are more easily detected. 
     The processing of steps S 68  to S 76  is then repeated with the image as scaled at step S 78 . The scaling function therefore reduces the size of the image so that each time the image is scaled larger candidates are detected, the scaling being applied to the closure of the image processed at the previous iteration. 
     The scaling and morphological closure with the structuring element B of step S 76  can be defined mathematically by equation (16): 
         J   s ( x,y )=[ J   s-1   •B ]√{square root over (2)} x ,√{square root over (2)} y )  (16)
 
     where • is morphological closure. 
     If it is determined at step S 74  that s is equal to 7 then at step S 77  candidate blot haemorrhages are determined by taking the maximum pixel value of the images D s  for s in the range 0 to 7 for each pixel of the image and determining if the resulting maximum value for a particular pixel is above an empirically determined threshold T to determine whether that pixel is to be considered to be a blot haemorrhage. A suitable value for T is 2.75 times the 90 th  percentile of the maxima. 
     At step S 78  a candidate haemorrhage is determined for each connected region consisting entirely of pixels having pixel value greater than T. For each of these regions, the pixels contained within the region are searched for the pixel which is darkest in the shade-corrected image J. At step S 78  the darkest pixel in the region is added to a set of candidates C. A pixel taken to indicate candidate haemorrhage is thus selected for each of the regions. For each pixel that it is determined at step S 77  that the maximum pixel value of the images D s  has a value less than T, at step S 79  the pixel is determined to not be a candidate. 
     Some example images showing stages in the processing of  FIG. 13  will now be described. 
     Referring to  FIG. 14 , five stages of processing a particular image according to  FIG. 5  are shown. Image (i) shows the original image portion which contains vessels  20 ,  21  and large dark areas  22 ,  23 ,  24 . It is desired to identify the large regions as candidate blot haemorrhages and to not identify the vessels as candidate blot haemorrhages. 
     Image (ii) shows an image D 1  created using the processing described above. The image area shown in the Image (ii) is the same as that of the Image (i). D 1  is the image processed at the smallest scale and it can be seen that only small regions have been identified. 
     Image (iii) shows the image −J 8 , that is the image at the largest scale after scaling and morphological closing with the structuring element B, and after inversion (as can be seen by the dark areas appearing bright and the relatively bright background showing as dark). At this largest scale (s equal to 8) only the largest dark area of the original image appears bright. 
     Image (iv) shows the result of combining D s  for all values of s and is the image to which thresholding is applied at step S 79 . It can be seen in the image (iv) that three areas  25 ,  26 ,  27  appear bright that correspond to the dark areas  22 ,  23 ,  24 . 
     The darkest pixels in the areas of the original image corresponding to bright areas such as areas  25 ,  26 ,  27  are added to a candidate set C. Region growing to identify the region of the original image is then performed from these candidates, and will now be described with reference to  FIG. 15 . 
     Referring to  FIG. 15 , at step S 85  a candidate c is selected from the candidate set C that has not previously been selected. At step S 86  a threshold t is set to a value of 0.1. At step S 87  a region C t  of the original image is determined such that C t  is the largest connected region (defined in a particular embodiment using orthogonal and diagonal adjacency) containing c and satisfying equation (18) shown below: 
         J ( p )≦ J ( c )+ t, ∀pεC   t   (18)
 
     where J is the normalised original image determined at step S 52  of  FIG. 11  and J(p) is pixel p of image J. 
     The area C t  determined according to the inequality of equation (18) is a collection of connected pixels of the original image in which each pixel in the area is less dark than the darkest pixel by no more than the value t. 
     At step S 88  it is determined if the number of pixels in the area C t  is less than 3000 pixels. If it is determined that the number of pixels is less than 3000 in the area C t  then at step S 89  area C t  is added to a set S and at step S 90  the threshold t is increased by a value of 0.1. Processing then continues at step S 86  as described above. 
     The loop of steps S 86  to S 90  identifies a plurality of increasingly large regions of pixels that are relatively dark when compared to the pixels that lie on the outside of the selected region. Each time the threshold t is increased, pixels that are connected to the region containing the seed pixel c that are less dark than allowed into the region by the previous value of t are included in the area C t . If it is determined at step S 88  that the number of pixels that are in the region determined by the threshold t is greater than 3000 then it is determined that the area allowed by the threshold t is too large and processing continues at step S 91 . 
     At step S 91  an energy function is used to determine an energy associated with a particular threshold t: 
         E ( t )=mean pεboundary(C     t     ) └grad( p ) 2 ┘  (19)
 
     Where: 
     boundary (C t ) is the set of pixels on the boundary of the region C t ; and 
     grad(p) is the gradient magnitude of the normalised original image at a pixel p. 
     It can therefore be seen that the energy for a particular threshold t is the mean of the square of the gradient of those pixels that lie on the boundary of the region C t . The processing of step S 91  produces an energy value for each threshold value t that was determined to result in a region C t  containing less than 3000 pixels, i.e. an energy value for each threshold resulting in a region C t  being added to the sets at step S 89 . 
     At step S 92  the values of E(t) are Gaussian smoothed which produces a smoothed plot of the values of energy values E(t) against threshold values t. A suitable value for the Gaussian smoothing function is 0.2, although any suitable value could be used. 
     At step S 93  the values of t at which the Gaussian smoothed plot of the values of E(t) produce a peak are determined and at step S 94  the areas C t  (referred to as regions r) for values of t for which the smoothed plot of E(t) produces a peak are added to a candidate region set R. Values of t at which E(t) is a peak are likely to be where the boundary of C t  separates a blot haemorrhage from its background as the peaks are where the gradient is at a maximum. This is so because the energy function takes as input the gradient at boundary pixels, as can be seen from equation (19). 
     At step S 95  it is determined if there are more candidates in C which have not been processed. If it is determined that there are more candidates in C then processing continues at step S 85  where a new candidate c is selected. If it is determined that all candidates c in C have been processed then at step S 96  the set of regions R is output. 
     Whilst it has been described above that the threshold is incremented by values of 0.1, it will be appreciated that other values of t are possible. For example increasing t by values smaller than 0.1 will give a larger number of areas C t  and therefore a smoother curve of the plot of values of E(t). The value of t may also be beneficially varied based upon the way in which normalization is carried out. Additionally, if it is determined that areas of an image that are possible blot haemorrhages may be larger or smaller than 3000 pixels, different values may be chosen for the threshold of step S 88 . 
     Some of the processing described below benefits from an accurate assessment of the properties of the background local to a particular candidate blot haemorrhage. First, it is necessary to determine a background region relevant to a particular blot haemorrhage.  FIG. 16  shows the processing carried out to determine the relevant background region, which is carried out at step S 55  of  FIG. 11 . At step S 100  a candidate blot haemorrhage pixel c is selected for processing. At step S 101 , a region W centred on the pixel c of dimension 121×121 pixels is determined. A gradient is computed for each pixel in W at step S 102 . A h-minima transform is then applied to the determined gradients at step S 103  to reduce the number of regions generated by subsequent application of a watershed transform as described below. A value of h for application of the h-minima transform is selected such that the number of minima remaining after application of the h-minima transform is between 20 and 60. 
     A watershed transform is then applied to the output of the h-minima transform at step S 104 . The watershed transform divides the area W into m sub-regions. A seed region for the next stage of region growing is then created by taking the union of all sub-regions which intersect the region r (determined at step S 94  of  FIG. 15 ) containing the pixel c at step S 105 . 
     At step S 106  a check is carried out to determine whether the created region is sufficiently large. If this is the case, processing passes to step S 107  where the created region is defined as the background surrounding r. Otherwise, processing continues at step S 108  a further sub-region is added to the created region, the further sub-region being selected from sub-regions which are adjacent the created region, and being selected on the basis that its mean pixel value is most similar to that of the created region. Processing passes from step S 108  to step S 109  where a check is carried out to determine whether adding a further sub-region would result in too large a change in pixel mean or standard deviation. This might be caused if a vessel is included in an added sub-region. If this is the case, processing passes to step S 107 . Otherwise, processing returns to step S 106 . 
     The region created at step S 107  represents a region of background retina surrounding the candidate blot haemorrhage and is denoted B. The region B is used to generate data indicative of the background of the candidate c. 
     A region identified as a candidate blot haemorrhage by the processing of  FIG. 15  may lie on a vessel or on a crossing of a plurality of vessels. In such a case it may be that the region is not a blot haemorrhage. Identifying vessels that are close to an identified candidate blot haemorrhage is therefore desirable. Processing to identify vessels will now be described with reference to  FIG. 17 , and this processing is carried out at step S 56  of  FIG. 11 . 
     Referring to  FIG. 17 , at step S 115  a region r in the set of candidate regions R identified by the processing of  FIG. 15  that has not previously been processed is selected. At step S 116  an area S surrounding the selected region r of the input image is selected. The region r that has been determined to be a candidate blot haemorrhage is removed from the area S for the purposes of further processing. At step S 117  a counter variable q is set to the value 5 and at step S 118  the area S of the image A is tangentially shifted by q pixels. At step S 119  the minimum of: S shifted by q pixels and; the inverse of S shifted by q pixels in the opposite direction; is determined according to equation (20) for all pixels so as to generate an image M τq . 
         M   τq =min(τ q ( S ),−τ −q ( S ))  (20)
 
     At step S 120  it is determined if q has the value 11, which value acts as an upper bound for the counter variable q. If it is determined that q has a value of less than 11 then at step S 121  q is incremented and processing continues at step S 118 . If it is determined at step S 120  that q is equal to 11 then it is determined that the image S has been tangentially shifted by q pixels for q in the range 5 to 11 and at step S 122  an image V is created by taking the maximum at each pixel across the images M τq  for values of q in the range 5 to 11. At step S 123  the image V is thresholded and skeletonised to produce a binary image containing chains of pixels. These chains are split wherever they form junctions so that each chain is a loop or a 2-ended segment. 2-ended segments having one end closer to c than about 0.05 DD (13 pixels) and the other end further than about 0.15 DD (35 pixels) from c are retained as candidate vessel segments at step S 124 , and this set is denoted U seg  with members u seg . Checking that the ends of a segment satisfy these location constraints relative to c increases the chance that the segment is part of a vessel of which the candidate haemorrhage, c, is also a part. All other 2-ended segments and all loops are rejected. 
     Each candidate vessel segment u seg  is classified at step S 125  as vessel or background according to the following features:
         1) Mean width of the candidate vessel segment region;   2) Standard deviation of the width of the candidate vessel segment region;   3) Width of the haemorrhage candidate at an orientation perpendicular to the mean orientation of the candidate vessel segment;   4) The mean of the square of the gradient magnitude along the boundary of the candidate vessel segment region;   5) The mean brightness of the vessel relative to the brightness and variation in brightness in background region B. The background region B is the region of retina surrounding the haemorrhage candidate determined by the processing of  FIG. 16 ;   6) The standard deviation of brightness of the vessel relative to the brightness and variation in brightness in background region B; and   7) The distance that the extrapolated vessel segment passes from the centre of the candidate haemorrhage.       

     Using a training set of candidate vessel segments classified as vessel or background by a human observer, a support vector machine is trained to classify test candidate vessel segments as either vessel or background based on the values evaluated for the above features. The support vector machine outputs a confidence that a candidate vessel is a vessel or background. For each candidate blot haemorrhage the maximum of these confidences is taken for all candidate vessel segments surrounding the candidate blot haemorrhage. 
     At step S 126  it is determined if there are more regions r in R that have not been processed. If it is determined that there are more regions in R then processing continues at step S 115 . 
     Referring now to  FIGS. 18A and 18B , two example candidate regions generated using the processing of  FIG. 15  are shown at four stages of the processing of  FIG. 17 . 
       FIG. 18A  shows a blot haemorrhage  30  and  FIG. 18B  shows an area  31  identified as a candidate blot haemorrhage which is in fact the intersection of a number of vessels and is not a blot haemorrhage. 
     Images (i) in each of  FIGS. 18A and 18B  show candidate blot haemorrhages  30 ,  31  outlined in the original images. Candidate  31  of  FIG. 18B  lies at the join of vessels  32 ,  33  and therefore appears as a dark area similar to a blot haemorrhage. 
     Image (ii) in each of  FIGS. 18A and 18B  shows the result of taking a tangential gradient (step S 118  of  FIG. 17 ) and image (iii) in each of  FIGS. 18A and 18B  shows the image V created at step S 122  of  FIG. 17 . The image (iv) of each of  FIGS. 18A and 18B  shows the original image with identified vessel segments shown as white lines. 
     The location of a candidate blot haemorrhage may be compared to detected vessel segments. Blot haemorrhages are often located on vessels as can be seen in  FIG. 18A . Here it can be seen that a genuine blot haemorrhage lies on a vessel. In this case, a high vessel confidence could cause wrong classification of the blot haemorrhage unless another feature is evaluated that can distinguish between haemorrhages located on vessels such as in  FIG. 18A  and vessel crossings as shown in  FIG. 18B , which may appear similar. Various parameters may be analysed as part of a process referred to as “discontinuity assessment” which allows candidate detections on vessels to be effectively distinguished as haemorrhage or not haemorrhage. 
     Discontinuity assessment is calculated for haemorrhage candidates which have one or more associated candidate vessel segments with a confidence, as calculated at step S 125 , greater than a threshold such as 0.5. Discontinuity assessment can be based upon three factors, calculated using the candidate vessel segments whose confidence, as calculated at step S 125 , is greater than aforementioned threshold. viz: 
     
       
         
           
             
               
                 
                   
                     stronger 
                     ( 
                     i 
                     ) 
                   
                   = 
                   
                     
                       z 
                       1.4 
                       2.8 
                     
                      
                     
                       ( 
                       
                         
                           E 
                           H 
                         
                         / 
                         
                           E 
                           
                             V 
                             i 
                           
                         
                       
                       ) 
                     
                   
                 
               
               
                 
                   ( 
                   21 
                   ) 
                 
               
             
             
               
                 
                   wider 
                   = 
                   
                     
                       z 
                       1.4 
                       2.3 
                     
                      
                     
                       ( 
                       
                         
                           W 
                           C 
                         
                         / 
                         
                           W 
                           
                             i 
                              
                             
                                 
                             
                              
                             n 
                           
                         
                       
                       ) 
                     
                   
                 
               
               
                 
                   ( 
                   22 
                   ) 
                 
               
             
             
               
                 
                   
                     junction 
                     = 
                     
                       max 
                        
                       
                         ( 
                         
                           
                             z 
                             110 
                             140 
                           
                            
                           
                             ( 
                             
                               α 
                               ij 
                             
                             ) 
                           
                         
                         ) 
                       
                     
                   
                    
                   
                     
 
                   
                    
                   
                     where 
                      
                     
                       : 
                     
                   
                    
                   
                     
 
                   
                    
                   
                     
                       
                         z 
                         L 
                         U 
                       
                        
                       
                         ( 
                         x 
                         ) 
                       
                     
                     = 
                     
                       { 
                       
                         
                           
                             0 
                           
                           
                             
                               
                                 if 
                                  
                                 
                                     
                                 
                                  
                                 x 
                               
                               ≤ 
                               L 
                             
                           
                         
                         
                           
                             
                               
                                 ( 
                                 
                                   x 
                                   - 
                                   L 
                                 
                                 ) 
                               
                               
                                 ( 
                                 
                                   U 
                                   - 
                                   L 
                                 
                                 ) 
                               
                             
                           
                           
                             
                               
                                 if 
                                  
                                 
                                     
                                 
                                  
                                 L 
                               
                               &lt; 
                               x 
                               &lt; 
                               U 
                             
                           
                         
                         
                           
                             1 
                           
                           
                             
                               
                                 
                                   if 
                                    
                                   
                                       
                                   
                                    
                                   U 
                                 
                                 ≤ 
                                 x 
                               
                               , 
                             
                           
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   23 
                   ) 
                 
               
             
           
         
       
     
     is a z-function of a type used in fuzzy logic,
 
E H  and E V  are “energies” of the blot haemorrhage candidate and vessel candidate respectively, meaning the mean squared gradient magnitude along the item boundary,
 
W H  is the mean width of the blot haemorrhage candidate,
 
W in  is the diameter of a circle inscribed in the union of all vessel segments after they have been extrapolated towards the blot haemorrhage candidate until the vessel segments intersect each other;
 
α ij  is the intersection angle in degrees between two vessel segments, indexed i and j.
 
A value for discontinuity assessment can be determined using equation (24):
 
     
       
         
           
             
               
                 
                   
                     Discontinuity 
                      
                     
                         
                     
                      
                     assessment 
                   
                   = 
                   
                     max 
                      
                     
                       ( 
                       
                         
                           min 
                            
                           
                             ( 
                             
                               
                                 1 
                                 - 
                                 junction 
                               
                               , 
                               
                                 
                                   min 
                                   i 
                                 
                                  
                                 
                                   ( 
                                   
                                     stronger 
                                      
                                     
                                       ( 
                                       i 
                                       ) 
                                     
                                   
                                   ) 
                                 
                               
                             
                             ) 
                           
                         
                         , 
                         wider 
                       
                       ) 
                     
                   
                 
               
               
                 
                   ( 
                   24 
                   ) 
                 
               
             
           
         
       
     
     Expression 24 takes a value in the range 0 to 1 where 0 represents a low confidence of continuity, meaning the candidate haemorrhage is likely to be part of the detected vessel segment(s) and 1 represents a high confidence of a discontinuity meaning the candidate haemorrhage is likely to be a haemorrhage intersecting a vessel. is calculated to indicate the relation between the width and contrast of the candidate blot haemorrhage and the identified vessels surrounding the candidate blot haemorrhage. 
     The vessel confidence of  FIG. 17  and discontinuity assessment based upon the vessel identification are passed to feature evaluation processing which will now be described with reference to  FIG. 19 . 
     Referring to  FIG. 19 , processing to evaluate the features of a candidate region is shown. The processing is intended to provide data indicating whether a candidate region is likely to be a blot haemorrhage or some other region, for example an intersection of vessels. 
     At step S 130  a candidate region r in the candidate region set R is selected that has not previously been processed. At step  131  a feature vector v r  is determined for the selected candidate region. The feature vector v r  is a vector determined from a number of features as set out in Table 1 below. 
     At step S 132  a check is carried out to determine whether further candidate regions remain to be processed. If this is the case, processing returns to step S 130 . Otherwise processing passes to step S 133  where a candidate vector is selected for processing. At step S 134  a check is carried out to determine whether the candidate vector relates to a candidate region located within 100 pixels of the fovea, which is located using processing of the type described above with reference to  FIG. 10 . 
     If the check of step S 134  is satisfied processing passes to step S 135  where the processed vector is added to a set of vectors associated with candidates within 100 pixels of the located fovea. Otherwise, processing passes to step S 136  where the processed vector is added to a set of vectors associated with candidates located more than 100 pixels from the located fovea. Processing passes from each of steps S 135  and S 136  to step S 137  where a check is carried out to determine whether further candidates remain to be processed. If it is determined that further candidates remain to be processed, processing passes from step S 137  back to step S 133 . 
     
       
         
           
               
               
             
               
                 TABLE 1 
               
               
                   
               
               
                 Feature 
                 Description 
               
               
                   
               
             
            
               
                 Area 
                 The number of pixels in r 
               
               
                 Width 
                 Twice the mean distance of pixels in the skeletonisation (maximal 
               
               
                   
                 morphological thinning) of r from the boundary of r 
               
               
                 Normalised 
                 Int(r)/Con(B) 
               
               
                 Intensity 
                 where: 
               
               
                   
                 Int(r) is the mean intensity of A within r; and 
               
               
                   
                 Con(B) is a measure of contrast in the background surrounding the 
               
               
                   
                 candidate, given by the mean of medium frequency energy present within the 
               
               
                   
                 area generated by the processing of FIG. 16. 
               
               
                 Relative Energy 
                 Mean gradient magnitude of A along the boundary of r divided by the 
               
               
                   
                 minimum of A within r 
               
               
                 Directionality 
                 A histogram of gradient directions, θ, within r is created, with each pixel 
               
               
                   
                 weighted by the local gradient magnitude. The histogram is convolved with a 
               
               
                   
                 filter, 1 + cos(θ). During this convolution, the histogram is treated as periodic. 
               
               
                   
                 Directionality is defined as the standard deviation of the resulting values 
               
               
                   
                 divided by their mean. 
               
               
                 Normalised 
                 (Int(r) − Int(r3))/Con(B) where Int(r3) is the mean intensity in r after 
               
               
                 Relative 
                 morphological erosion by a circle of radius 3. 
               
               
                 Intensity 
               
               
                 Vessel 
                 Maximum of confidences of candidate vessel segments as described with 
               
               
                 confidence 
                 reference to FIG. 17, or zero if no candidate vessel segments were 
               
               
                   
                 detected. 
               
               
                 Discontinuity 
                 Discontinuity assessment of candidate relative to vessel segments as 
               
               
                 Assessment 
                 described above or zero if no candidate vessel segments had a confidence 
               
               
                   
                 higher than the threshold for inclusion within the discontinuity assessment 
               
               
                   
                 evaluation. 
               
               
                   
               
            
           
         
       
     
     When all candidates have been processed in the manner described above, processing passes from step S 137  to step S 138  where vectors associated with candidate regions within 100 pixels of the fovea are processed to identify at most one processed region as the fovea. Candidates which are not identified as the fovea at step S 138 , together with candidates located more than 100 pixels from the expected fovea position, are then input to a support vector machine at step S 139  to be classified as either a blot haemorrhage or background. 
     If the candidate region is within 100 pixels of the fovea, then the blot haemorrhage candidate may in fact be foveal darkening. If a classifier trained to output a confidence of being a fovea or of being a blot haemorrhage returns a higher result for fovea, for one or more haemorrhage candidates, then one of these candidates may be removed from a set of candidate blot haemorrhages. If there is a choice of candidates to be removed then the one nearest to the fovea location, as previously determined, should be removed. The blot haemorrhage candidates should then be classified as blot haemorrhage or background based on their feature vectors. The classification described above may be carried out by a support vector machine trained using a set of candidates generated from a set of training images in which each generated candidate has been hand classified as a fovea, haemorrhage or background by a human observer. 
     A training set of candidate blot haemorrhages are hand-classified as blot haemorrhage or background and the support vector machine is trained using these hand-classified candidates, such that on being presented with a particular feature vector, the support vector machine can effectively differentiate candidate areas which are blot haemorrhages from those which are not. 
     The preceding description has been concerned with the identification of blot haemorrhages. This identification is important, because it is known that the presence of blot haemorrhages on the retina is an indicator of diabetic retinopathy. As such, the techniques described above find application in automated processing of images for the detection of diabetic retinopathy. Blot haemorrhages can also be indicative of other disease conditions. As such, the techniques described above can be used to process images to identify patients suffering from other diseases of which blot haemorrhages are a symptom. 
     It is also known that exudates are indicative of disease states. As such, it is also useful to process retinal images to detect the presence of exudates. 
     Referring now to  FIG. 20 , processing to identify exudates in an image is shown. At step S 150  an image A corresponding to image  2  of  FIG. 1  is input to the computer  5  of  FIG. 1  for processing. At step S 151  the image A is normalised in the same way as previously described with reference to  FIG. 12 . In colour images, exudates are usually most visible in the green colour plane of the image, and as such the processing of  FIG. 12  carried out at step S 151 , and indeed most processing described below, if a colour image is being used, is carried out on the green colour plane. 
     At step S 152  the optic disc is detected. The optic disc is a highly reflective region of the eye and it and the area surrounding it can therefore be falsely detected as exudate. Location of the optic disc is carried out using processing described above with reference to  FIG. 8 . A circular region of the image A of diameter 1.3 DD centred on the optic disc centre is excluded from further analysis. 
     At step S 153  the normalised image is processed to detect candidate exudates as described in further detail below with reference to  FIG. 21 . The processing of step S 153  returns a single pixel in the image A for each detected candidate exudate. At step S 154  the candidate exudates identified at step S 153  are subjected to region growing to determine the region of the image A that is a possible exudate region corresponding to the identified candidate pixel. A suitable procedure for region growing is that described above with reference to  FIG. 15 . 
     At step S 155  watershed region growing is applied as described above with reference to  FIG. 16 . Watershed region growing finds regions of retina that are not vessels or other lesions and these regions are processed to determine some of the features that are evaluated to generate a feature vector at step S 156 , the feature vector is created to include parameters indicative of a candidate exudate. The feature evaluation of step S 156  is described in further detail below. 
     At step S 157  each candidate exudate is processed to determine a confidence that the candidate is exudate, drusen or background. The determination is based upon the feature vector determined at step S 156 . 
     The detection of candidate exudates is now described with reference to  FIG. 21 . Some steps of the processing of  FIG. 21  are very similar to equivalent steps in the processing of  FIG. 13 , and as such are only briefly described. 
     At step S 160  the input image is smoothed in a process similar to that applied at step S 65  of  FIG. 13 . At step S 162  a counter variable n is initialised to a value of 0. At step S 163  a linear structuring element is defined, using the function used at step S 70  in the processing of  FIG. 13 , the function being shown in equation (18). At step S 164  the linear structuring element defined at step S 163  is used in a morphological opening operation of similar form to the operation carried out at step S 71  of  FIG. 13 . At step S 165  a check is carried out to determine whether the counter n has a value of 7. If this is not the case, processing passes from step S 165  to step S 166  where the value of n is incremented before processing continues at step S 163 . When it is determined at step S 165  that the value of n is 7, processing passes to step S 167 . 
     The loop of steps S 163  to S 166  acts in a similar way to that of steps S 70  to S 73  of  FIG. 13  to perform morphological opening with a series of eight structuring elements arranged at different orientations. Each image output from one of the opening operations includes only linear structures extending at a particular orientation. 
     At step S 167  an image D s  is created by subtracting, for each pixel, the maximum value for that pixel across all images M n . As explained with reference to step S 74  of  FIG. 13 , this has the effect of removing linear structures from the image, though, in the case described here with reference to exudate detection, the linear structures removed are brighter than the surrounding retina. 
     Processing passes from step S 167  to step S 168  where a check is carried out to determine whether the value of s is equal to eight. If this is not the case, processing passes to step S 169  where the value of s is incremented, before processing continues at step S 170  where the image is scaled, relative to the original image, by a scaling factor based upon s, more particularly the scaling factor 2 s-1  described with reference to  FIG. 13 . Processing passes from step S 170  to step S 162 . 
     When it is determined at step S 168  that the value of s is equal to 8, processing passes to step S 171 . At step S 171 , a check is carried out for a particular pixel to determine whether the maximum value for that pixel across all images D s  is greater than a threshold, determined as described below. If this is the case, a candidate region associated with the pixel is considered to be candidate exudate at step S 172 . Otherwise, the candidate region is not considered to be a candidate exudate at step S 173 . 
     The threshold applied at step S 171  is selected firstly by fitting a gamma-distribution to the distribution of heights of the regional maxima in D s . The threshold is placed at the point where the cumulative fitted distribution (its integral from −∞ to the point in question, with the integral of the whole distribution being 1) is 1-5/n, where n is the number of maxima in D s . Only those maxima in D s  which are less than this threshold are retained. 
     Referring to  FIG. 22 , processing to evaluate the features of a candidate region is shown. The processing is intended to provide data indicating whether a candidate region is likely to be exudate, drusen or background. Some steps of the processing of  FIG. 22  are very similar to equivalent steps in the processing of  FIG. 19 , and as such are only briefly described. 
     At step S 175  a candidate region r in the candidate region set R is selected that has not previously been processed. At step S 176  a feature vector v r  is determined for the selected candidate region. The feature vector v r  is a vector determined from a number of features as set out in Table 2 below. 
     
       
         
           
               
               
             
               
                 TABLE 2 
               
               
                   
               
               
                 Feature 
                 Description 
               
               
                   
               
             
            
               
                 Area 
                 The number of pixels in r 
               
               
                 distance from 
                 The distance of the candidate from the nearest Microaneurysm. 
               
               
                 MA 
                 Microaneurysm detection is described below with reference to FIG. 13. 
               
               
                 Normalised 
                 (L r  − L bg )/C bg   
               
               
                 Luminosity 
                 where L r  is the mean luminosity in the region r of the normalised image 
               
               
                   
                 generated by the processing of FIG. 12; 
               
               
                   
                 L bg  is the mean luminosity in the local background to the candidate 
               
               
                   
                 determined within the area generated by the processing of FIG. 16; and 
               
               
                   
                 C bg  is the mean contrast in the local background to the candidate determined 
               
               
                   
                 within the same area. 
               
               
                 Normalised sd 
                 sd(L r )/C bg   
               
               
                 of Luminosity 
                 where L r  and C bg  are as described above. 
               
               
                 Normalised 
                 Calculated as the mean gradient magnitude of the region r along the 
               
               
                 Boundary 
                 boundary of r divided by C bg   
               
               
                 Gradient 
               
               
                 Spread 
                 The spread of the region r evaluated as ({square root over (N r )}/d − 3{square root over (π)})/2 
               
               
                   
                 where N r  is the number of pixels in r; and 
               
               
                   
                 d is the mean distance of pixels in r from its boundary. The spread value has 
               
               
                   
                 a minimum of 0 for a circle. 
               
               
                 Standardised 
                 A good quality well-exposed retinal image is chosen as a standard. 
               
               
                 Colour Features 
                 Histogram standardised planes are generated by applying a strictly 
               
               
                   
                 increasing transformation to the red and green colour planes of I so that 
               
               
                   
                 each result has a histogram which is similar to that of the corresponding 
               
               
                   
                 plane of the standard image. The mean is taken of the standardised red and 
               
               
                   
                 green planes over r. 
               
               
                   
               
            
           
         
       
     
     At step S 177  a check is carried out to determine whether further candidate regions remain to be processed and if this is the case, processing returns to step S 176 . Otherwise processing passes to step S 178  where a candidate vector is selected for processing. 
     A basic support vector machine is able to perform binary classification. To allow classification as either exudate, drusen or background, each of the classes are compared to each of the other classes using three one-against-one support vector machines and the mean is taken of the results. At step S 179  the selected vector is processed by a support vector machine to classify the candidate as either exudate or drusen. At step S 180  the selected vector is processed by a second support vector machine to classify the candidate as either exudate or background and at step S 181  the selected vector is processed by a third support vector machine to classify the candidate as either drusen or background. Each support vector machine outputs a likelihood that a candidate is each of the two categories that the support vector machine is trained to assess. The likelihood for both categories sums to 1. At step S 181  the mean of the likelihoods output from the three support vector machines for each class is taken. It will be appreciated that the resulting likelihoods calculated by taking the mean in the manner described above for the three categories will also sum to 1. 
     At step S 182  a check is performed to determine if there are more candidates to be evaluated. If it is determined that there are more candidates to be evaluated then the processing of steps S 178  to S 182  is repeated. Otherwise at step S 183  the processing of  FIG. 22  ends. 
     A training set of candidate exudates are hand-classified as exudate, drusen or background and each support vector machine is trained upon these hand-classified candidates, such that on being presented with a particular feature vector, each support vector machine can effectively differentiate candidate areas which the particular support vector machine is intended to classify. 
     Referring now to  FIG. 13 , processing carried out to detect microaneurysms is described. Detection of microaneurysms is required in order to evaluate the feature “distance from microaneurysm” shown in Table 2. 
     Candidate microaneurysms are located using a method similar to that of  FIG. 13 , although the input image is processed only at a single scale. That is, the processing of  FIG. 13  is performed without the loop provided by step S 74 , and consequently without repeated scaling of the image as carried out at step S 76 . A set of candidate microaneurysms is created, however, as discussed with reference to steps S 77  and S 78  above. At step S 77 , the threshold used to determine whether a processed region is a candidate microaneurysm can suitably be selected to be 5 times the 95th percentile of pixels in D. 
     Each candidate microaneurysm, represented by a respective pixel, is subjected to region growing as described with reference to  FIG. 15  above so as to create a candidate area for each microaneurysm. Watershed region growing, as described above with reference to  FIG. 16  is also carried out to allow characteristics of the background of a candidate microaneurysm to be determined. More particularly, an estimate of background contrast: the standard deviation of pixels in the normalised image after high pass filtering within the region obtained from watershed retinal region growing can be determined and denoted BC. 
     A paraboloid is then fitted to the 2-dimensional region generated by the processing of  FIG. 15 . From the fitted paraboloid, the major- and minor-axis lengths are calculated as well as the eccentricity of the microaneurysm candidate. 
     Features used to determine whether a particular candidate microaneurysm is in fact a microaneurysm may include:
         1. The number of peaks in energy function E, where the energy function has a form similar to equation (19) above;   2. Major and minor axis lengths determined as described above;   3. The sharpness of the fitted paraboloid (or alternatively the size of the fitted paraboloid at a constant depth relative to its apex can be used since this is inversely proportional to the sharpness of the paraboloid);   4. Depth (relative intensity) of the candidate microaneurysm using the original image and the background intensity estimated during normalisation;   5. Depth of the candidate microaneurysm using the normalised image and the fitted paraboloid divided by BC;   6. Energy of the candidate microaneurysm, i.e. the mean squared gradient magnitude around the candidate boundary divided by BC.   7. The depth of the candidate microaneurysm normalised by its size (depth divided by geometric mean of axis lengths) divided by BC.   8. The energy of the candidate microaneurysm normalised by the square root of its depth divided by BC.       

     Using a training set, a K-Nearest Neighbour (KNN)-classifier is used to classify candidates. A distance metric is evaluated between a feature vector to be tested and each of the feature vectors evaluated for a training set in which each of the associated candidate microaneurysms was hand-annotated as microaneurysm or not microaneurysm. The distance metric can be evaluated, for example as the sum of the squares of differences between the test and training features. A set is determined consisting of the K nearest, based on the distance metric, training candidate feature vectors to the test candidate feature vector. A candidate is considered to be a microaneurysm if L or more members of this set are true microaneurysms. For example, a candidate microaneurysm would be considered to be a true microaneurysm for L=5 and K=15 meaning 5 out of 15 nearest neighbours are true microaneurysms. 
     The method of detecting blot haemorrhages described above has been tested on 10,846 images. The images had been previously hand classified to identify blot haemorrhages present as follows: greater than or equal to four blot haemorrhages in both hemifields in 70 images; greater than or equal to four blot haemorrhages in either hemifield in 164 images; macular blot haemorrhages in 193 images; blot haemorrhages in both hemifields in 214 images; and blot haemorrhages in either hemifield in 763 images. 
     Receiver Operator Characteristic (ROC) curves for each of these categories are displayed in  FIG. 23 . A line  101  shows data obtained from images including four or more blot haemorrhages in both hemifields. A line  102  shows data obtained from images having four or more blot haemorrhages in either hemifield. A line  103  shows data obtained from images including macular blot haemorrhages. A line  104  shows data from includes including blot haemorrhages in both hemifields, while a line  105  shows data from images including blot haemorrhages in either hemifield. 
     Since the images with blot haemorrhages were drawn from a larger population than images without blot haemorrhages, data was rated to adjust to the prevalence of blot haemorrhages in the screened population of images, estimated to be 3.2%. High sensitivity and specificity are attained for detection of greater than or equal to 4 blot haemorrhages in both hemifields (98.6% and 95.5% respectively) and greater than or equal to four blot haemorrhages in either hemifield (91.6% and 93.9% respectively). 
     The method of detecting exudates as described above has been tested on a set of 13,219 images. Images had been previously classified manually for the presence of exudates and drusen as follows: 300 with exudates less than or equal to 2 DD from the fovea, of which 199 had exudates less than or equal to 1 DD from the fovea; 842 images with drusen; 64 images with cotton-wool spots; 857 images with other detectable bright spots. 13.4% (1825) of the images with exudates contained one of the other categories of bright objects. 
       FIG. 24  shows ROC curves for exudate detection less than or equal to 2 DD from the fovea) ( FIG. 24A ) and less than or equal to 1 DD from the fovea ( FIG. 24B ). Images with referable or observable exudates (less than or equal to 2 DD from the fovea) were recognised at sensitivity 95.0% and specificity 84.6% and images with referable exudates (less than or equal to 1 DD from the fovea) were recognised at sensitivity 94.5% and specificity 84.3%. 
     Although it is necessary to check the performance of the automated system by comparison with a human observer, it should be recognised that opinions confirming the disease content of retinal images can differ substantially. In studies comparing automated exudate detection with human expert detection, a retinal specialist attained 90% sensitivity and 98% specificity compared to a reference standard and a retinal specialist obtained 53% sensitivity and 99% specificity compared to a general ophthalmologist. The latter of these results is close to the ROC curve in  FIG. 24 . 
     The methods described above can be applied to retinal images to enable effective detection of blot haemorrhages and exudates. It is known, as indicated above, that the presence of blot haemorrhages and exudates in retinal images is indicative of various disease. Thus, the methods described herein can be effectively employed in the screening of retinal images by an automated, computer-based process. That is, a retinal image may be input to a computer arranged to carry out the methods described herein so as to detect the presence of blot haemorrhages and exudates within the image. Data indicating the occurrence of blot haemorrhages and exudates can then be further processed to automatically provide indications of relevant disease, in particular indications of diabetic retinopathy or age-related macular degeneration. 
       FIG. 25  is a schematic illustration of a suitable arrangement for providing indications of whether a particular image includes indicators of disease. An image  200  is input to a computer  201 . The image  200  may be captured by a digital imaging device (e.g. a camera) and provided directly to the computer  201  by an appropriate connection. The computer  201  processes the image  200  and generates output data  202  (which may be displayed on a display screen, or provided in printed form in some embodiments). The computer  201  carries out various processing. In particular, a blot haemorrhage detection process  203  is arranged to process the image in the manner described above with reference to  FIG. 11  to determine whether the image includes blot haemorrhages. An exudate detection process  204  is arranged to process the image in the manner described above with reference to  FIG. 20  to identify exudates within the image  200 . Data generated by the blot haemorrhage detection process  203  and the exudate detection process  204  is input to a disease determination process  205  which is arranged to generate the output data  202  discussed above. 
     The computer  201  can conveniently be a desktop computer of conventional type comprising a memory arranged to store the image  200 , the blot haemorrhage detection process  203 , the exudates detection process  204  and the disease determination process  205 . The various processes can be executed by a suitable microprocessor provided by the computer  201 . The computer  201  may further comprise input devices (e.g. a keyboard and mouse) and output devices (e.g. a display screen and printer). 
       FIG. 26  shows a schematic illustration of an alternative arrangement to that of  FIG. 25  for providing indications of whether images taken from a particular patient includes indicators of disease, based upon an image of the right eye of the patient  206  and an image of the left eye of the patient  207 . Both the right and left eye images  206 ,  207  are input to the computer  208 . The images  206 ,  207  may again be captured by a digital imaging device (e.g. a camera) and provided to the computer  208  by an appropriate connection such as a connection from a patient management system or a connection directly from the camera. The computer  208  processes the images  206 ,  207  and generates output data  209  which may be displayed on a display screen, or provided in printed form. The computer  208  carries out various processing on different portions of the images. 
     In particular, for each image  206 ,  207 , a microaneurysm detection process  210  is arranged to detect microaneurysms in the whole of each of the images in the manner described above with reference to  FIG. 13 . More specifically, the processing of  FIG. 13  is performed without the loop provided by step S 74 , and consequently without repeated scaling of the image as carried out at step S 76 . 
     A blot haemorrhage detection process  211  is arranged to process an area within 1 DD of the centre point of the fovea of each of the images in the manner described above with reference to  FIG. 11  (i.e. a circular area centred on the centre of the fovea having a radius of 1 DD). The fovea can be detected in each of the images using the processing of  FIG. 10  or using any other suitable method for detection of the fovea. 
     An exudate detection process  212  is arranged to process the area within 1 DD of the centre point of the fovea described above in the manner described previously with reference to  FIG. 20  so as to detect the presence of exudates. A further exudate detection process  213  is arranged to process a larger area within 2 DD of the centre point of the fovea of each of the images, again in the manner described above with reference to  FIG. 20 . 
     It will be appreciated that the further exudate detection process  213  can use the data from the exudate detection process  212  and search only the area formed from the circular area with a radius of 2 DD centred on the centre of the fovea not searched by the exudate detection process  212  which operates on a circular area of radius 1 DD centred on the centre of the fovea. It will be appreciated that reducing the area of the image that is processed to detect blot haemorrhages and exudates lessens required processing and can decrease the time taken to process each pair of images in the described manner. 
     Data generated by each of detection processes  210  to  213  for each of the images is passed to a disease determination process  214 . Each of the processes  210  to  214  may be initiated in any convenient way. For example the processes may be initiated through a patient management system or through interaction with the camera used to acquire the images. The processing of disease determination process  214  will now be described with reference to  FIG. 27 . 
     At step S 200  a set of values L Left  and a value MA Left  are generated from the data generated from the image of the left eye  207  and at step S 201  a set of values L Right  and a value MA Right  are generated from the data generated from the image of the right eye  206 . Processing to generate the sets of values L Left  and L Right  is described in further detail below with reference to  FIG. 28 . At step S 202  a value D is generated from the sets of values L Left  and L light  and the values MA Left  and MA Right  according to equation (25) below: 
     
       
         
           
             
               
                 
                   D 
                   = 
                   
                     
                       2 
                        
                       
                           
                       
                        
                       
                         max 
                          
                         
                           ( 
                           
                             
                               MA 
                               Left 
                             
                             , 
                             
                               MA 
                               Right 
                             
                           
                           ) 
                         
                       
                     
                     + 
                     
                       
                         ∑ 
                         
                           i 
                           = 
                           1 
                         
                         3 
                       
                        
                       
                         [ 
                         
                           
                             α 
                             i 
                           
                            
                           
                             max 
                              
                             
                               ( 
                               
                                 
                                   L 
                                   i 
                                   Left 
                                 
                                 , 
                                 
                                   L 
                                   i 
                                   Right 
                                 
                               
                               ) 
                             
                           
                         
                         ] 
                       
                     
                   
                 
               
               
                 
                   ( 
                   25 
                   ) 
                 
               
             
           
         
       
     
     where MA Left  and MA Right  are the number of microaneurysms detected by the microaneurysm detection process  210  of  FIG. 26  in the left and right eyes respectively, L i   Left  and L i   Right  are the ith values in the sets L Left  and L Right  respectively, max(L i   Left ,L i   Right ) is the maximum of L i   Left  and L i   Right , and α i  is a weight associated with the ith value in each of the sets L Left  and L Right . 
     Processing to generate the sets L Left  and L Right  will now be described with reference to  FIG. 28 . The processing of  FIG. 28  is repeated for each of the left and right patient eye images  206  and  207 , to generate the set L Left  and L Right  respectively. 
     At step S 207  a set BH is received from the blot haemorrhage detection process  211  of  FIG. 26 . Each element of the set BH represents a candidate blot haemorrhage within 1 DD of the centre of the fovea, and each candidate blot haemorrhage has an associated confidence value indicating a confidence that the candidate represents a blot haemorrhage. At step S 208  a value L 1  is determined by summing the three largest confidences associated with members of the set BH. 
     At step S 209  a set Ex 1  is received from the exudate detection process  212  of  FIG. 26 . Each element of the set Ex 1  represents a candidate exudate within 1 DD of the centre of the fovea, and each candidate exudate has an associated confidence indicating a confidence that the candidate represents an exudate. At step S 210  a value L 2  is determined by summing the three largest confidences associated with members of the set Ex 1 . At step S 211  a set Ex 2  is received from the exudate detection process  213  of  FIG. 26 . Each element of the set Ex 2  represents a candidate exudate within 2 DD from the centre of the fovea. At step S 212  a value L 3  is generated from the set Ex 2  by summing the two largest confidences associated with members of the set Ex 2 . 
     Whilst it has been described that the disease determination processing  214  of  FIG. 26  takes input from the detection processes described previously with reference to  FIGS. 10 ,  11 ,  20  and  23 , it will be appreciated that the disease determination process depends only upon the sets of values input to it, and any suitable detection methods may be used in detection processes  210  to  213  of  FIG. 26 . 
     The weights α i  of equation (25) can be determined using a training set of manually classified pairs of patient eye images. Each pair of images from a patient is classified as either indicating disease or not indicating disease. Each patient right eye image and left eye image in the training set has associated values corresponding to the sets of values L Right  and L Left . A linear classifier is then trained to predict disease by optimising the weights α i  applied to the elements of the sets L Left  and L Right  according to equation (25). More particularly, values for the weights a, are generated by choosing those values that maximise the area under a Receiver Operator Characteristic (ROC) curve. The ROC curve plots sensitivity against 1-specificity for a test designed to detect referable retinopathy (the specificity and sensitivity being determined with reference to data provided by a clinically trained observer) as a threshold applied to D in equation (25) is varied. Other suitable methods for determining weights can of course be used. 
       FIG. 29A  shows how the area under respective ROC curves varies as values for weights α i  are varied at a coarse scale between values of 0 and 10.  FIG. 29B  shows how the area under respective ROC curves varies as values for the weights α i  are varied between 0 and 2 at a finer scale. Each line shows variance of an area under an ROC curve as a respective weight α i  is varied while all other weights are fixed at their respective optimum values. The areas are calculated based upon a training set described in further detail below. It can be seen that the scale on the y-axis is such that the area under the ROC curve varies relatively little (little more than 0.01) as the values of the weights vary between 0 and 10. 
     A line  215  in  FIG. 29A  shows how the area under a ROC curve varies as the value of a weight α 4 , associated with a further value L 4  associated with each training image, indicative of a confidence value for four or more blot haemorrhages in either hemifield of the image of the eye, is varied. It can be seen from  FIG. 29A  that the area under the ROC curve decreases sharply as the value of α 4  increases from a value of 0 to 10. The line  215  in  FIG. 29B  shows that as the values of the weight α 4  are varied between 0 and 2 the area under the ROC curve drops sharply as the weight increases from a value of 0 indicating that the optimal value for the weight α 4  is 0. That is, the value L 4  need not be included in equation (25). 
     A line  216  in  FIG. 29A  shows how the area under a ROC curve varies as the value of the weight α 3  described above with reference to equation (25) is varied between 0 and 10. It can be seen from  FIG. 29A  that the area under the ROC curve increases as the weight α 3  increases from a value of 0 to a value of approximately 4. The area under the ROC curve remains substantially constant while the weight α 3  has value between 4 and 7 and then the area under the ROC curve reduces for values of α 3  greater than 7. While the line  216  is shown in  FIG. 29B , the range of values of α 3  shown in  FIG. 29B  are not of relevance to its determination. 
     A line  217  in  FIG. 29A  shows how the area under a ROC curve varies as the value of the weight α 2  described above with reference to equation (25) varies. It can be seen from  FIG. 29A  that the area under the ROC curve remains substantially constant for values of the weight α 2  between 0 and 7 and then the area under the ROC curve reduces for values greater than 7. While the line  217  is shown in  FIG. 29B , the range of values of α 2  shown in  FIG. 29B  are not of relevance to the determination of the value of the weight α 2 . 
     A line  218  in  FIG. 29B  shows how the area under a ROC curve varies as the value of the weight α 1  described above with reference to equation (25) varies. It can be seen from  FIG. 29B  that the area under the ROC curve remains substantially constant for values of the weight α 1  from 0 to 0.8 and reduces for values greater than 0.8. While the line  218  is shown in  FIG. 29A , the range of values of α i  shown in  FIG. 29B  are of greater relevance to its determination. 
     Whilst  FIGS. 29A and 29B  show variance of the area under ROC curves for each weight, in one embodiment the actual values of the weights α i  were selected based upon a local maximum for the area under a ROC curve based as all weights are varied together. The values were determined based upon a training set of images taken from 200 patients randomly selected who had been manually classified as having disease and 400 patients manually classified as not having disease. The values of α i  determined on the training set are shown in Table 3 below, together with the portion of the photograph used to determine each respective confidence. 
     The values of α i  applied to the values L 1  to L 3  shown in Table 3 are intended to be generally applicable and not be specific to the particular arrangement of equipment used in the study set out above. However, use of different detection algorithms, or even particular arrangements of equipment, could mean different values of α i  applied to the values L 1  to L 3  are more effective than those shown. In this case, suitable values of α i  to be applied to the values L 1  to L 3  can be determined based upon images obtained using the particular arrangement of equipment and processed using the particular detection algorithms. 
     
       
         
           
               
               
               
               
             
               
                   
                 TABLE 3 
               
               
                   
                   
               
               
                   
                   
                 Portion of photograph 
                 Weight 
               
               
                   
                 Image-based lesion measure 
                 used 
                 (α i ) 
               
               
                   
                   
               
             
            
               
                   
               
            
           
           
               
               
               
               
            
               
                 L 1   
                 Confidence value for blot 
                 ≦1 DD from the centre 
                 0.4 
               
               
                   
                 haemorrhage ≦1 DD from the 
                 of the fovea 
               
               
                   
                 centre of the fovea 
               
               
                 L 2   
                 Confidence value for exudates 
                 ≦1 DD from the centre 
                 6 
               
               
                   
                 ≦1 DD from the centre of the 
                 of the fovea 
               
               
                   
                 fovea 
               
               
                 L 3   
                 Confidence value for exudates 
                 ≦2 DD from the centre 
                 5 
               
               
                   
                 ≦2 DD from the centre of the 
                 of the fovea 
               
               
                   
                 fovea 
               
               
                 L 4   
                 Confidence value for four or 
                 Whole photograph 
                 0 
               
               
                   
                 more blot haemorrhages in 
               
               
                   
                 superior or inferior hemifield 
               
               
                   
               
            
           
         
       
     
     Given that the weight α 4  associated with the value L 4  is zero, the value of L 4  is not included in equation (25) as can be seen by the upper limit on the summation. 
     The value D produced at step S 202  of  FIG. 27  is thresholded to give a Boolean disease/no disease determination. The value of the threshold applied to the value D may be chosen to give a required specificity for detection of observable/referable retinopathy. In particular, the threshold may be chosen based upon a training set of images as classified by a clinically trained observer. 
     The further value L 5  associated with each training image, indicative of a confidence value for four or more blot haemorrhages in either hemifield of the image of the eye, was included in sets L Left  and L Right  for the purposes of verifying that blot haemorrhages need only be located within 1 DD of the centre of the fovea (it will be recalled that the value L 2  is associated with blot haemorrhages within 1 DD of the centre of the fovea). As set out above with reference to  FIGS. 29A and 29B , the area under an ROC curve decreases sharply as the value of α 5  associated with the value L 5  increases. The optimal value of the weight α 5  associated with the value L 5  was therefore determined to be zero, indicating that a lesion configuration of four or more blot haemorrhages in either or both hemifield does not contribute to the overall assessment of disease. This confirms that searching the whole of the image for blot haemorrhages is not necessary and as such searching for blot haemorrhages can be limited to an area of each image with a diameter of 1 DD centred on the fovea, as has been described in the processing set out above with reference to  FIGS. 27 and 28 , and as is represented by the data L 2 . 
     The automated method for determining an indication of disease described with reference to  FIG. 26  above was tested on a set consisting of 1253 patients manually graded as having observable/referable retinopathy and 6333 patients manually graded as having mild or no retinopathy. The median age of the test set was 65 years with inter-quartile range 19. Reference standard grading was performed by a clinical research fellow using the Scottish Diabetic Retinopathy Grading Scheme 2005. 
       FIG. 30  shows a Receiver Operator Characteristic curve (ROC) for the detection of lesions of observable/referable retinopathy using the process described above. The data was weighted to correct for the higher prevalence of observable/referable retinopathy in images in the test set (12.9%) than in the screened population, found to be 3.2% in an earlier study.  FIG. 31  shows the number of patients referred and not referred for manual grading by the automated system for each grade of observable/referable retinopathy. Percentages are included with 95% confidence intervals (CI). 
     Referring to  FIG. 31 , headings enclosed by braces  220   a  and  220   b  indicate classifications assigned to images by manual classifiers. Images classified as “No retinopathy” and “Mild retinopathy” are considered to not warrant referral to a specialist ophthalmologist. Images classified as “Technical failure” indicate that the image quality is not high enough or the image could not be processed for any other reason. These images should be referred to manual graders by an automated system. Column  221  indicates a total of all other image classifications which are considered to be observable/referable images, that is the total number of images classified as one of the categories enclosed by brace  220   b . Column  222  indicates a total of those images considered to indicate referable retinopathy. Rows  223  and  224  indicate the numbers of images not referred and referred respectively by the automated system described above with reference to  FIG. 26 . Row  225  indicates the percentage of images in each manual image classification referred by the automated system. 
     In a known manual grading system patient images are graded twice by disease/no disease graders. The first grading passes images with any sign of disease to a second grading. Workload reduction with automation, based on a 37.9% manual first grading referral rate and a prevalence of referable cases in the screened population of 4.9%, was calculated for the test set according to equation (26) below. 
     
       
         
           
             
               
                 
                   100 
                   - 
                   
                     100 
                     * 
                     
                       
                         ( 
                         
                           
                             4.9 
                             × 
                             sensitivity 
                           
                           + 
                           
                             
                               ( 
                               
                                 100 
                                 - 
                                 4.9 
                               
                               ) 
                             
                             × 
                             
                               ( 
                               
                                 1 
                                 - 
                                 specificity 
                               
                               ) 
                             
                           
                         
                         ) 
                       
                       
                         [ 
                         
                           100 
                           + 
                           37.9 
                         
                         ] 
                       
                     
                   
                 
               
               
                 
                   ( 
                   26 
                   ) 
                 
               
             
           
         
       
     
     The value sensitivity is a value between 0 and 1 indicating the proportion of referable images that the automatic system refers. The value (1−specificity) is a value between 0 and 1 indicating the proportion of non-referable images that the automatic system refers. 
     The numerator of the fraction of equation (26) is the number of manual gradings per 100 patients when automation is used, given that automation replaces the manual first grading step. The value 4.9× sensitivity gives the number of referable images referred by the automatic system. The value (100−4.9)×(1−specificity) gives the number of non-referable images referred by the automated system. The denominator of the fraction [100+37.9] is the number of manual gradings per 100 patients in a fully manual system. This assumes that in a completely manual system, all patients will be first graded and 37.9% will be second graded. 
     The workload reduction was calculated to be 59.1%. Six patients classified by the reference graders as having proliferative retinopathy were missed by the automated system. Three of these were thought to have new vessels at the disc, all of which had microaneurysms and one had haemorrhage within one disc diameter of the centre of the fovea. The remaining three were thought to have new vessels elsewhere, of which one had microaneurysms. The average time to process each image was 320 seconds on a 3 GHz processor. This would allow up to 390,000 images to be processed annually on a computer unit with 4 processor cores. 
     Although specific embodiments of the invention have been described above, it will be appreciated that various modifications can be made to the described embodiments without departing from the spirit and scope of the present invention. That is, the described embodiments are to be considered in all respects exemplary and non-limiting. In particular, where a particular form has been described for particular processing, it will be appreciated that such processing may be carried out in any suitable form arranged to provide suitable output data.