Florence Cloppet
email: florence.cloppet@mi.parisdescartes.fr

A Boucher
email: arnaud.boucher@mi.parisdescartes.fr

Segmentation of complex nucleus configurations in biological images

Keywords: Nucleus segmentation, prior information, watershed-based segmentation, biomedical imaging

This paper presents a new segmentation method of complex nucleus configurations. The specificity of this work is to introduce prior information about the usual shape of cells nuclei, in order to optimize the selection of markers from which the flooding will start, during the watershed-based segmentation.

Introduction

Over the past fifteen years, the field of optical imaging has been developing rapidly with the advent of a wide range of new probes (especially fluorescent ones), with the advances in microscopy (immuno-fluorescence light microscope, confocal microscope, video microscopy...) and image analysis. All these advances have revolutionized our understanding and knowledge of the dynamic structures and events related to living cells.

As cell nuclear morphology is a hot topic research now in biology, the biological background of this paper is based on alterations in nuclear morphology. The knowledge of the cellular structure that underlies nuclear morphometric changes, with respect to a non pathological shape, is still limited. But such changes are frequently observed in various nuclei, including several human diseases as cancer, laminopathies, a subset of diseases caused by mutations in Lamins A/C or LMNA gene involved in the architectural proteins (Lamins A and C) of the cell nucleus membrane. Several biological functions seem to be related to major changes in the geometry of the nucleus. Nuclear size and shape are good visual descriptors of abnormal nuclear behavior, and could lead to a quick and efficient diagnosis of nucleus normality.

The development of large scale quantitative biological research or high throughput drug screens can lead to greater understanding of the mechanisms that underlie nuclear morphometric changes, and thus might help to find a treatment. However, the increasing quantity and complexity of image data render manual analysis unreasonably time consuming. Therefore, automatic techniques to analyze nuclei shape in a way that yields both quantitative and objective measurements are of considerable interest. To deal with such large screening tests, it is important to develop methods that allow highly reliable shape description and characterization [START_REF] Cloppet | Angular Bisector Network, a simplified generalized Voronoï diagram : application to processing complex intersections in biomedical images[END_REF]. This involves computer vision problems for which efficient solutions have not yet been found, such as extracting in a highly precise way the contour of each nucleus, even if the nuclei are aggregated or overlapping.

Section 2 gives some details about the typical problems encountered in the images and deals with the problem of segmentation, while section 3 describes our contribution to segment complex nucleus configurations. Section 4 illustrates the relevance of the developed method through results presented and discussed. Finally, section 5 is devoted to concluding remarks.

Segmentation

In this study, digitized pictures of immunostained nuclei (see Figure 1) were sampled using a conventional epifluorescent microscope (Leica DMR) coupled to a Princeton-Roper camera. As mentioned above, the presence of the mutated lamin protein resulting in abnormal nuclear shape, is proved by the immunodetection of lamin A/C (primary antibody directed against lamin A/C, secondary antibody coupled to the fluorescence probe FITC: Fluoresceine Iso Thio Cyanat). Other fluorescent staining, with DAPI (Di Aminido Phenyl lndol) can also be used to visualize DNA.

Fig. 1. Nuclei of cells in immunofluorescence images:

Images of lamin repartition (left) and DNA (right)

In fact, the major problem with such 2D images is that they are disturbed by out-of-focus light mostly due to the thickness of the specimen. Nuclei are sometimes occluded by overlapping material from higher or lower focal planes (see Figure 2).

Fig. 2. Occluded nuclei due to overlapping material from higher focal plane

This problem of overlapping areas could be solved by using a confocal microscope, but this makes the acquisition of images much slower. Indeed, at the moment, this is hardly compatible with high throughput drug screens. Furthermore, both sources of images (confocal or light microscope) present similar disturbances due to the auto-fluorescence of the medium, the acquisition noise induced by the camera detector, and various imperfections of the optical system.

Another source of problem is the density of the observed nuclei. In fact, in high density population, nuclei are clustered together with attaching borders (see Figure 3).Therefore, the segmentation step developed to extract each nucleus, in complex nucleus configurations, such as overlapping nuclei and/or clustered nuclei, has to be done with high precision. Segmentation remains a major challenge in image analysis. It refers to the task of detecting boundaries of objects of interest in an image. It should yield a partitioning of the image into disjoint uniform regions (according to a feature such as gray level, color or texture). It can rely both on the uniformity of the feature inside the region or on edge evidence. Segmentation has been widely studied as shown by [START_REF] Pal | A review on image segmentation techniques[END_REF][START_REF] Pal | A review on image segmentation techniques[END_REF][START_REF] Freixenet | Yet another survey on image segmentation: region and boundary information integration[END_REF]al., 2002. In this section, we will consider methods that can deal only with objects of complex topologies, and that have been proven to be powerful especially in biomedical imaging. Several approaches have been proposed, and they are mainly divided into two categories: energy-driven segmentation [START_REF] Kass | Snakes: active contours models[END_REF][START_REF] Zhu | Region competition : unifying snakes, region growing and Bayes/MDL for multiband image segmentation[END_REF][START_REF] Caselles | Geodesic active contours[END_REF] and watershed-based segmentation [START_REF] Vincent | Watershed in digital spaces, an efficient algorithm based on immersion simulation[END_REF][START_REF] Beucher | The morphological approach of segmentation: the watershed transformation[END_REF].

In the first category, the segmentation is obtained as a result of the minimization of an energy function, composed of a data-driven term and a regularization term. To prevent the problems encountered by snake-based methods with topological changes, Osher and Sethian (1988) proposed a level set approach, which represents a deformable two-dimensional contour as a level curve of a three dimensional function that is deformed by internal and external forces. The multiphase level set method [START_REF] Luminata | A multiphase level set framework for image segmentation using the Mumford and Shah model[END_REF] allows the representation of complex topologies. In the context of biological object segmentation, these methods usually failed to separate objects that are in contact or overlap. Recent methods, based on deformable contours [START_REF] Zimmer | Coupled Parametric Active Contours[END_REF] or level sets [START_REF] Zhang | Tracking fluorescent cells with coupled geometric active contours[END_REF], prevent the contours of the nuclei that are in contact from merging. They introduce in the energy functional an overlap penalty between each pair of contours. They can extract regions in contact, but cannot manage to segment overlapping regions. These methods are very useful in applications of isolated cells tracking, as it can be assumed that each cell contour is not very different between two consecutive images of a video, that is the result contours of the previous image (at time t-1) can be use as initial contours of the image at time t. The only initialization step done manually is the first one made on image at t=1. But it remains a major drawback for high throughput drug screens, which is the biological context of the study.

In the second category, the image is processed as a height function that describes a landscape [START_REF] Beucher | Use of watersheds in contour detection[END_REF]. Among the watershed-based segmentation methods, the watershed transform [START_REF] Vincent | Watershed in digital spaces, an efficient algorithm based on immersion simulation[END_REF] has proven to be a powerful and fast technique for both contour detection and region-based segmentation. An intuitive description of this transform can be done as follows: consider the image is a topographical relief, where the height of each point is directly related to its gray level, and imagine rain is falling on this relief: the watersheds are the lines that separate the catchments basins. It is generally computed on the gradient of the original image, so that the watersheds are located at high gradient points.

Many interesting properties, among which simplicity, speed and complete partition of the image, are useful for many different image segmentation applications. It has been applied to biological images [START_REF] Costa | Cell nuclei segmentation in noisy images using morphological watersheds[END_REF][START_REF] Chen | Automated segmentation, classification, and tracking of cancer cell nuclei in time-lapse microscopy[END_REF]) and has given good results, especially on blurred gray-level images with low contrast and weak boundaries. However, there are some drawbacks such as over-segmentation, or great sensitivity to noise…. In fact, fluctuations in the gradient magnitude image as well as noise can be considered as local minima, and lead to over-segmentation. The watershed segmentation can be improved, in particular by the introduction of knowledge. [START_REF] Nguyen | Watersnakes: energy driven watershed segmentation[END_REF]Ji, 2006, andHamarneh andLi, 2009 have introduced knowledge into the evolution of region boundaries during the merging process applied to limit over-segmentation A more commonly used solution is the use of knowledge to select good region markers, that represent the number and approximate location of the desired regions. [START_REF] Grau | Improved watershed transform for medical image segmentation using prior information[END_REF]al. 2004 andLefevre 2007, have shown that both markers spatial positions and content can be used as a knowledge source in order to decrease over-segmentation. The algorithm of Lefevre 2007, requires a manually made spatial initialization to determine the initial positions of the markers, and involves a supervised pixel classification from which the knowledge used to select the best markers is built. Obviously this step increases the computational cost. The method of [START_REF] Grau | Improved watershed transform for medical image segmentation using prior information[END_REF], developed for anatomical structures in medical images, used a prior model of the anatomical structure that is expected to be segmented in the form of a statistical atlas. This atlas is used to generate the markers automatically. In this paper, we propose a fully automated method of segmentation for complex nucleus configurations, especially in aggregated and/or overlapping nucleus configurations where no satisfactory solution has yet been found. This method introduces prior information based on geometrical properties, via configurations nucleus templates, and gray-level criterion that gives a semantic understanding of these configurations. This high level interpretation of the biological information will allow an automated detection of the desired markers that are connected with the number of nuclei implied in such complex nucleus configurations.

Segmenting overlapping and/or clustered nuclei

Figure 3 shows the flow chart of the full automated cell segmentation procedure. It is a three steps approach. First clusters of cells nuclei are segmented form the background, with an automated thresholding. Then the isolated cells nuclei are extracted. Finally, complex nucleus configurations, such as overlapping or clustered nuclei, are processed. The main difficulty is to first detect whether a connected component extracted on the image corresponds to a single nucleus or to a group of several nuclei. This information will be helpful for selecting the right markers in order to avoid over or sub-segmentation of cells nuclei.

Markers Selection with use of prior information

Overlapping/clustered nuclei have intersecting contours in points of high concavity. To select those involved in the segmentation process among the whole set of concavity points, prior information about normal or abnormal nuclei shape has to be used. For example, normal nuclei have a slightly elongated but regular shape with a quite a small area (see Figure 4a), while abnormal nuclei might be much bigger (see Figure 4b) and/or lobulated (See Figure 4c). In fact, points of significant concavity that are extracted can correspond to clustered or overlapping nuclei areas, but also to lobulated areas of abnormal nuclei. (See Figure 5 In order to identify these three types of areas, segments linking significant concavity points (see Figure 6a) are explored for selecting the one(s) that correspond(s) to true separation segment(s) (see Figure 6b). For this purpose, two templates of high biological plausibility (based on visual inspection of about 300 nuclei) have been established:

-overlapping template, -aggregated template.

The overlapping template (see Figure 7a), is based on the hypothesis that if two nuclei overlap, the connected component has to be elongated enough to contain two nuclei.

[AB] is the potential segment of separation that connects two significant concavity points A and B, I is the middle point of this segment, [CD] At the moment, the distances CD and AB are equal as well as EF and GH. If the nuclei were circles, we could have fixed equal values to AB, CD, EF and GH. But as they are not, EF has to be smaller than AB and has been evaluated to the third of distance AB after observing a high number of nuclei.

To select the segment [AB] as a true separation segment between overlapping nuclei, all the points C, D, E, F, G, H of the template have to be inside the connected component (see Figure 7b). The problem with this template is that in case of very big and lobulated nuclei, potential segments of separation might fit the overlapping template, though there is no overlapping area (see Figure 8). Therefore, in order to distinguish these two configurations, a study of gray levels intensity along [AB] is performed. In case of overlapping nuclei, the gray level intensity is higher in the area of overlapping than the mean intensity of the connected component. In case of a big lobulated nucleus, no significant variation of gray levels can be observed, and then the potential segment of separation is rejected (see Figure 8c).

The aggregating template (see Figure 9a) illustrates the fact that in aggregating areas, there is a narrow pathway where the two nuclei are in contact. This template is more restrictive than the overlapping template. If the aggregating template fits (see Figure 9b), the segment [AB] is validated as a line of separation between nuclei in contact (see Figure 9c). As the template could be too rigid, some flexibility is introduced, and a rule is applied as follows: "Only one point among the points C, D, E, F, G, H can be outside the connected component".

Watershed transform

If the connected component includes several nuclei, they have to be separated. The watershed transform is then used. The critical step of this transform is the setting of the markers which are the starting points of flooding.

The marker image used for watershed segmentation is a binary image consisting of either single marker points (generally the ultimate point) or larger marker regions, where each connected marker is placed inside an object of interest. Each initial marker has a one-to-one relationship to a specific watershed region, thus the number of markers will be equal to the final number of watershed regions.

In the developed method, marker regions, based on knowledge extracted at the previous step about nuclei configurations, are used rather than ultimate points, which failed in case of complex nucleus configurations to extract the good number of nuclei. These markers regions are usually extracted as low gradient zones in the image, which is the case in the middle of the nucleus. Indeed, the selected marker regions are the segments [OM] and [NP], which are by definition positioned approximately at the middle of each region located on both sides of the validated segment of separation (see Figure 10), and which width is fixed to 10 pixels to prevent them to be trapped in a local minimum.

I: middle point of [AB] (segment of separation). J: middle point of [AJ]. K: middle point of [IB]. d: perpendicular line to [AB] at point I. E, F: points of intersection between the line parallel to d at point J and the contour of the connected component or another validated segment of separation. G, H: points of intersection between the line parallel to d at point K and the contour of the connected component or another validated segment of separation. M, N: middle point of [EJ] and [JF] respectively. O, P: middle point of [GK] and [KH] respectively.

Fig. 10. Definition of marker regions

If the nuclei have been identified at the previous step as aggregating nuclei, a marker is set in the middle (see Figure 10) of the regions located on both sides of the validated segment of separation. If the nuclei overlap, there will be one more marker in the overlapping area along the segment [AB] (see Figure 11). Then the gradient image is considered as a topographical surface that is flooded from the markers. Finally, a post-processing is applied to connect each region to the corresponding nucleus (see Figure 12): regions 1 and 3 form the first nucleus and regions 2 and 3 form the second one.

Results

To test the efficiency of the proposed method, 50 digitized pictures (1392x1040 pixels) of immunostained nuclei, sampled using a conventional epifluorescent microscope (Leica DMR) coupled to a Princeton-Roper camera, were tested and evaluated. It is not possible to obtain a ground truth based on experts manual segmentation on a great number of nuclei Therefore, on these images containing 712 nuclei, a visual analysis of the results was done by three biological experts, in order to assess the relevance of the extracted nuclei in terms of number of nuclei, but they also have to check that no minor deformation of nuclei shape is missed. As these minor deformations have a fundamental influence in drawing the right conclusions about shape normality, it is important to well extract them during the segmentation step. This expert analysis gives in results the following nuclei repartition (See Table 1 The results of the proposed system (see Table 2) show that 97.47% of nuclei are well segmented in terms of number of nuclei, according to biological experts. Furthermore, no minor deformation is missed. All isolated nuclei are well extracted, so that the 2.53% of nuclei that are not well segmented are implied in complex nucleus configurations. The rate of over-segmentation (i.e. individual nuclei that are split in more than one particle) is very low (0.28%), and show the relevance of the extracted markers. The remaining cases of over-segmentation are observed in some very lobulated nuclei, where two high concavity points match the aggregating template (see Figure 13a) and satisfy the gray-level criterion.

Sub-segmentation (i.e. nucleus configurations are not split in the right number of individual nuclei) rate is slightly higher. According to Table 3, it can be shown that most of the errors of sub-segmentation are observed during the segmentation of overlapping nuclei. They are due to the fact that only one point of high concavity is detected (see Figure 13 b), therefore no segment of separation can be extracted. The average computational cost of the proposed method has been measured around 3 sec per aggregated or overlapping nucleus configurations, with a Java-based implementation on a Pentium T8100 2.1 GHz / 1 GB RAM laptop. For a whole image containing approximately 30% of nuclei involved in complex nucleus configurations, the average computational cost has been measured around 30 sec per image.

Conclusion

In this paper we tackle the problem of segmentation of complex nuclei configuration, such as aggregated or overlapping nuclei, where no satisfactory solution was found in the past. The developed method is based on the marker based watershed algorithm which is one of the most widely accepted solutions in biomedical imaging. This kind of algorithm requires a spatial initialization made either automatically or manually to define the initial position of the catchment basins (i.e. the markers).

We propose a method based on an automatic generation of markers introducing prior information based on geometrical properties by the way of configuration nuclei templates, and on gray-level criterion. Prior information gives a semantic understanding of the nucleus configurations that is useful to determine the right number of markers and their approximate location. This method is not time consuming and yield to a very low rate of over-segmentation.

Future works will be devoted to the introduction of a formal learning step in order to obtain optimal aggregated/overlapping templates based on biological knowledge extracted as rules or image samples. That will make the general application of the proposed method even more reliable and flexible as it will be able to adapt to new nucleus configurations.

Regularization of watershed lines [START_REF] Nguyen | Watersnakes: energy driven watershed segmentation[END_REF][START_REF] Tai | Level set methods for watershed image segmentation[END_REF] can also be introduced in order to smooth the extracted contours, so as to limit inaccuracies in subsequent measurements of shape features.

Fig. 3 .

 3 Fig. 3. Block diagram of the segmentation system

Fig. 5 .

 5 ) a-slightly elongated but regular bbig c-big and lobulated Fig. 4. Domain knowledge about normal (a) and abnormal (b-c) nuclei shape a b c Points of significant concavity in overlapping nuclei (a) clustered nuclei (b) and lobulated abnormal nucleus (c)

Fig. 6 .

 6 Fig. 6. Selection of a true separation segment [AB] (b) among all potential segments linking two significant concavity points (a) in various nuclei configurations

  is a segment perpendicular to [AB] at point I, [EF] and [GH] are perpendicular to [AB] respectively at points A and B.

Fig. 7 .

 7 Fig. 7. Use of the overlapping template (a) to select (b) the true separation segment (c) in overlapping nuclei configurations

Fig. 8 .

 8 Fig. 8. In case of very big and lobulated nucleus potential segments of separation (a) can fit the overlapping template (b) but don't satisfy the gray-level criterion (c)

  It differs by the angles the consistent segments [EA][EF] and [GB][GH] make with the potential segment of separation [AB]. The distances CD, EA, AF, GB, BH (with equal values for EA, AF, GB, BH) are also evaluated a high number of nuclei. The angles are measured using left and right tangent vectors at concavity points A and B.

Fig. 9 .

 9 Fig. 9. Use of the Aggregating template (a) to select (b) the true separation segment (c) in aggregated nuclei configurations

Fig. 11 .

 11 Fig. 11. Initializing step of the watershed transform in case of overlapping nuclei

Fig. 12 .

 12 Fig. 12. Result of overlapping nuclei segmentation

  Example of over-segmentation (a) and, sub-segmentation (b)

Table 1 : Nucleus distribution according to biological experts

 1 ).

	Isolated Nuclei	Aggregated Nuclei	Overlapping Nuclei	Total
	553 (77.67%)	89 (12.50%)	70 (9.83%)	712

Acknowledgements

The authors would like to thank N. Levy, P. Cau, C. Navarro and S. Perreira (Inserm U910 "Génétique Médicale et Développement", Faculté de Médecine de la Timone, 13385 Marseille Cedex 05) for providing images, and for their valuable suggestions and inspiring conversations.


A Boucher

F Cloppet

N Vincent
email: nicole.vincent@mi.parisdescartes.fr

P Jouve
email: pjouve@fenics.fr

Visual perception driven registration of mammograms

Keywords: registration, control points, mammogram, anatomical features I

This paper aims to develop a methodology to register pairs of temporal mammograms. Control points based on anatomical features are detected in an automated way. Thereby, image semantic is used to extract landmarks based on these control points. A referential is generated from these control points based on this referential the studied images are realigned using different levels of observation leading to both rigid and non-rigid transforms according to expert mammogram reading.

INTRODUCTION

Breast cancer is the most prevalent cancer among women [START_REF] Suri | Recent Advances in Breast Imaging, Mammography, and CAD of Breast Cancer[END_REF]. The majority of new breast cancers are diagnosed as a result of an abnormality seen on a mammogram, it is at the moment the best way to detect breast cancers in early stages. This has yielded to a resultant improvement in survival rates. Mammograms are X-ray images, where 3D-objects such as pectoral muscle, nipple, and adipose tissue, are projected in 2D dimensions (see Fig. 1). Several mammograms are done during a woman life. Such temporal information is used by radiologists. To provide an efficient solution to the problem of comparison of digital mammograms, a registration step is necessary. It is useful to reduce differences due solely to positioning the breast with respect to the device while preserving differences in mammograms caused by breast changes due to the apparition or progression of a disease. This is an important issue before a computed-aided diagnosis can be performed.

The purpose of this paper is to develop an automated method to register a couple of temporal mammograms in mediolateral-oblique (MLO) view. The method is inspired by experts reading methods and based on the patient anatomy. After presenting the problem of mammogram registration in section 2, our contribution will be detailed in section 3. The results are discussed in section 4.

II. MAMMOGRAM REGISTRATION

The mapping of information given by several multi or mono modal images has been an active research topic in the biomedical field. If no registration is applied, it is difficult to localize and compare regions of interest on two images, as differences of size, orientation or other distortions may have occurred. Brown [START_REF] Brown | A survey of image registration techniques[END_REF] defines the transform associated with a registration as follows:

t 0 = argmax t ∈ T S ( I , t(J) ) (1)
where I and t(J) are the referent and transformed images, compared using a similarity measure S.

It consists in an estimation of an optimal transform t belonging to T, which allows the mapping of two images I and J; and argmax represents a resolution strategy. As illustrated by Figure 2, a transform t can be rigid (translation, rotation), or non rigid (affine : global linear transformation ; deformation field). On the one hand, rigid and affine transforms are widely used to recover global deformations, e.g. for intra-patient registration. However, they do not cope with local deformations. On the other hand, non-rigid transforms address the problem of local deformations; however they match contours of structure with local irregularities due to the high number of degrees of freedom. Various methods have been developed to correlate areas in mammograms in order to detect differences that are likely to be associated with an abnormality. A survey of these methods can be found in [START_REF] Suri | Recent Advances in Breast Imaging, Mammography, and CAD of Breast Cancer[END_REF].

Most of the rigid methods rely on extraction of corresponding landmarks (points, curves or surfaces). These landmarks quality plays an important role in the alignment accuracy. Tedious is manual landmarks generation and it is a time-consuming task, therefore a large number of methods are based on an automated extraction of landmarks in mammograms. According to [START_REF] Wirth | A non-rigid approach to medical image registration : Matching image of the breast[END_REF], the only robust corresponding landmarks for mammogram registration are anatomical features such as breast boundary, pectoral muscle, nipple, crossing points of horizontal and vertical structures, fibrous tissue, ducts and vessels, the latest are less robust in case of controlateral mammograms registration, due to the natural dissymmetry between left and right breasts. A lot of methods have been developed to extract robust anatomical features [START_REF] Kwok | Automatic pectoral muscle segmentation on mediolateral oblique view mammograms[END_REF][START_REF] Mirzaalian | Pectoral muscle segmentation on digital mammograms by nonlinear diffusion filtering[END_REF][START_REF] Mustra | Breast border extraction and pectoral muscle detection using wavelet decomposition[END_REF]. These robust anatomical features are not homogeneously spatially distributed. Therefore a rigid registration of mammograms needs an intelligent use of these anatomical features.

Most of non-rigid registration methods are iterative and minimize a cost or an energy function defined in terms of the intensity difference between images. The similarity criterion may involve mutual information (2) that implies entropy computation (3):

IM(I, J) = H(I) + H(J) -H(I,J) (2) HሺI ୩ ሻ ൌ െ ∑ p ୧ େ . logሺ p ୧ ሻ (3)
where I and J are the processed images, C is the set of n gray levels, and i is comprised between 1 and n.

Other similarity measures can be used such as the sum of absolute differences, the sum of squared differences, a correlation coefficient or the Wood's criteria [START_REF] Woods | Rapid Automated Algorithm for Aligning and Reslicing PET Images[END_REF]. As these methods imply an iterative minimization of an energy function, they are quite time consuming. To decrease the iteration number, some methods use a rigid approach to initialise the search in a non-rigid method [START_REF] Chiou | Hybrid Registration of Corresponding Mammogram Images for Automatic Detection of Breast Cancer[END_REF]. Moreover, some biomechanical models [START_REF] Rajagopal | The breast biomechanics reference state for multi-modal image analysis[END_REF][START_REF] Lee | Correlation of breast image alignment using bio-mechanical modeling[END_REF] have been proposed to improve non-rigid registration, with adaptation to breast compression that implies deformation of the soft tissue. Another drawback of non-rigid methods is the potential modification of abnormalities evolving over time.

In the image mammograms registration, very important is to preserve differences caused by breast changes or new disease appearance or progress of existing disease. To preserve as much as possible these mammograms differences not caused by positioning, we propose a fully automated coarse-to-fine registration method. It is based on robust anatomical features linked to intensity and on semantical interpretation. We have been inspired by expert's perception and reasoning way.

Then local non-rigid transforms are possible. The joint use of a global rigid registration followed by a local affine method will allow us to simulate a deformation field while maintaining the processing speed of rigid methods.

III. VISUAL PERCEPTION DRIVEN METHOD

Most of the registration methods look for matching images either at a global level or at a local level according to a function to be optimized. Of course this is not consistent with expert reasoning and human vision processing. Different observation levels have to be taken into account. A global view focuses on anatomical elements that can easily be matched and that give sense to the different parts of the breast. Indeed, texture must not be used at this level as it is too much involved in the defaults that are the core of the mammogram reading. In a second step, some semantical meaning is associated with each region, with possibly more local processes. A raw matching is then done based on a rigid anatomic referential (a primary model) associated with each breast.

To simulate the estimation of a deformation field some independent affine transforms are applied to regions linked to this referential, as the experts practice. This approach allows simulation of elastic deformation and maintains a low computation time. First, the anatomical elements enable the building of a referential associated with each breast image. The breast boundary, the pectoral muscle, the nipple, the lower (fold between the belly and the breast) and higher limits of the breast (see Fig. 4) can be extracted in a robust way. They are referential points commonly used by experts.

As the paper aim is not anatomical features extraction, we only here describe briefly this step. This line is used as the Y axis of the referential to be built. The line perpendicular to Y axis and going through the nipple centre of gravity is chosen as X axis (see Fig. 5). The lower and higher limits of the breast are localized under and above the nipple axis where the breast thickness (distance from boundary breast to Y axis) is the smallest. The referential highlights two parts in the breast, the upper Z u and the lower part Z l .Whereas they have the same texture, they are not deforming in the same way according to the acquisition process.

The registration is based on a patient referential and a two-step matching process. The first matching, scale independent, is a rigid transform. A translation, matches the origins of the referentials and a central rotation superimposes the axis. Let F be this transform. In the second step, only local processes are considered. The regions are inspired by the way the experts are reading the images. The experts vertically scan the image, therefore the Y axis is promoted. Then the scanning is modeled by way of horizontal rectangular zones (Z i ), more precisely the zones are one pixel high. The images I and J are defined as:

I = ڂ Z ୪୍ ୪∈୍ J = ڂ Z ୪ ୪∈ (4)
where l is the index of the zone Z along the Y axis.

As noticed above, the stretching along the Y axis is not the same in upper part Z u and in lower part Z l . Note these two pressures V Zu and V Zl . V a (Z lI ) will be registered on Z lJ with a taking Z u or Z l value.

The transformation between two segments with equal y coordinates is modeled by a linear scaling, in the X axis direction. The scale changes are not identical all along the breast, figuring the different pressure the breast receives. We note H b the scale change applied on the horizontal with b ordinate.

Let S be a horizontal segment in image I :

S' = V a [ H b (S) ] (5)
where a designates Z u or Z l and b is the y coordinate of S in the referential of image J. The breast referential allows to reduce differences due solely to positioning (orientation) of the breast during the acquisition step (see Fig. 6), this is the coarse registration step. Then, during the fine registration step, the deformations are locally defined. Evaluating registration algorithms is difficult due to the lack of gold standard in most clinical procedures [START_REF] Betrouni | Multimodality medical images registration : From conception to validation[END_REF]. The bronze standard is a real-data based statistical method that provides an alternative registration reference through a computationally intensive image database registration procedure. It is a statistical procedure whose purpose is to estimate the registration algorithms performances by considering the exact result as a hidden variable estimated from a large enough data set. The developed method was tested on 30 couples of temporal mammograms in MLO view. Experts have selected manually five landmarks (micro calcifications, macro calcifications …) present on both images. Distance before registration is computed by superimposing, in an absolute and global way, the two images of each couple, independently of their contents. In order to evaluate our results (see Table 1 and Fig. 7), we compare them with those obtained with a global translation matching the nipples in the two images, as practiced, in a mental way, by radiologists.

The results show a diminution of the mean distance between the selected landmarks. Our method gives better results than a rigid registration either for mean distance between landmarks or standard deviation. Maximum error is also reduced (Fig. 7). The processing time of our method is higher than with translation registration but is still reasonable, it is far lower than the processing time of complex non-rigid methods, which is about 4 minutes in the most favorable case with parallel implementation [START_REF] Stefanescu | A Grid Service for the Interactive use of a Parallel Non-Rigid Registration Algorithm of Medical Images[END_REF]. We propose a fully automated method to register mammograms. This method is a coarse-to-fine approach. It is a globally rigid method based on robust anatomical features that minimizes differences caused by positioning and compression as much as possible. Some deformation based on linear scalings can then be applied locally, as some semantical meaning has been associated with each region at the previous step. This method preserves in mammograms differences really caused by breast changes and new disease appearance or progress of existing disease.

Distance

Future works will be devoted to define more precisely the non-rigid transforms. Indeed, the different lines are not independent and we will involve some texture elements to have an even more realistic vision.

Figure 1 -

 1 Figure 1 -Example of a MLO view mammogram

Figure 2 -

 2 Figure 2 -Types of transformations

Figure 3 -

 3 Figure 3 -Block diagram process

Figure 4 -

 4 Figure 4 -Anatomical features extracted Thresholding enables the breast boundary extraction. The nipple detection is based on a patch learning process.

Figure 5 -

 5 Figure 5 -Orthogonal system used as referential

Figure 6 -

 6 Figure 6 -Registration process of image I in image J IV. RESULTS

Figure 7 -

 7 Figure 7 -Distance between couple of selected landmarks grouped in quartiles

Table 1 -

 1

		Before	Translation	Our
		registration	registration	method
	Mean	1.02	0.68	0.62
	Standard Deviation	0.58	0.40	0.33
	Mean execution time (sec)	-	3.05	3.95

Comparative results: mean distance (cm) and execution time (sec) for a couple of temporal mammograms

ACKNOWLEDGMENTS

We would you like to thank the ANRT, the French National Agency of Research and Technology for its financial support. We also want to thank Dr A. Bremont for his valuable suggestions and inspiring conversations.


Michal Valko
email: michal@cs.pitt.edu

Gregory F Cooper

Amy Seybert
email: seyberta@upmc.edu

Shyam Visweswaran

Melissa Saul

Milos Hauskrecht

Conditional anomaly detection methods for patient-management alert systems

Keywords: anomaly detection, alert systems, monitoring, health-care applications, metric learning

come

Introduction

Anomaly detection methods can be very useful in identifying interesting or concerning events. Typical anomaly detection attempts to identify unusual data instances that deviate from the majority of examples in the dataset. Such instances indicate anomalous (out of ordinary) circumstances, for example, a network attack [START_REF] Eskin | Anomaly detection over noisy data using learned probability distributions[END_REF] or a disease outbreak [START_REF] Wong | Bayesian network anomaly pattern detection for disease outbreaks[END_REF]. In this work, we study conditional anomaly detection [START_REF] Hauskrecht | Evidence-based anomaly detection[END_REF] framework that extends standard anomaly detection by identifying partial patterns in data instances that are anomalous with respect to the remaining data features. Such a framework is particularly promising for identifying unusual patient-management decisions or patient outcomes in clinical environment [START_REF] Hauskrecht | Evidence-based anomaly detection[END_REF].

Our conditional anomaly detection approach is inspired by classification model learning. Let x defines a vector of input attributes (representing the patients state) and y defines the output attribute (representing the target patient-management decision). Our goal is to decide if the example (x, y) is conditionally anomalous with respect to past examples (patients) in the database. In other words, we ask if the patient management decision y is unusual for the patient condition x, by taking into account records for past patients in the database. Our anomaly detection framework works by first building a discriminative measure d(•) that reflects the severity with which an example differs from conditional (input-to-output) patterns observed in the database. All anomaly calls are then defined relative to this measure. To construct d we rely on methods derived from classification model learning.

In particular, our method exploits discriminant functions often used to make classification model calls. We investigate and experiment with discriminative measures derived from two classification models: the Naïve Bayes model [START_REF] Domingos | On the optimality of the simple bayesian classifier under zeroone loss[END_REF] and the support vector machines [START_REF] Vapnik | The nature of statistical learning theory[END_REF].

The anomaly detection call for the current instance (patient) can be made with respect to either all patients in the database or their smaller subset. In this work we pursue instance-based anomaly detection approach. The instance-based methods do not try to learn a universal predictive model for all possible patient instances at the same time, instead the model is optimized for every data instance (patient) individually. The instance-specific model M x may provide a better option if the predictive model is less complex and the dataset is small [START_REF] Aha | Instance-based learning algorithms[END_REF].

An instance-specific methods typically rely on a distance metric to pick the examples most relevant for the comparison. However, standard distance metrics such as Euclidean or Mahalanobis metrics are not the best for the anomaly detection task since they may be biased by feature duplicates or features that are irrelevant for predicting the outcome y. Thus, instead of choosing one of the standard distance metrics we investigate and test metric-learning methods that let us adapt predictive models to specifics of the currently evaluated example x.

We investigate two metric-learning methods that were originally used for building non-parametric classification models. The first method is NCA [START_REF] Goldberger | Neighbourhood components analysis[END_REF]. The method adjusts the parameters of the generalized distance metric so that the accuracy of the associated nearest neighbor classifier is optimized. The second method, RCA [START_REF] Bar-Hillel | Learning a mahalanobis metric from equivalence constraints[END_REF] optimizes mutual information between the distribution in the original and the transformed space with restriction that distances between same class cases do not exceed a fixed threshold. We test the methods and show their benefits on two real-world problems: identification of unusual patient management decisions for (1) patients suffering from the community acquired pneumonia, and (2) post-surgical cardiac patients on the Heparin therapy.

Methodology

Conditional anomaly detection

The objective of standard anomaly detection is to identify a data example a that deviates from all other examples E in the database. Conditional anomaly detection [START_REF] Hauskrecht | Evidence-based anomaly detection[END_REF] is different. The goal is to detect an unusual pattern relating input attributes x and output attributes y in the example a, that deviates from patterns observed in other examples in the database. To assess the conditional anomaly of a we propose to first build (learn) a one-dimensional projection d(•) of the data that reflects the prevailing (or expected) conditional pattern in the database for y given x. The projection model d is then used to analyze the deviations of a's to determine the anomaly. We say that the case a is anomalous in the output attribute(s) y with respect to input x, if the value d(y|x) falls below certain threshold. Our conditional anomaly detection framework can be used for a number of purposes. Our objective here is to use it detect anomalous patient-management decisions. In this case the input attributes x define the patients condition and the output attribute y corresponds to the patient-management decision we want to evaluate.

Discriminative projections

In our work we consider two methods for building discriminative projections d(•). Both of these methods are derived from the models used frequently in classification model learning: the Naïve Bayes model [START_REF] Domingos | On the optimality of the simple bayesian classifier under zeroone loss[END_REF] and the support vector machines (SVM) [START_REF] Vapnik | The nature of statistical learning theory[END_REF]. The fact that we use classification models is not a coincidence. Classification models attempt to learn conditional patterns in between inputs x and class outputs y from the past data and apply them to predict the class membership for the future inputs. In our case, we aim to model the relation between input x and output patterns y and apply it to detect pattern deviations in the new example (x, y). In both cases the model learning attempts to capture the prevailing conditional patterns observed in the dataset and the difference is in how the learned patterns are used in the two frameworks.

Naïve Bayes model

A Naïve Bayes classifier [START_REF] Heckerman | A tutorial on learning with bayesian networks[END_REF] is a generative classification model used frequently in machine learning literature and comes with excellent discriminative performance on many ML datasets. The Naïve Bayes model is a special Bayesian belief network [START_REF] Pearl | Probabilistic reasoning in intelligent systems: networks of plausible inference[END_REF][START_REF] Lauritzen | Local computations with probabilities on graphical structures and their application to expert systems[END_REF]) that defines the full joint probability of variables x and the class variable y as: P (x, y) = P (y)P (x|y) = P (y)

k i=1 P (x i |y)
The model is fully defined by the following set of pa-rameters: (1) prior distribution on class variable and (2) class-conditional densities for all features x. This decomposition reflects the major assumption behind the model: all features (attributes) of x are independent given the class variable y we would like to predict. We note that any probabilistic calculation can be performed once the full joint model is known. The parameters of the Naïve Bayes model can be learned using the maximum likelihood or the Bayesian approaches from the training data. We adopt the Bayesian framework to learn the parameters of the model and compute any related statistics. Let M define the Naïve Bayes model. In such a case the parameters θ M of the model M are treated as random variables and are described in terms of a density function P (θ M |M ). To simplify the calculations we assume [START_REF] Heckerman | A tutorial on learning with bayesian networks[END_REF] (1) parameter independence and (2) conjugate priors. In such a case, the posterior follows the same distribution as the prior and updating reduces to updates of sufficient statistics. Similarly, many probabilistic calculations can be performed in the closed form. The Naïve Bayes model predicts the class y by calculating the class posterior P (y|x). If one model is used then the class posterior is calculated as:

P (y|x) = P (y)P (x|y) P (x) ∝ P (y) k i=1 P (x i |y)
The Naïve Bayes model can be adopted for the anomaly detection purposes by defining the discriminative projection of an example (x, y) to be equal to the class posterior, that is: d(y|x) = P (y|x). In this case the projection has an intuitive probabilistic interpretation: an example (x, y) is anomalous if the probability of the decision y with respect to its input attributes x and past examples in the database is small. Moreover, the smaller is the probability, the more likely is the anomaly. We note that the Naïve Bayes model described here can easily extend to more complex generative models based on the Bayesian belief networks.

Support vector machines

The support vector machine (SVM) [START_REF] Vapnik | The nature of statistical learning theory[END_REF][START_REF] Burges | A tutorial on support vector Conditional anomaly detection methods for patient-management alert systems machines for pattern recognition[END_REF] is a discriminative machine learning model very popular in the machine learning community primarily thanks to its ability to learn highquality discriminative patterns in high-dimensional datasets. In our work we adopt the linear support vector machine algorithm to build the conditional projection d for the anomaly detection purposes.

The linear support vector machine learns a linear decision boundary that separates the n-dimensional feature space into 2 partitions corresponding to two classes of examples. The boundary is a hyperplane given by the equation

w T x + w 0 = 0,
where w is the normal to the hyperplane, and w 0 is the distance separating the "support vectors"a set of representative training examples from each class which are most helpful for defining the decision boundary. The parameters of the model (w and w 0 ) can be learned from the data through quadratic optimization using a set of Lagrange parameters [START_REF] Vapnik | The nature of statistical learning theory[END_REF]. These parameters allow us to redefine the decision boundary as

w T x + w 0 = i∈SV αi y i (x T i x) + w 0 ,
where only samples in the support vector set (SV ) contribute to the computation of the decision boundary.

To support classification tasks, the projection defining the decision boundary is used to determine the class of a new example. That is, if the value

w T x + w 0 ≥ 0 is positive then C(x) belong to one class, if it is neg-
ative it belongs to the other class. However, in our conditional anomaly framework we use the projection itself for the positive class and the negated projection for the negative class to measure the deviation:

d(y|x) = y(w T x + w 0 ), where y ∈ {-1, 1}
In other words, the smaller the projection is the more likely is the example anomalous. We note that the negative projections correspond to misclassified examples.

Instance-specific models

Discriminative models used for anomaly detection purposes can be of different complexity. However, if the dataset used to learn the model is relatively small, a more complex model may become very hard to learn reliably. In such a case a simpler parametric model with a smaller number of parameters may be preferred. Unfortunately, a simpler model may sacrifice some flexibility and its predictions may become biased towards the population of examples that occurs with a higher prior probability. To make more accurate predictions for any instance, we resort to instance-specific predictive methods and models [START_REF] Aha | Instance-based learning algorithms[END_REF]. The models in instance-based methods are individually optimized for every data instance x. To reflect this, we denote the predictive model for x as M x . The benefit of the instance-based models is its more accurate fit to any data instance; the limitation is that the models must be trained only on the data that are relevant for x.

Choosing the examples that are most relevant for training the instance-specific model is the bottleneck of the method. We discuss methods to achieve this later on.

Selecting relevant examples

3.1. Exact match.

Clearly, the best examples are the ones that exactly match the input attributes of the instance x. However, it is very likely that in real-world databases none or only few cases match the target case exactly so there is no or a very weak population support to draw any statistically sound anomaly conclusion.

Similarity-based match

One way to address the problem of insufficient population available through the exact match is to define a distance metric on the space of attributes C(x) that lets us select the examples closest to the target example x. The distance metric defines the proximity of any two cases in the dataset, and the k closest matches to the target case define the best population of size k. Different distance metrics are possible. An example is the generalized distance metric r 2 defined:

r 2 (x i , x j ) = (x i -x j ) T Γ -1 (x i -x j ), (1)
where Γ -1 is a matrix that weights attributes of patient cases proportionally to their importance. Different weights lead to a different distance metric. For example, if Γ is the identity matrix I, the equation defines the Euclidean distance of x i relative to x j . The Mahalanobis distance [START_REF] Mahalanobis | On the generalized distance in statistics[END_REF] is obtained from (1) by choosing Γ to be the population covariance matrix Σ which lets us incorporate the dependencies among the attributes.

The Euclidean and Mahalanobis metrics are standard off-shelf distance metrics often applied in many learning tasks. However, they come with many deficiencies. The Euclidean metric ignores feature correlates which leads to "double-counting" when defining the distance in between the points. The Mahalanobis distance resolves this problem by reweighting the attributes according to their covariances. Nevertheless, the major deficiency of both Mahalanobis and Euclidean metrics is that they may not properly determine the relevance of an attribute for predicting the outcome attribute y.

The relevance of input attributes for anomaly detec-tion is determined by their influence on the output attribute y. Intuitively, an input attribute is relevant for the output y if is able to predict or help to predict its changes. To incorporate the relevance aspect of the problem into the metric we adapt (learn) the parameters of the generalized distance metric with the help of examples in the database.

Metric-learning

The problem of distance metric learning in context of classification tasks has been studied by [START_REF] Goldberger | Neighbourhood components analysis[END_REF] and [START_REF] Bar-Hillel | Learning a mahalanobis metric from equivalence constraints[END_REF]. We adapt these metric learning methods to support probabilistic anomaly detection. In the following we briefly summarize the two methods. [START_REF] Goldberger | Neighbourhood components analysis[END_REF] explores the learning of the metric in context of the nearest neighbor classification. They learn a generalized metric:

d 2 (x 1 , x 2 ) = (x 1 -x 2 ) T Q(x 1 -x 2 ) = (x 1 -x 2 ) T A T A(x 1 -x 2 ) = (Ax 1 -Ax 2 ) T (Ax 1 -Ax 2 )
by directly learning its corresponding linear transformation A. They introduce a new optimization criterion (NCA), that is, as argued by the authors, more suitable for the nearest-neighbor classification purposes. The criterion is based on a new, probabilistic version of the cost function for the leave-one-out classification error in the k-NN framework. Each point i can now select any other point j with some probability p ij defined as softmax function over distances in the transformed space:

p ij = exp(-||Ax i -Ax j || 2 ) k =i exp(-||Ax k -Ax j || 2 )
A linear transformation A is then sought to maximize the expected number of correctly classified cases (with k-NN):

arg max A g(A) = arg max A i j∈Ci p ij
where C i is the set of cases that belong to the same class as i. Intuitively, the criterion aims to learn a generalized distance metric by shrinking the distance between similar points to zero, and expanding the distance between dissimilar points to infinity.

The algorithm and the metric it generates was shown to outperform other metrics for a number of learning problems. The method climbs the gradient of g(A), which is (x ij being x i -x j ): [START_REF] Bar-Hillel | Learning a mahalanobis metric from equivalence constraints[END_REF] and [START_REF] Shental | Adjustment learning and relevant component analysis[END_REF] define a different optimization criterion based on the mutual information. The advantage of their method (relevant component analysis -RCA) is the existence of the closed form (efficient) solution. Under the mutual information criterion, the class information is incorporated and optimized by computing the averages of class covariance matrices. The resulting matrix is obtained by

∂g ∂A = 2A i   p i k p ik x ik x T ik - j∈Ci p ij x ij x T ij   (Bar
Σ RCA = k i=1 Σi A = Σ -1 2
where Σi is the sample covariance matrix of class i and A is the resulting transformation for the data. The disadvantage of the method is that it assumes Gaussian distribution for the classes.

Experimental evaluation

We test anomaly detection framework and its the instance-based methods on the problem of identification of anomalous patient-management decisions for two real-world clinical datasets.

Pneumonia PORT dataset

The Pneumonia PORT dataset is based on the study conducted from October 1991 to March 1994 on 2287 patients with community-acquired pneumonia from three geographical locations at five medical institutions. [START_REF] Kapoor | Assessment of the variantion and outcomes of pneumonia: Pneumonia patient outcomes research team (PORT) final report[END_REF][START_REF] Fine | A prediction rule to identify low-risk patients with community-acquired pneumonia[END_REF]. The original PORT data were analyzed by [START_REF] Fine | A prediction rule to identify low-risk patients with community-acquired pneumonia[END_REF], who derived a prediction rule with 30-day hospital mortality rate as the outcome. The authors developed a logistic regression model, which helped to identify 20 attributes that contribute the most to the mortality rate of pneumonia. To explore the anomaly detection methods, we have experimented with a simpler version of the PORT dataset that records, for every patient, only the attributes identified by Fine's study [START_REF] Fine | A prediction rule to identify low-risk patients with community-acquired pneumonia[END_REF]. The attributes are summarized in Figure 1. The output attribute corresponds to the hospitalization decision.

To study our anomaly detection methods in PORT dataset, we used 100 patient cases (out of a total of 2287 of cases). The cases picked for the study consisted of 21 cases that were found anomalous according to a simple Naïve Bayes detector (with detection threshold 0.05) that was trained on all cases in the database. The remaining 79 cases were selected randomly from the rest of the database. Each of the 100 cases was then evaluated independently by a panel of three physicians. The physicians were asked whether they agree with the hospitalization decision or not. Using panel's answers, the admission decision was labeled as anomalous when (1) at least two physicians disagreed with the actual admission decision that was taken for a given patient case or (2) all three indicated they were unsure (gray area) about the appropriateness of the management decision. Out of 100 cases, the panel judged 23 as anomalous hospitalization decisions; 77 patient cases were labeled as not being anomalous. The assessment of 100 cases by the panel represented the correct assessment of unusual hospitalization decisions.

HIT dataset

Heparin-induced thrombocytopenia (HIT) [START_REF] Warkentin | Heparininduced thrombocytopenia: recognition, treatment, and prevention: the seventh accp conference on antithrombotic and thrombolytic therapy[END_REF]) is a transient pro-thrombotic disorder induced by Heparin exposure with subsequent thrombocytopenia and associated thrombosis. HIT is a condition that is life-threatening if it is not detected and managed properly. The presence of HIT is tested by a special lab assay: Heparin Platelet factor 4 antibody (HPF4).

The HIT dataset used in our experiment was built from de-identified data selected from 4273 records of post-surgical cardiac patients treated at one of the University of Pittsburgh Medical Center (UPMC) teaching hospitals. The data for the was obtained with University of Pittsburgh Institutional Review Board approval. The data collected for patients was obtained from the MARS system, which serves as an archive for much of the data collected at UPMC. The records for individual patients included discharge records, demographics, all labs and tests (including standard and all special tests), two medication databases, and a financial charges database. For the purpose of this experiment the data were preprocessed and used to build a dataset of 45767 patient state examples for which the HPF4 test-order decision (order vs. no-order) was considered and evaluated. The patient states were generated automatically at discrete time points marked by the arrival of a new platelet result, a key feature used in the HIT detection. A total of 271 HPF4 orders were associated with these states (prior of a test order is 0.59%) Each data-point generated consisted of a total of 45 features that included recent platelets, platelet trends, platelet drops from nadir and the first platelet value, a set of similar values for hemoglobin and hemoglobin trends, whether a transfusion was done in last 48 hourse, an indicator of the ongoing Heparin treatment and the total time on Heparin.

To study the performance of our anomaly detection methods in the HIT dataset, we used 60 patient state cases (out of a total of 45767 of cases). The cases picked for the study consisted of 30 cases with the HPF4 order and 30 cases without HPF4. Each of these 60 cases was evaluated for appropriateness of HPF4 order by a pharmacy expert. 28 were found anomalous.

Experiments

All the experiments followed the leave-one-out scheme. That is, for each example in the dataset of patient cases (100 for PORT and 60 for HIT) evaluated by humans, we first learn the metric. Next, we identified the cases in E most similar to it with respect to that metric. The cases chosen were either the some number of closest cases (40 for PORT and 100 for HIT), or all the other cases (2286 for PORT or 45766 for HIT) in the dataset. We then learned the NB model or SVM and calculated the projection.

The target example was declared anomalous if its projection value fell below the detection threshold. The anomaly calls made by our algorithms were compared to the assessment of the panel and the resulting statistics (sensitivity, specificity) were calculated. To gain insight on the overall performance of each method we varied its detection threshold and calculated corresponding receiver operating characteristic (ROC). For the hospital deployment no all thresholds are acceptable. Consequently, for the evaluation we selected only that part of the ROC curve that corresponds to specificity equal or greater than 95%. The 95% specificity limit means that at most 1 in 20 normal cases analyzed may yield a false alarm.

Results

Tables 1 and2 show the ROC statistics for the feasible detection range. We see that for both datasets and models, using the NCA metric and selecting the closest patients outperformed all other methods (except for NB for PORT where it ended up second best). Moreover, in most of the other cases local models (using only close patients) achieved superior performance over their global counterparts. Close patients let us fit better the predictive model to the target patient, while taking all samples into the consideration biases the population. Regarding local models, performances of Naïve Bayes and SVM projections are comparable. For the HIT dataset we also show traditional area under ROC for the full specificity range and PR (precision-

Conclusions

Conditional anomaly detection is a promising methodology for detecting unusual events that may correspond to the medical errors or unusual outcomes.

We have proposed a new anomaly detection approach that uses the discriminative projection techniques to identify anomalies. The method generalizes previously proposed probabilistic anomaly detection framework [START_REF] Hauskrecht | Evidence-based anomaly detection[END_REF]. The advantage of the method is that it performs fully unsupervised and with the minimum input from the domain expert.

The new method was tested on the new Heparin induced thrombocytopenia dataset with over 40k patient state entries. The experiments demonstrated that our evidence-based anomaly detection methods can detect clinically important anomalies very well, with the detector based on the NB or SVM projections.

Despite initial encouraging results, our current approach can be further refined and extended. For example, instance-based (local) models tested in this paper always used a fixed number of 40 or 100 closest patients (or more, if the distances were the same). However, the patient's neighborhood and its size depend on the patient and data available in the database. We plan to address the problem by developing methods that are able to automatically identify and select only patients that are close enough for the case in hand.

Table 1 .

 1 Figure1. Attributes from the Pneumonia PORT dataset used in the anomaly detection study. PORT dataset: Area under the ROC curve in the feasible range of 95% -100% specificity. Please note that the baseline value for the random choice is 2.5%, maximum is 100 %.

		PORT dataset
	metric	model #cases	area
	any	NB	2286 11.6 %
	metric model #cases	area
	NCA	NB	40 16.8 %
	Mahalanobis	NB	40 17.6 %
	RCA	NB	40 17.6 %
	Euclidean	NB	40 16.4 %
	metric	model #cases	area
	any	SVM	2286 12.1 %
	metric model #cases	area
	NCA	SVM	40 19.0 %
	Mahalanobis	SVM	40 11.9 %
	RCA	SVM	40 10.4 %
	Euclidean	SVM	40 11.2 %

Table 2 .

 2 HIT dataset: Area under the ROC curve in the feasible range of 95% -100% specificity. Please note that the baseline value for the random choice is 2.5%, maximum is 100 %.

		HIT dataset
	metric	model AU-ROC AU-PR
	any	NB	57.8 %	50.9 %
	metric	model AU-ROC AU-PR
	NCA	NB	90.6 %	90.8 %
	Mahalanobis	NB	84.9 %	80.5 %
	RCA	NB	84.9 %	80.5 %
	Euclidean	NB	85.3 %	78.9 %
	metric	model AU-ROC AU-PR
	any	SVM	87.3 %	86.6 %
	metric	model AU-ROC AU-PR
	NCA	SVM	90.8 %	90.6 %
	Mahalanobis	SVM	87.6 %	82.9 %
	RCA	SVM	87.6 %	82.9 %
	Euclidean	SVM	90.4 %	90.8 %

Table 3 .

 3 HIT dataset: Area under the Receiver Operating Characteristic and Precision-Recall curves. recall) curve in table 3. The results in table 3 are qualitatively equivalent to those in table 2.

Acknowledgements

The research presented in this paper was funded by the grants R21-LM009102-01A1 and R01-LM06696 from the National Library of Medicine and the grant IIS-0325581 from National Science Foundation. The authors would like to thank Michael Fine who allowed us to use the PORT data.