to do functional neuroanatomy in the brain is to relate a discrete area , or network of areas , to a specific function , or set of functions .
the strongest evidence for such a relation is the observation of disruption of a function following disruption of its putative anatomical substrate .
unfortunately , such evidence is difficult to obtain in the human brain because disrupting its activity can be done experimentally only transiently , and only for accessible regions of cortex .
we must therefore rely on data derived from patients with focal brain lesions of natural , or incidental surgical , causes .
now human lesion data is difficult to use for the purposes of functional neuroanatomy for two conflicting reasons .
first , delineating the precise extent of the lesion  at a resolution commensurate with the underlying anatomical architecture  has hitherto been done manually , by trained operators , making the process extremely time - consuming and susceptible to operator bias .
second , the resolving power of a lesion study is generally limited not by the physical resolution of the images but by the resolution of the lesion sampling of the brain : effectively the anatomical scale at which the effects of the absence or presence of damage can be reliably determined .
for example , if within a set of patients under study whenever a given voxel is hit a cluster of other voxels are always also hit the resolution of the resultant lesion map is not limited by the dimensions of the voxel but by the size of the cluster of invariantly affected other voxels .
this limit depends not only on variations in the frequency of damage to locations across the brain but also on the multivariate pattern of damage in the population of lesions , a factor that is hard to quantify owing to the likely complexity of what is a very high - dimensional multivariate distribution ( nachev and husain , 2007 ) .
a further complication of lesion studies is the dynamic nature of the consequences of focal injury on the operation of what is inevitably a distributed , plastic network .
acutely , areas remote from the site of injury may be transiently affected in ways that do not necessarily reflect the functional contribution of the target .
chronically , remote reorganisation may abnormally compensate for a deficit , camouflaging the target 's true role in the normal state . to obtain a synoptic picture of the role of a given area we therefore need both acute and chronic lesion studies , with image processing methodology optimised for each .
to realize in practice the power lesion - mapping has in theory we thus need methodology that permits much larger datasets to be generated ; inevitably , in a fully - automated manner .
this requires the development of unsupervised algorithms for the two critical steps in the processing of lesion images : distinguishing damaged from normal brain ( lesion segmentation ) and determining the anatomical labels of the damaged areas ( lesion registration ) .
although a number of satisfactory algorithms exist for the latter ( crinion et  al . , 2007 ; andersen et  al . , 2010 ; nachev et  al . , 2008 ) no comprehensive solution exists for the former
first , we show how a simple recently - described general measure of anomaly can be used to perform lesion segmentation theoretically in any imaging modality where inter - subject registration to a set of reference normal images is possible .
second , we describe and proceed to evaluate an algorithm based on this approach that is optimised for the characterisation of acute ischaemic lesions as shown by diffusion - weighted magnetic resonance imaging ( mri ) ( dwi ) .
third , we apply the new method to acute ischaemic lesions involving the occipital lobe so as to generate a map of the patterns of damage to the region , facilitating predictions about the resolvability with lesion - mapping of specific subareas within this region . in describing the general approach and our specific application
we need to make a set of general points about lesion segmentation agnostically of the lesion type and imaging modality , and a set of specific points pertinent to acute stroke . taking the general points first , any comprehensive method for lesion segmentation has to deal with five fundamental problems .
first , for any given imaging modality , the signal at any specific point in the brain will usually vary from one normal individual to another in a way that is difficult to parameterise : the population distribution is often not only not gaussian , but multimodal .
our method is therefore non - parametric ( lao et  al . , 2008 ) .
second , deciding whether or not a region is abnormal often depends on the signal not just in one imaging sequence but several different ones : where an abnormality is not replicated across more than one type of sequence it may merely reflect noise or artefact .
our method is therefore potentially multispectral ( prastawa et  al . , 2004 ) .
third , although the signal properties of normal tissue may be definable , they are often not for lesions , simply because it is in the very nature of pathology to be heterogeneous in signal .
our method is therefore agnostic of the specific properties of the lesion signal : it identifies everything that is anomalous in relation to the normal reference ( prastawa et  al . , 2004 ;
fourth , whether or not the signal at any given locus is interpreted as normal or damaged often depends on the signal in its immediate anatomical vicinity .
our method therefore incorporates local information , in a high - dimensional manner , when determining the anomaly of each point in the brain .
fifth , the optimal properties of an image on which to perform lesion segmentation are opposite to those of an image on which to perform lesion registration : this is so because in the former normal tissue contrast interferes with the lesion contrast one needs to distinguish normal from damaged brain , whereas in the latter lesion contrast interferes with the normal tissue contrast one needs to determine the anatomical labels of the lesion .
our method therefore uses different imaging sequences for each task : one optimised for lesion segmentation and another optimised for lesion registration . since
clinical scans invariably use at least two different sequences this does not limit the application .
the core of our method is a simple measure of the anomaly of an unknown test datum in relation to a reference set of data already known to be normal . to determine the anomaly of a single datum
one may simply compare it to the k instances within the reference set that resemble it most closely : its k nearest neighbours ( cover and hart , 1967 ) . where the datum is a single value , a scalar ,
this is simply a matter of finding the k points that are closest to it on a linear scale . where the datum has n variables describing it , a vector , this is some distance measure in n dimensional space
 most simply the euclidean  of the datum to its k nearest neighbours . to derive a unitary measure
, one can take the mean of the distances of the test datum to each of the k nearest neighbours , a measure known as gamma (  ) ( harmeling et  al . , 2006 ) .
the attraction of gamma is that it is indifferent to the shape and number of modes of the reference population distribution , and it is determined only by the reference set , requiring no prior knowledge of the properties of anomalous data .
the problem with gamma , however , is that it is sensitive to variations in the density of the reference population , that is , the degree of similarity between normal points within the reference set .
1a , where a heterogeneous synthetic two - dimensional dataset is labelled by the gamma score of each point : the diameter of each point is drawn in proportion to the score .
it is easy to see that members of the smaller , denser cluster will generally have a smaller value of gamma than members of the larger , sparser cluster even though each belongs equally strongly to its respective cluster .
thus , if anomaly were to be determined by a fixed threshold of gamma , either too many members of the sparse cluster will be labelled as anomalous or too few of the denser one .
it is thus easy to see that the score is confounded by local differences in density .
so for gamma to be useful we need to find a way of correcting for such differences . as originally proposed in the context of computer network intrusion detection ( rieck and laskov , 2006 ) , an elegant way of doing this is to correct it by the mean of the distances between the members of the  clique  of k nearest neighbours , effectively a measure of local density .
thus , where the k nearest neighbours are sparsely distributed , gamma will have to be larger to be judged to be anomalous , and vice versa .
the same dataset relabelled by this corrected score , zeta (  ) ( rieck and laskov , 2006 ) , is shown in fig . 
we can now see that differences in local density no longer affect the score to the same degree . in short
, the zeta score embodies the simple intuition that you are an outsider if those closest to you are closer to each other than they are to you . to apply zeta to the task of lesion segmentation
, we calculate the zeta score independently for each voxel within a test image that has been brought into registration with a large set of images known to be free of lesions .
the test datum is thus the signal at each voxel , and the reference data from which we choose the k nearest neighbours is the signal in each of the voxels at the same anatomical location across the reference set .
rather than using just the signal at a given voxel , a scalar value , we can also include the signal of the 26 immediately adjacent voxels ( a 3    3    3 cube centred on each voxel ) , resulting in a vector of 27 values for each voxel .
this allows us to capture information about local variations in the signal at any given point , a 
approach analogously used in multivariate functional imaging ( kriegeskorte et  al . , 2006 ) .
the result is a map of the brain where the anomaly of each voxel is labelled by its zeta score ; appropriately thresholded this yields a binary lesion map of the brain .
turning to our specific implementation of the approach , here we focus on ischaemic lesions as commonest aetiology of lesion used in neuroanatomical studies and we choose diffusion - weighted mri as the imaging modality with the highest resolution and lesion contrast within the acute period of stroke . the relative difference in lesion contrast between diffusion - weighted and other sequences in the acute period of ischaemic damage
allows us to make good use of different sequences for inter - subject registration and lesion segmentation as argued to be desirable in the foregoing .
the complex distribution of the normal diffusion signal across the brain  clearly multimodal in areas of high artefact  allows us to test the approach on data bound to be ill - suited to conventional parametric methods .
lesion - function mapping is commonly regarded as being limited primarily by the anatomical resolution of the images and the frequency at which any given area is sampled .
the former depends on the product of the physical resolution of the images and the size of the inter - subject registration error , the latter on the sample size relative to the frequency of damage to the area under investigation .
however , there is another , perhaps more important limit : the resolution of lesion sampling .
this is not the frequency at which any given area is sampled , but the volume of brain that can not be dissociated from any given area because it is always damaged together with it , at least in the sample under consideration .
this will always arise as a theoretical problem where lesion size is substantially larger than the sampling resolution : almost invariably the case in lesion studies .
the extent to which it is a problem in practice depends on correlations across voxels : essentially , the high - dimensional distribution or pattern of injury , which is unknown for vascular lesions .
our aim here is therefore to take a large sample of lesions , comparable to the largest lesion studies in the literature , and to show for the subset that involves occipital cortex how the resolvability of distinct subareas is limited by this effect .
researchers can use this map to determine whether or not a large sample of vascular occipital lesions is likely to allow them to resolve between any two subareas with the aid of lesion - mapping .
the zeta score , then , provides a simple , intuitive , and general way of segmenting lesions from images of the brain .
we proceed to give a formal description of the algorithm , and go on to evaluate its performance in the context of acute ischaemic damage imaged with diffusion - weighted mri .
finally , we apply the method to generate a probabilistic map of the damage to the brain in cases where the occipital lobe is affected , demonstrating the utility of the approach .
we have seen that the gamma score of a given point , x , is the average distance to its k nearest neighbours ( nn ) , viz : ( harmeling et  al . , 2006).(1)k(x)=1ki=1k(xnni(x ) )
the zeta score is the difference between gamma  an index of the anomaly of the point in relation to its neighbours  and the average inner - clique distance of its neighbours  an index of the density of the neighbourhood clique ( rieck and laskov , 2006):(2)k(x)=1ki=1k(xnni(x))1k(k1)i=1kj=1k(nni(x)nnj(x ) )
a key part of the elegance of this approach is that  being adaptive to the data  it has only few parameters to adjust .
so to use the zeta score in the imaging domain , we need to take only two key decisions .
the optimal value of k is determined by the smallest scale of any heterogeneity in the reference population , effectively the size of the smallest distinct cluster of values differing from the main population yet still within the normal range . since the clustering of data has no established measure no rigid guide to the selection of k
here we arbitrarily set the value of k at 15% of the size of the number of cases within the reference set ( 15 images in our case ) . the second decision to take
is how best to parameterize the signal at each voxel : since we do not need to use a scalar  we can have a vector of values  we can incorporate as few or as many parameters per voxel as is computationally practicable .
for example , if one were using diffusion tensor imaging one could use each of the eigenvalues per voxel , corresponding to the magnitude of anisotropy in each plane .
similarly , one could incorporate the values of several different types of imaging sequence , coregistered so that the voxels are all matched for anatomical location , enabling multispectral lesion segmentation . additionally or alternatively , as we do here
, one can add information about the signal of adjacent voxels , thereby incorporating local pattern anatomical information in the parameterisation of the signal at each voxel .
the distance measure we compute at each voxel here is therefore the euclidean distance in 27 dimensions  the signal at the index voxel and at the 26 voxels ( a 3    3    3 cube ) immediately adjacent to it  between each voxel in the test image and the homologous voxels and their adjacents in a set of reference , normal images coregistered with the index image so that anatomically homologous regions are brought into alignment .
the result is a voxel - wise map of zeta values from which a binary image ( lesioned vs non - lesioned ) can be created by choosing an appropriate threshold zeta value .
thus though we use only one sequence to derive zeta , the b1000 image , we parameterise each voxel as the multivariate pattern of a 27 voxel cluster . in order to evaluate the performance of zeta we need to generate a set of zeta maps that we can compare with a  ground truth  of what is really damaged and what is not .
for this we naturally need two kinds of data : a set of target images of focally lesioned brains ( in one or more imaging modality ) and a set of reference images of non - lesioned brains with characteristics otherwise similar to the target set . for the former here we used dwi data derived from 38 patients with a diagnosis of acute stroke , and for the latter dwi data from 95 patients investigated for possible stroke but with no radiological evidence of that or any other focal abnormality . the data for each patient consisted of standard clinical axially - acquired diffusion - weighted echoplanar imaging ( b0 and b1000 ) sampled at 1  mm    1  mm    6.5  mm resolution , and obtained on a ge genesis signa 1.5 t mri scanner in a single session for each patient .
the b0 images satisfactorily distinguish between cerebrospinal fluid , grey , and white matter in line with their t2-weighting while being relatively unaffected by acute ischaemic lesions in the first few hours .
the b1000 images , by contrast , show little normal tissue differentiation but marked differences in signal between damaged and undamaged tissue .
the lesioned brain images were unselected except for a minimum lesion size of 216  mm .
note this value is considerably smaller than most of the lesions investigated in clinical studies , making our test more stringent than required .
the patients with lesions were unselected except for the presence of a visible lesion on dwi , with a mean age of 62  years [ standard deviation ( sd )  =  20 ] , equal sex ratio , and were all scanned between 1 and 17  days from symptom onset ( mean  =  4.1  days , sd  =  3.8  days ) .
the 95 reference images were unselected except for the absence of a visible lesion on dwi , with a mean age of 67 ( sd  =  15 ) and sex ratio of 62 men to 33 women . the period from symptom onset to scan in the patients with lesions ( 117  days ) is broader than the conventional notion of  acute  .
for some patients , the lesion will therefore have been visible on the b0 image used in the normalisation step , where the abnormal signal may have degraded the quality of registration to some degree as previously discussed .
so as to test the real - world usefulness of the algorithm , however , where a single time point within the hyperacute stage of b0 negativity is rarely practicable , we chose a spread of time points , including those in which the b0 image will have been contaminated to some degree .
to the extent to which the quality of normalisation impacts on the performance of the algorithm better results might be achievable with a set of scans acquired only in the hyperacute stage .
the 75 control subjects were matched to the lesioned brains by being drawn from a similar clinical population : patients attending a transient ischaemic attack ( tia ) clinic who had not been found to have any acute lesions or any chronic lesions large enough to be visible on dwi .
they had a mean age of 65  years ( sd  =  16 ) and a sex ratio of 43 men to 32 women .
the target and a separate set of normal images from the control population were used to create a synthetic dataset of 2850 artificially lesioned brains where each of the 38 lesions was  transplanted  into the 75 normal brains ( brett et  al . ,
this was done by extracting the lesion signal  as defined by manual segmentation  and inserting it into the homologous region of the reference image ( see below ) .
the point of this manoeuvre was to evaluate the algorithm 's ability to distinguish lesioned from non - lesioned signal at the same location .
unless one is fortunate enough to have a pre- and post - lesion scan for several patients  which we are not  the only way of performing such an analysis is by using  chimaeric  images of the kind we describe , an approach others have recently followed ( ripolls et  al . , 2012 ) . if we had no grounds to suspect that lesion discriminability varied with location an analysis of this kind would have been unnecessary .
however , it is clear that this can not be assumed to be so , especially with dwi where some locations ( e.g. , frontal and temporal poles ) are much more prone to artefactually abnormal signal than others . the data used to apply the method so as to generate a probabilistic map of damage to the occipital lobe consisted of dwi volumes of 435 unselected patients with acute stroke admitted to university college london hospitals with a lesion detectable on dwi .
the computation of zeta and any other voxel - wise statistic that compares a test image against a reference set of images implies knowledge of the spatial correspondence between the two . in other words ,
the computation must be preceded by an inter - subject co - registration step in which spatially homologous voxels in the test and reference brain images are brought into alignment : a process conventionally referred to as normalisation ( brett et  al .
, 2001 ; ashburner and friston , 1997 , 1999 ; friston et  al . , 1995 ) .
normalisation is commonly done by finding a set of linear and non - linear transformations of the target image such that differences in the signal at each voxel between the target and reference images are minimized at the cost of some plausible level of distortion . in the presence of a lesion , normalization
can therefore be difficult since the signal within the lesioned region will inevitably tend to interfere with the signal matching process ( shen et  al . , 2007 ) .
one can minimize this effect by obtaining more than one imaging sequence and performing the normalization on the sequence that is least sensitive to the lesion , as we do here ( see below ) . alternatively or additionally
, 2001 ) or fill it in with signal from the homologous region in the contralateral hemisphere ( nachev et  al . , 2008 ) , but this of course implies knowing the spatial parameters of the lesion . another option is to combine the normalization and tissue segmentation ( white and grey matter , and cerebrospinal fluid ) into a single generative model , explicitly modelling the lesioned area as falling within a third , abnormal class .
this procedure may be used iteratively , with the lesioned class derived from one run serving as a prior for the next ( seghier et  al . , 2008 ) .
since the departure of the signal at each voxel from the expected value for grey and white matter respectively determines the extent to which it is likely to disturb the matching process , this approach would seem excellent at minimizing the impact of the lesion . as a means of lesion segmentation , however , it can only be expected to be satisfactory in the special case where the optimal normalization modality and the optimal lesion identification modality are the same .
clearly this is neither inevitable nor desirable since the more the lesion disturbs the tissue contrast the worse the normalisation , and the greater the tissue contrast the more likely it is to interfere with lesion segmentation as already argued .
our approach here is therefore to derive the normalization parameters by applying the unified normalization and segmentation procedure implemented in statistical parametric mapping 5 ( spm5  www.fil.ion.ucl.ac.uk/spm ) to the b0 image  which has low lesion contrast and good grey / white matter contrast  and then to apply the resultant normalisation parameters to the corresponding b1000 image , thereby allowing much better normalisation than could be achieved from the b1000 image alone .
note that the success of the normalisation procedure is not dependent on the b0 image 's being free of any lesion contrast because spm5s routine will automatically eliminate most such voxels during the segmentation component as already explained ( seghier et  al . , 2008 )
we have previously shown that , in keeping with their t2-weighting , clinical b0 images of this resolution can be satisfactorily normalised ( nachev et  al . , 2008 ) . clearly , this step is only practicable if the b0 and b1000 images are already in registration .
since each pair is taken in a single run , the images are bound to be close ; nonetheless we first rigid - body coregistered the b0 to the b1000 using spm5 's standard co - registration routine before performing the normalization outlined above .
note that this step is robust to the difference in lesion signal between the two scans because it is a rigid - body transformation .
all further operations were performed on the normalized images , which were resliced to 2    2    2  mm voxels using 7th order spline interpolation .
the set of artificially lesioned brain images was generated by replacing in each non - lesioned image and sequence ( b0 and b1000 ) the signal intensities in the corresponding voxels within each of the lesions , resulting in 38    75  =  2850 chimaeric images .
for the chimaeric images to be realistic the non - lesioned tissue in each donor / recipient pair must be broadly similar .
since unlike ct scans mri images do not have a standardized intensity scale ( bergeest and jager , 2008 ; nyl and udupa , 1999 ) , each transplantation step must be preceded by equating the global intensity of the donor and recipient images .
this was done by calculating the mean whole brain signal  excluding voxels falling within the lesion  of the lesioned , donor image and the mean whole brain signal of the recipient image , and multiplying the signal value at each voxel within the latter image by the ratio of the two means ( brett et  al . , 2001 ) .
this procedure ensured that the transplanted signal was as close as possible to what abnormal signal might have looked like had it been present in the recipient ; without it , unrealistic global variations in the signal would have clouded our assessment of the algorithm 's performance .
we confirmed that this procedure did not artificially introduce an artefactual difference between the lesion signal and the signal in the immediate vicinity of the lesion by comparing the ratios between intralesion ( 5 voxels inside the lesion boundary ) and perilesion ( 5 voxels outside the lesion boundary ) signal for all original lesion images and their corresponding sets of chimaerics ( paired t test p value < .01 for each set ) . in order to evaluate our segmentation algorithm
we require a  ground truth  or a standard that defines what is truly lesioned .
since even if we could have a histological reference the post - mortem specimen would be anatomically too distorted to be useable in mapping , the only standard here can be manual segmentation by eye : this is the approach every other major study in the field has taken ( brett et  al .
, 2001 ; andersen et  al . , 2010 ; bhanu prakash et  al . , 2008 ; gupta et  al . , 2008 ;
each lesion was therefore segmented by hand by a trained clinician ( ym ) with the aid of mipav 's gradient , magnitude and direction live - wire edge - detection tool ( barrett and mortensen , 1997 ; chodorowski et  al .
the manual segmentation was performed on normalized versions of each image ( see below ) .
with all the images in alignment we proceeded to calculate a voxel - wise zeta map for each lesioned image ( native and chimaeric ) .
zeta may be based on any distance measure of the signal at each voxel : here we chose the simplest : the euclidean . so as to capture not just the signal at a specific voxel but its variation in relation to its local anatomy , we parameterized the signal at each voxel as centred within a 3 by 3 by 3 cube from the b1000 image , giving 27 dimensions for each distance measure at each voxel . thus , the zeta map consists of a single scalar volume for each test subject where the value of each voxel is the zeta value for the vector of the signal at that location in the b1000 image and its 26 neighbours , the corresponding vectors in the normal set being the reference . a zeta map
one could use it directly to probe lesion  function relationships using either multivariate or mass univariate inference in the conventional way .
for the purposes of validation , however , we need to compare it to manual lesion segmentation , which necessarily generates a binary mask : we must therefore find an appropriate threshold .
clearly , if our algorithm is to be unsupervised the threshold needs to be set automatically .
the following procedure was used to find the optimal threshold for binarizing the zeta map into lesioned and non - lesioned tissue so as to allow comparison with the  gold standard  manual segmentation . note
that this is inevitably a heuristic process just as the manual process itself is heuristic , though of course the criterion is fixed .
first , we generated a whole brain map of the reference set data only , sensitive to variations in the normal signal diversity across different parts of the brain , what we call a diversity map .
this was done by iteratively calculating the zeta value of each voxel within each brain image within the reference set only , taking the remainder of the set as the reference , in essence treating each member of the set as a test volume with the remainder being the reference .
one map per image in the reference set  giving us a set of empirical distributions of zeta values at each voxel within the reference set .
these distributions were readily parameterised as generalized extreme value distributions [ equation ( 3 ) ] , with parameters location (  ) , scale (  ) , and shape (  ) fitted using maximum likelihood estimation.(3)f(x;,,)=1[1+(x)](1)1 this diversity map allowed us to index the degree to which the normal data at each anatomical location would exhibit the characteristics of an outlier , thereby allowing better discrimination of true outliers .
we thus used the diversity map to correct each test zeta image by subtracting it from the median of the diversity map , calculated independently for each voxel location .
second , this corrected zeta map was thresholded at 0 , for any negative value would imply the signal at a voxel being less anomalous than that of any within the normal reference set .
third , this coarsely binarized map was then submitted to a noise removal step where clusters with fewer than 5 voxels and/or fewer than 2 voxels with a zeta score of 70% of the maximum zeta value across the brain were removed .
fourth , the surviving clusters , taken one by one , were then subjected to a further thresholding procedure , where starting at 0 and incrementally ascending the function , the zeta value was identified at which the slope of the relation between zeta and the volume of the cluster identified as anomalous at that threshold became greater than 10 voxels / zeta ( for clusters smaller than 15  000 voxels ) or 100 voxels / zeta ( for clusters larger than 15  000 voxels ) . this cluster level zeta threshold was derived and applied independently for each cluster .
fifth , the denoising step above was repeated . finally , a  hole filling 
morphological operation was applied ( since a completely disconnected island of normal tissue can not possibly be functional ) and linear slice artefacts were automatically removed ( based on their wholly non - anatomical spatial and signal characteristics ) , resulting in a final binary lesion image .
ideally , this standard should be perfect . in the case of images of human brain lesions perfection
is not possible because there is normally no independent way of verifying the extent of the lesion but by the appearance of the lesion on the image alone .
the standard here is therefore the manual tracing of each lesion by a trained operator : this is considered the  gold standard  in the field .
what we have to determine is therefore the correspondence between two binary volume images . in line with established practice , ( anbeek et  al . , 2004 ;
, 1994 ; dice , 1945 ) here we use the following summary measures , sensitivity = truepositivestruepositives+falsenegativesspecificity = truenegativestruenegatives+falsepositivessimilarityindex=2(truepositives)2(truepositives)+falsepositives+falsenegativeswhere true positives are voxels correctly identified as lesioned , true negatives are voxels correctly identified as healthy , false positives are voxels incorrectly identified as lesioned , and false negatives are voxels incorrectly identified as healthy . for the chimaeric images ,
for the native images , since subject and lesion are the same no such separation can be made . to generate a map of lesions that involve the occipital lobe we processed 435 diffusion - weighted images of acute ischaemic lesions and selected from this set those that involved the occipital lobe as defined by the mask supplied with mricron ( brodmann areas 17 , 18 , 19 , including underlying white matter ) .
involvement was defined by a minimum of 30 voxels of the lesion falling within the volume of the mask which was satisfied by 120 lesions .
summary statistics were performed in a voxel - wise manner using custom matlab scripts , and the data was visualised with mricron ( http://www.mccauslandcenter.sc.edu/mricro/mricron ) .
the evaluation produces a set of measures for each lesion natively , and each lesion  subject combination chimerically . in the native case
this analysis shows a median sensitivity of .9602 , specificity of .9979 and similarity index ( si ) of .7342 , with corresponding standard error of the median of .0074 , .0007 and .0145 respectively ( fig .  5 ) .
in the chimaeric case , we can obtain estimates of the variability of each measure partitioned either by lesion or by subject . in the lesion
- wise analysis , there are 38 groups , corresponding to each lesion , with replications within each group corresponding to each subject . in the subject - wise analysis
there are 75 groups , corresponding to each subject , with replications within each group corresponding to each lesion .
this analysis shows an overall median for sensitivity of .9362 , specificity of .9986 and si of .7840 . as illustrated in fig . 
6a , the range of medians for the 38 groups in the lesion - wise analysis was .75781.000 , .9807.9999 and .5724.9281 for sensitivity , specificity and si respectively .
the median sensitivity , specificity and si for the above range of values were .9400 ( se  =  .0069 ) , .9988 ( se  =  .0002 ) and .7829 ( se  =  .0118 ) respectively .
this narrow variability suggests a robustness of the method to variations in lesion spatial and signal characteristics .
the subject - wise assessment showed median values of .9399 , .9988 and .7927 , with standard errors of the median now .0049 , .0001 and .0084 for sensitivity , specificity and si respectively .
6b shows the median values for the test parameters , across all lesions for each of the 75 recipient subjects .
once again , the minimal variability in the values  sensitivity .8296.9601 , specificity .9899.9992 , si .5252 .8102  suggests a resistance to noise and artefact .
finally the total computation time required to process a test brain image  spatially normalize the brain volume , create a zeta indexed map and segment into a binary volume  was in the order of 8  min .
7 shows the frequency distribution of damage in patients where the occipital lobe is affected within our dataset .
note that areas remote from the occipital lobe are frequently coincidentally damaged , and that the defining border  between occipital and other areas  is not clearly reproduced in the map .
this has implications for those wishing to use lesion - mapping so as to isolate functions peculiar to the occipital lobe . so as to clarify the nature of this effect we calculated the percentage of the total volume of each lesion involving the occipital lobe at the set criterion ( 30 voxels or more ) as a function of lesion volume ( fig . 
note that no lesion fell wholly within the occipital lobe and that smaller lesions were not more likely to be confined to the occipital lobe than larger ones .
indeed , the distribution of fractional occipital volume for lesions that involve the occipital lobe ( fig . 
a four - dimensional ( 4d ) version of the map ( incorporating the subject dimension ) , allowing researchers to determine the resolvability of areas of specific interest to them within the occipital lobe can be obtained from the authors on request .
we have proposed a simple , general , unsupervised method for segmenting brain lesions derived from a recently - described measure of anomaly and theoretically applicable where a reference set of normal data is available .
evaluation of its fidelity against manual segmentation in the context of diffusion - weighted imaging of ischaemic lesions has shown good performance across a range of subjects and lesion parameters ; indeed no better results have been reported for diffusion - weighted imaging in stroke .
the method is theoretically easy to adapt to other lesion types and imaging modalities , either alone , or in multispectral combination , and creates the possibility of high - throughput , fully - automated image processing pipelines for conducting large scale lesion  function studies .
first , we have shown that the performance metrics of zeta segmentation exceed those of other published algorithms for diffusion - weighted imaging , approaching the current gold standard 
manual segmentation  as closely as one could reasonably expect of any automated algorithm ( table 1 ) .
critically , there is little variability across subjects and lesions within which a bias could conceivably emerge , indeed less variability than one typically finds between observers ( fiez et  al . , 2000 ) .
clearly , since manual segmentation is the gold standard here we can not conclude that zeta is better than manual segmentation but it may be more consistent .
second , the naturally adaptive nature of the core anomaly metric theoretically allows zeta segmentation to handle highly heterogeneous data , adapting to any idiosyncrasy of a particular sequence , in our case with the prominent regional artefact seen on diffusion - weighted imaging .
zeta is able to distinguish artefactually high signal in such areas from true lesion signal elsewhere .
third , zeta segmentation can be applied to any dataset without prior knowledge of the spatial or signal features of the lesions .
it requires no priors for the number of lesions , their pattern of distribution , their signal intensity , or other features automated segmentation algorithms frequently demand from the user .
since lesions are both relatively rare  and so heterogeneous that their characteristics can not easily be parameterized  it is advantageous that the method does not require information that is difficult or impossible to acquire .
fourth , the high - dimensional nature of zeta makes it relatively easy to incorporate local anatomical information in the decision whether or not to label a voxel as abnormal by deriving the anomaly value not only from the signal at each voxel but also from voxels that are anatomically adjacent .
the method thus becomes sensitive to variations in the local pattern of the signal , not just its point value , just as a human operator would naturally perform the task .
fifth , zeta segmentation can easily be made multispectral : all that is required is for the different imaging modalities to be in spatial register .
since co - registering different modalities of a single subject is generally substantially easier than co - registering a single modality to a standard template this is simple to implement .
for example , in the case of stroke it is possible to compensate for signal drop out in older lesions by adding another modality , such as a t2 * sensitive sequence , in the distance matrix .
the anomaly of each point will then be calculated on the basis of both sequences .
indeed , the only limits to the number of modalities one can add are computational , and these are likely to be modest .
sixth , zeta segmentation naturally produces a real - numbered image where the anomaly of each voxel is given a continuous index .
since there are rarely any physiological grounds for drawing a sharp threshold at which neural tissue ceases to be active this seems to us a more appropriate instrument than any kind of binary measure . the zeta map can be directly used in making lesion - outcome inferences or relating patterns of damage to other factors of interest .
seventh , zeta has only one parameter to adjust , the size of the clique , k. the lower the value of k the smaller the scale of signal inhomogeneity to which the algorithm is sensitive .
we have not found a need to explore a wide range of k , but if this is necessary a univariate search would be quick and simple to implement .
it is also theoretically possible to adjust k for each voxel independently in response to sequence - specific voxel - wise variations in the pattern of inhomogeneity , perhaps because of distortion or artefact .
if binary maps are required other parameters come into play depending on the chosen heuristic for thresholding ; this is not necessary if the zeta maps are used directly , as we would suggest is the best approach . finally , zeta is computationally economical , taking approximately 8  min per image running on commodity hardware .
first , our method may be taken to assume that there is a monotonic relation between signal abnormality in a region and the probability or degree of dysfunction .
naturally , it may not be so , but until we have an independent means of establishing the link between physiology and the mri signal in a given sequence , the signal intensity is all we have .
lesion segmentation is a segmentation of images , not of the brains they imperfectly reflect .
second , the zeta anomaly score is a continuous variable with no a priori criterion on which one could discretize it . where binary maps are required a heuristic method of selecting a threshold
however , the performance of the heuristic we have adopted here is close to that of a trained operator , whose one criterion is by its very nature heuristic .
whether it would also do so for other modalities is a matter to be established empirically .
given the nature of the zeta score , it is to be expected that data with similar lesion contrast - to - noise ratios would perform similarly , but naturally this is not something that can be guaranteed .
although binarized data simplifies the business of inference  labelling every part of the brain as either completely lesioned or completely healthy 
it involves assumptions that are neither justified empirically nor plausible a priori . while the csf - filled centre of a chronic lesion clearly can not have any function whatsoever
, there is no reason to suppose that our arbitrary labelling of the margins of the lesion inevitably corresponds to a critical level of deterioration of physiological function , at least not for a great many lesions .
indeed , the habit of using discrete maps is arguably an artefact of the traditional way of segmenting lesions  by hand - drawn line . in any event
it is a feature of the general approach to making inferences using lesions in the brain , not of this particular method .
third , in common with any voxel - wise algorithm , zeta segmentation relies on the images already being spatially in register . until recently , this was a major obstacle because without knowledge of the spatial characteristics of the lesion it is difficult to minimize its impact on the normalization process .
however , spm5 's combined normalization and segmentation routine has been shown successfully to overcome this problem ( crinion et  al .
, 2007 ; andersen et  al . , 2010 ) . critically , that study used images for normalisation with very strong lesion contrast , maximally testing the algorithm 's capacity to deal with lesioned images , and so one would expect their findings to extend to other modalities .
the improvement in performance is roughly four times that of the preceding gold standard ( cf .
one can further minimize the impact of this problem by deriving the normalization parameters from a separate imaging sequence with low lesion contrast , as we have done here , and as we suggest is the theoretically optimal approach .
fourth , as with other voxel - wise methods ,  islands  of normal tissue completely disconnected from the rest of the brain by damaged tissue will nonetheless be scored as normal . while in the special case of completely surrounded islands
morphological operation we employ here , it would be difficult to construct a method that deals robustly with intermediate cases .
once again , this is not something that any automated algorithm could easily solve without detailed knowledge of the connectivity of each brain area and the location of each connecting tract : information that we do not yet have to any degree of precision . rather than making ad hoc decisions in each particular case  as is implicit in manual segmentation  it is perhaps best not to attempt an automated ( or indeed manual ) solution to this problem . given that most lesions are more or less ellipsoidal , any distortion resulting from such effects is likely to be of relatively minor importance .
fifth , zeta requires a set of normal images to use as a standard reference . while the relatively modest size employed here appears to be sufficient
clearly , the performance of the algorithm will be dependent on how closely the standard set is matched to the test images : ideally , the only difference between them should be the presence or absence of lesions of the type being sought .
a researcher considering our approach has two choices : either to use the standard reference set ( which can be obtained from the authors ) or to create a local reference set derived from normal images taken on the same scanner .
these do not need to be of normal subjects but merely of patients without demonstrable acute ischaemic lesions .
indeed , it would be desirable to choose a reference cohort that is matched to the test population in all respects except the presence of acute ischaemic lesions , following the basic principles of  control  selection in experiments generally .
we have used the new method to generate a probabilistic map of damage to the brain whenever the occipital lobe is affected .
the map shows that isolated damage to the occipital lobe is rare , and that coincident damage to neighbouring areas makes the task of dissociating their function from that of occipital cortex easier for some than for others .
we have provided a map that allows researchers to determine their likely ability to dissociate a given neighbouring area in a lesion - mapping study .
we have devised a simple unsupervised lesion segmentation algorithm based on zeta , a recently - described anomaly score .
the algorithm places minimal demands on the experimenter and the data , and has wide potential applicability .
tested against a dataset of vascular lesions captured by dwi , it compares favourably against manual segmentation across a range of lesion sizes , locations and morphologies , and exceeds the performance of any other published algorithm .
we have applied the algorithm to a set of images of vascular injury involving the occipital lobe , and derived a probabilistic map showing the resolvability of regions within occipital cortex from its neighbours . a matlab ( matlab 7.11 .
natick , ma , 2000 ) implementation of the algorithm , including a reference dataset for dwi , is available from the authors on request .
 making robust inferences about the functional neuroanatomy of the brain is critically dependent on experimental techniques that examine the consequences of focal loss of brain function . 
 unfortunately , the use of the most comprehensive such technique  lesion - function mapping  is complicated by the need for time - consuming and subjective manual delineation of the lesions , greatly limiting the practicability of the approach . here 
 we exploit a recently - described general measure of statistical anomaly , zeta , to devise a fully - automated , high - dimensional algorithm for identifying the parameters of lesions within a brain image given a reference set of normal brain images . 
 we proceed to evaluate such an algorithm in the context of diffusion - weighted imaging of the commonest type of lesion used in neuroanatomical research : ischaemic damage . 
 summary performance metrics exceed those previously published for diffusion - weighted imaging and approach the current gold standard  manual segmentation  
 sufficiently closely for fully - automated lesion - mapping studies to become a possibility . 
 we apply the new method to 435 unselected images of patients with ischaemic stroke to derive a probabilistic map of the pattern of damage in lesions involving the occipital lobe , demonstrating the variation of anatomical resolvability of occipital areas so as to guide future lesion - function studies of the region . 
Introduction
Methods
Results
Discussion
Conclusion
Role of funding source