R I Macphail

M A Courty

A Gebhardt

Soil micromorphological evidence of early agriculture in north-west Europe

Introduction

Various soil features, such as ard marks, spade marks, cultivation ridges, lynchets, and other colluvial deposits resulting from plough-wash, have been interpreted [START_REF] Fowler | Plough marks, lynchcts and early fields[END_REF]Courty et al. 1989: 130-3, pl. IVg) as evidence of early agriculture. Other environmental information (such as snails, charred seeds, pollen) has also been used to support such interpretations. Soil micromorphology is another potential tool, one which can give great detail and thereby clarify field evidence. Also soil micromorphology can help establish a past record of cultivation where none is obvious in the field, and can often place such evidence into broader sequences, which may include phases of clearance, occupation, abandonment, and the regeneration ofnatural vegetation cover. In this paper we introduce the scientific basis of soil micromorphology and show, through experimental approaches and the study of archaeological sites, how it can be used to investigate soils and ancient agriculture. Soil micromorphology, the microscopic study of soils in thin section, is based on the same principles as petrography. Unlike rocks or pottery, however, soils are unconsolidated and therefore it is necessary to take undisturbed samples carefully, sometimes supported by a metal box (Courty et al. 1989: 4&3), to ensure that the arrangement of soil materials, sand, silt, clay and organic matter etc., and the gas/water spaces (voids: Bullock et al. 1985: 43-5), are preserved. Oven or acetone 'dried' soils are impregnated by resins (Murphy 1986: 55-64), which eventually permits a cured block to be cut up, and the slices to be mounted on glass and made into thin sections some 20-30 p.m thick, and up to 15 x 7cm in size. Such large thin sections can give a very good representation of a buried soil, especially when numerous sampling columns are taken across an old land surface. Problems of soil heterogeneity can be tackled in this way (Macphail, in press). Soil micromorphology has developed its own complex descriptive terminology [START_REF] Kubiena | Micropedology[END_REF]Brewer 1964;[START_REF] Bullock | Handbookfor Soil Thin Section Description[END_REF]. Although used mainly in the study ofrecent soils, the technique has also been employed in the detection of ancient agriculture in soils buried beneath archaeological monuments (Cornwall 1958: 215;Manby 1963;Ashbee 1966;[START_REF] Romans | Soils and archaeology in Scotland[END_REF]1983). The interpretation of soil microfabrics is based on detailed pedogenic studies and on experiments on soils. Soil microfabrics are interpreted on a hierarchical basis, the effect of each soil process, including those relating to cultivation, being superimposed on, or disrupting earlier microfeatures [START_REF] Bullock | Evolution of a palaeo-argillic brown earth (Paleudalf) from Oxfordshire, England[END_REF]Fedoroff 1982: 295-8;Courty et al. 1989: l38-9).

Investigative approaches

Before any micromorphological evidence of ancient agriculture can be recognized, it is necessary to understand both the character of the natural soil and the precise changes which result from disturbances. Experimental studies of undisturbed soils, and the changes consequent upon cultivation, are therefore necessary. All human activity transforms the soil. The problem is to separate the effects of tillage from other forms of disruption, such as clearance of vegetation (see below) or trampling (Courty et al. 1989: 124). Modern agricultural science provides basic data on the differences between cultivated and non-cultivated soils, but other information can be obtained from simple experiments in 'ancient agriculture', mainly using non-mechanized tillage methods. The results from such experiments are, however, limited to the effects of generally short-term tillage. The effects of other transformations, relating to abandonment, revegetation and biological reworking, or burial, are not present. One major aim of the experiments reported here was to recognize features which could be used as indicators, and could be sought in archaeological soils. The aim was not, therefore, to reproduce exactly all the features which could be found in ancient cultivated soils; this would be impossible, in practical terms, since in all probability examples of undisturbed mid Holocene forest soils no longer exist. In fact, the full effect of tillage on 'ancient soils' has to be interpreted from the actual study of sites of prehistoric cultivation, and these also supply data on post-tillage soil transformations (Macphail, in press) and the effects of burial [START_REF] Allen | Micromorphology and magnetic susceptibility studies: their combined role in interpreting archaeological soils and sediments[END_REF].

Clearance features

Clearance of vegetation is commonly considered a necessary precursor to cultivation, although this, as shown below for the Hambacher Forest experiment, is not always essential. No experimental soils data on clearance by ancient techniques is available, but the results of studies of soils affected by 'blow-downs', or tree-throw caused by high winds (Lutz and Griswold 1939;[START_REF] Denny | Microrelief resulting from fallen trees[END_REF], and modern deforestation (Courty et al. 1989: 28&9O) are accessible. In addition, microfabrics of asymmetrical subsoil hollow infills in palaeosols, interpreted as resulting from the toppling of large trees in prehistory (Macphail 1986: fig. 2;1987;Macphail and Goldberg, in press) are also considered. The major features that result from this form of profile turbation are listed in Table 1. Ifa tree hollow is left to infill naturally, most of the microfabric heterogeneity may be homogenized by biological activity (Macphail 1987). Rapid burial by a monument, or through postclearance cultivation effects (burial and infilling by plough wash) may on the other hand allow this heterogeneity to be preserved, in the same way as pit and ditch fill microfabrics can be protected by dumping and sedimentation. Where young secondary woodland or shrubs are cleared, the impact on the soil profile is likely to be far less and therefore less readily recognizable in the field. The resulting heterogeneous microfabrics, characterized by dusty clay infills (Bullock ct al. 1985: 96), are, however, typical of disturbed soils (Macphail and Goldberg, in press), particularly when associated with indicators of vegetation burning (abundant charcoal fragments and phytoliths; Table 1). The in situ burning of toppled trees may permit areas of soil that are sufficiently strongly burned to be preserved as coarse red soil fragments. The presence of coarse wood charcoal in the surface horizons of a Neolithic soil buried by a barrow at Daladies, Angus, Scotland, led other workers [START_REF] Romans | Soils and archaeology in Scotland[END_REF] to infer 'slash and burn' clearance. Spontaneous fires could provide natural soil charcoal, but the natural burning of deciduous trees in the mid-to late Holocene of north-west Europe is thought unlikely.

The effects of cultivation on soil microfabrics

Much of our knowledge of the mieromorphology of cultivated soils comes from the work of [START_REF] Jongerius | Some morphological aspects of regrouping phenomena in Dutch soils[END_REF][START_REF] Jongerius | Micromorphology in agriculture[END_REF] and Kooistra (1987), Some of the basic soil microfabric changes which occur in natural soils after they are cultivated are given in Table 2. Many of the large-scale changes and typical features occurring in modern cultivated soils relate to the effects of pressure and shear (and weight and vibration) caused by modern machinery. Such changes cannot be directly compared to those brought about by early cultivation. Nevertheless, those phenomena in modern cultivated soils which arise from the fact that, in comparison to natural vegetated soils, cultivated ones remain bare, or are irregularly protected, are the same. These phenomena are the breaking up and loss of structure of the surface soil aggregates under rain splash impact, and the consequent slaking and mobilization of water saturated soil. This is followed by washing and sedimentation; holes between tillage clods are infilled, and surface crusts (Plsl and 2) form (sometimes even between furrows), which may contain strongly oriented particles [START_REF] Jongerius | Micromorphology in agriculture[END_REF][START_REF] Boiffin | Dynamique de formation des routes superficielles: rapport de l'analyse microscopique[END_REF]. Trapped air often leads to the formation of vesicles in fully saturated horizons. These effects are also linked to changes in structural stability, in part governed by organic matter content, when forest and grassland are converted to arable use (Table 2; Imeson and Jungerius 1976;[START_REF] Grieve | Some contrasts in soil development between grassland and deciduous woodland sites[END_REF]. Soils that are puddlecl (Jongerius 1983: 113) through a combination of high moisture contents and trampling (e.g. by cattle) may also develop similar fabrics [START_REF] Beckman | Micromorphological changes in surface soils following wetting, drying and trampling[END_REF]Courty et al, 1989: 130). The pressure and shear forces produced by ancient tillage in soils where the moisture content is sufficiently high can also produce features indicating puddling and slaking (Pl, 3). Trapped water within the plough zone may also cause internal (subsurface) slaking. The soil separates and particles move downwards to become sorted or unsorted textural pedofeatures (Bullock et al. 1985: 109-115). [START_REF] Jongerius | Some morphological aspects of regrouping phenomena in Dutch soils[END_REF] coined the term 'agricutan" for these coatings and infills that form down-profile as the result of surface or internal slaking. He noted that in light textured soils, these agricutans are mostly built up ofvarious silt and clay and impure clay (Bullock et al. 1985: 111, 114) layers, because of variations in sorting when all the components and size fractions of the ploughed horizon (Ap) are in movement after tillage disturbance. In heavier textured soils the coatings and infillings of voids are homogeneous, and principally made up of clay and microscopic fragments of organic matter and opaque material; hence the term 'dusty clay' (Pls 4, 5 and 6). In contrast, 'limpid clay' coatings (Bullock et al. 1985: 111, 114) are more generally associated with weathering, and clay translocation in alfisols (Soil Survey Staff 1975) under undisturbed broadleaved woodland [START_REF] Fedoroff | Comparative micromorphology of two late Pleistocene paleosols (in the Paris Basin)[END_REF]. In undisturbed alfisols these coatings are generally absent from surface horizons because of the homogenizing effect of biological activity, but they accumulate in the subsoil. Tillage also leads to the mixing of soil material and to the washing of pieces of plough soil down cracks in the profile [START_REF] Jongerius | Some morphological aspects of regrouping phenomena in Dutch soils[END_REF]1983). The physical action of tillage results in smaller aggregates accumulating, and void space diminishing, at the base of the plough zone (see Pl. 12), whereas large aggregates and coarse voids concentrate near the surface [START_REF] Dexter | Internal structure of tilled soil[END_REF]1979).

Experiments in 'ancient agriculture'

The microfabrics of four uncultivated soils were compared with those that developed after cultivation by hoe, ard, spade and motor-cultivator (Gebhardt, in press). Exact replication of ancient tillage methods, crops, or soil types was not possible, but as stated earlier, this was not believed to be critical for these experiments. Such experiments, however, cannot be carried out just anywhere, and so two sites with well-documented land use and soil characteristics at the modern experimental farm of Grignon, Yvelines, France (lnstitut National Agronomique), namely the Grignon forest and Dehérain plots [START_REF] Morel | Les essais defertilisation de longue durée de la station agronornique de Grignon[END_REF], and two 'ancient farms' (Butser Hill, Hampshire, England: Reynolds 1981; and Hambacher Forst, Elsdorf, Rhinelands, Germany: Lüning and Meurers-Balke 1980; Meurers-Balke 1985) were selected (Table 3). Thus in this preliminary series of soil micromorphological experiments in 'ancient agriculture', both the parent materials and soil types (cambisols and alfisols: Soil Survey Staff I975) at the 'ancient farms' were replicated at a soil science establishment. Sampling for micromorphology took place either immediately or a few days after tillage, using three or four tin monoliths and extending to depths of 30-50cm. Virgin or uncultivated surface horizons, and sometimes their subsoils (e.g. at Hambacher Forest), were similarly sampled.

The uncultivated soils are rather compact with good subangular blocky structures, except for the Hambacher soil which was very weakly structured because of its fine sandy/coarse silty texture. The Deherain soil is the most open structured (Pl. 7). All, however, are characterized by large amounts of included organic matter, and biological activity, which are responsible for the intra-pedal porosity. In all cases cultivation led to the soil becoming looser, through an 'explosion' of fine structures (Pl. 8) which relates to inherent weaknesses in soil peds caused by wetting and drying. These weaknesses fractured under cultivation impact, a phenomenon also known from modern agriculture. The size of the newly-formed structures reflects the nature of the cultivation implement: medium angular peds, Butser (ard); coarse angular peds, Deherain (spade). Cultivation by motor cultivator at Dehérain also led to the smaller peds concentrating at around 10-20cm depth to produce a weakly compacted layer. Manured cultivation soils at Butser and Dehérain showed enhanced structural stability and biological activity. As a result, continuing good structural development at Butser, even after tillage, made it difficult to ascertain the depth limit of ard cultivation, but the more coarse effects of spade cultivation at Dehérain was established down to 12cm. No depth limit to hoeing was recognized at Grignon Forest, because of structural instability (see below). At all the sites the nature of the uncultivated soil could often be recognized as humic soil fragments within the biggest clods, but the increase in porosity and numbers of structures is far more poorly expressed in the lighter textured soils at Grignon and Hambacher Forest.

For example, hoeing at Grignon Forest was accompanied by intensive watering, which eventually resulted in structural collapse. At Hambacher Forest the cultivated soil contains numerous rounded peds of humic topsoil and more mineral subsoil. In these weakly structured soils, the peds were apparently compacted into the loose soil matrix by ard tillage. At around 6cm depth, at the base of the tilled soil at Hambacher Forest, a lcm thick weakly-compacted silt pan formed, and at about 12-14cm depth some dusty clay accumulated. The latter was superimposed over limpid clay coatings that characterize the uncultivated 'forest' soil.

Data from ancient cultivated soils

Only a few archaeological sites with independent (non-pedological) evidence of cultivation have been investigated by soil micromorphology. The examples cited (Table 4) provide the best range available of microfabrics resulting from cultivation activity, varying according to topography, parent material and soil type. Soil texture strongly influences the fabric type, structure and tcxtural features that result from cultivation [START_REF] Jongerius | Some morphological aspects of regrouping phenomena in Dutch soils[END_REF], and so the examples are ordered from the coarsest soils to the finest. In the examples, the ancient cultivated soils were well sealed by a number of agencies, both natural (alluvium at Phoenix Wharf and blown sand at Lodbjerg) and human (colluvium at Ashcombe, walls at Le Wasserwald and a long barrow at Kilham). Any investigation of ancient soils should ascertain the degree of preservation, and soil micromorphology is probably the best technique for making this estimation. For example, biological reworking by earthworm or root penetration through the overburden, can be readily recognized, and alteration of organic matter and some other bio-chemical effects can also be identified (Macphail 1986: 269-70;[START_REF] Allen | Micromorphology and magnetic susceptibility studies: their combined role in interpreting archaeological soils and sediments[END_REF]. In fact, sampling upwards into the colluvium or mound above a buried landsurface is common practice; for instance, to make sure that clay coatings within the surface of a buried soil have not come from the construction of an archaeological structure [START_REF] Romans | The general effects of early agriculture on the soil[END_REF]. The identification of post-cultivation events, such as abandonment or occupation before burial, can be made through the analysis of microfeatures on an hierarchical basis. At Phoenix Wharf, London (Table 4), where ardmarks were interpreted, the buried plough soil comprises sand grains coated by a fine fabric of silt and organic matter (see Table 2; Pls 9 and 10). The underlying sand contains downwashed fine material, but because of the total absence of clay this subsoil has only rare and poorly formed coatings. On the more coarse windblown sands at Lodbjerg, west coast of Jutland [START_REF] Courty | Comparison between buried uncultivated and cultivated Iron Age soils on the west coast of Jutland, Denmark[END_REF], tillage (ardmarks) produced an Ap horizon similar to that of Phoenix Wharf, but no subsoil features. There is good organic matter preservation of the fine fabric at both sites through waterlogging (Phoenix Wharf) and acidity (Lodbjerg). This may help account for the well developed Ap horizons, which are much thicker than the 6cm depth impact of ard tillage. Also, plough soil accretion has been enhanced at the sites by alluviation and sand blowing, respectively, and the additional presence of much domestic waste may indicate manuring of the cultivation soils as well. Pollen analysis of the soil buried under the cairn at Chysauster, Cornwall (Scaife pers. comm.) shows a phase of woodland clearance, followed by cultivation and soil acidification. The turbated subsoil microfabric (Table 4) is good evidence of woodland clearance. This clearance probably reduced both the earthworm population and the stability of the topsoil and, with cultivation, may have led to surface soil acidification. The later Iron Age field systems at Chysauster are characterized by well developed lynchets, and within these colluvial deposits are a number of fragmented crusts and coatings, consistent with surface and internal slaking in cultivated soils [START_REF] Jongerius | Some morphological aspects of regrouping phenomena in Dutch soils[END_REF][START_REF] Jongerius | Micromorphology in agriculture[END_REF]); these indicators have been preserved by podzolic cementation. The colluvium also contains rounded soil fragments, both of the relict brown subsoil and the contemporary podzol cover, as evidence of the widespread and deep erosive effects of long term cultivation. Sampling of the soil outside the ancient field system at Le Wasserwald (Bas-Rhin, France) found only limpid clay void coatings, indicating the earlier presence of a forest cover in the area. In contrast, Gallo-Roman cultivation, as evidenced by plough marks cutting into the sandstone parent material, seems to have led to an accumulation of fine fabric at the soil surface. Translocation of coarse soil, similar to agricutans [START_REF] Jongerius | Some morphological aspects of regrouping phenomena in Dutch soils[END_REF]) was confined to the cultivated area, and was found to be superimposed over the limpid clay coatings. On the somewhat finer terrace alluvium at Strathallan, Perthshire, local pollen and field features ('cultivation ridges') suggested Neolithic cultivation [START_REF] Barclay | Sites of the third millennium BC to the first millennium AD at North Mains, Strathallan, Perthshire[END_REF]). Semi-quantitative numerical analysis of the vertical distribution of textural features was carried out on soils beneath the henge, and through the mound (see Pl. 4) and buried soil of the barrow site, including samples through the 'cultivation ridges' [START_REF] Romans | The general effects of early agriculture on the soil[END_REF]. The highest concentrations of textural pedofeatures were found near the soil surface (4-5cm) of the cultivation ridges under the barrow mound, whereas under the henge the greatest concentration of textural features was around 12cm depth. [START_REF] Romans | The general effects of early agriculture on the soil[END_REF] thought that the greatest concentration of textural features would reflect the soil depth most disturbed by shear and pressure forces produced by the impact of the tillage implement; they inferred, therefore, that there had been shallow hoeing of the 'cultivation ridges' at the mound site, and deeper ard tillage in the area of the henge. The lack of earthworm fabrics and the low organic status of the buried soils (even after post-burial oxidation was taken into account) indicated to Romans and Robertson that Neolithic cultivation had led to a diminution of soil organic matter, resulting in poorly stable soils (Courty et al. 1989: plate IVe). A land-use history comprising two cultivation phases separated by a period of stable grassland was inferred from soil pollen in the shallow (20cm) soil under the Neolithic long barrow at Kilham, Yorkshire [START_REF] Dimbleby | Pollen and land snail analysis of calcareous soils[END_REF]. The microfabric comprises a partially homogenized soil, with root channels coated by dusty clay (see Pl 5 and 6). These features may reflect either a later land-use history of grassland followed by cultivation (as inferred from the pollen), or that the root channels were the result of crop growth and were coated by the effects of continuing tillage. On similar, but colluvial, silt loam parent materials in the dry valley at Ashcombe Bottom, Sussex, Beaker period cultivation was apparently recorded by ardmarks [START_REF] Allen | Ashcombe Bottom Excavation AB84[END_REF]. The plough soil colluvium contains rounded soil fragments eroded from the relict alfisol cover that had formed from the weathered chalk and 'Clay-with-Flints'. Also pieces of clay rich crust and associated charcoal (see Plates 1 and 2), presumably resulting from the surface slaking of alfisol material in cultivated soils, were mixed into the Ap horizon by later tillage and colluviation. Some compact soil areas of the upper plough soil were characterized by planar voids (shear planes: see Pl. 3), presumably from 'plough' shear. In addition, the dense prismatic structures were broken up in places into finer subangular blocky structures by the impact of tillage, and biological activity has been concentrated in these more open soil areas (Pl. 11). In contrast, towards the base of the plough soil, the soil is much more compact, with textural features (Plate 12) indicating that the resulting poorer drainage encouraged internal slaking.

Discussion

The identification and interpretation of ancient cultivated soils requires the differentiation of precultivation features from those generated by the phase or phases of agriculture itself, and those caused by later events such as abandonment and burial. In the analysis of soil thin sections, relict late glacial features may have to be considered separately before early Holocene forest soil (alfisol) formation (e.g the limpid clay at Le Wasserwald) can be interpreted. Similarly, the turbating effects on the soil of shrub, scrub or forest clearance will need to be discriminated from later mixing by cultivation. For example, at Neolithic Maiden Castle, Dorset, and Beaker period Brean Down, Somerset (Macphail and Goldberg in press), vertical sampling of buried soils produced evidence of soil disruption at the base of the soil profiles which pre-dated microfabrics consistent with the effects of cultivation. This succession of features clearly indicated that woodland clearance preceded agricultural activity. Some workers (Liming and Meurers-Balke 1980;[START_REF] Meurers-Balke | Experiments zum Anbau und zur Verarbciterung priihistoriseher Getreidcarten[END_REF]Romans, pers. comm), however, have suggested that Neolithic farmers did not necessarily clear the land of trees, but ploughed between them. It is evident that soils will respond differently to tillage according to their texture, base status, original organic matter content and levels of biological activity. As predicted by [START_REF] Jongerius | Some morphological aspects of regrouping phenomena in Dutch soils[END_REF]1983), sands will lose their weak structures and these will be replaced by various intergranular fine fabrics; such effects are seen in the buried soils at Lodbjerg, Phoenix Wharf (Pls 9 and 10) and Le Wasserwald. The low pedality of silt loam soils can lead to loss of structure (as at Grignon Forest) and to weakly compacted pan formation (as at Hambacher Forest). The finer, high base status soils at Butser and Dehérain reacted to tillage by an explosion offine structures (Pls 7 and 8), that were generally stable, especially where manuring had enhanced levels of organic matter and biological activity. This break-up of coarse structures by tillage was found at Ashcombe (Pl.11), as was the movement down profile of fine soil peds (as noted at Dehérain), although the poorer structural stability at Ashcombe led to compaction (Pl. I2). At Strathallan the soils readily slaked (Pl. 4) under the impact of tillage, with supposed ard cultivation causing the concentration of textural features at l2cm [START_REF] Romans | The general effects of early agriculture on the soil[END_REF]. The actual depth limit of ard tillage at Hambacher Forest was only 6cm, but dusty clay coatings resulting from this cultivation were accumulating at l2-14cm. The abundant presence of dusty clay coatings in voids in upper soil horizons beneath the Neolithic long cairn at Hazleton, Gloucestershire, (Saville 1984;Macphail etal. I987: figs. 7 and 8), has been used to infer cultivation on this site. As at Strathallan, the impact of tillage on soils with poor structural stability is thought to have caused slaking and the development of textural features. One reason for this poor structural stability could be lack of organic matter. At both sites earthworm activity, which is abundant in the surface horizons of uncultivated forest soils, is very low, and may be another reflection of the low organic matter levels that could have resulted from early cultivation practices. At Phoenix Wharf and Lodbjerg, manuring was inferred, and at Hazleton it was noted that the final phase of cultivation either moved to a domestic refuse-rich midden area (Macphail, in press), or that manuring was undertaken. Beneath the immediate surface horizons of the buried soils at Maiden Castle, Hazleton and Kilham, biological features (peds or root channels) were found to be coated with dusty clay (Pls 5 and 6), indicating that during ancient agriculture, biological activity (e.g. rooting by crops/weeds; faunal working) could accompany the effects of surface soil slaking and down profile soil translocation, induced by tillage. At Chysauster, on the other hand, clearance and tillage led to a decrease in the earthworm population through acidification. Lastly, changes in the soil landscape may also be inferred from eroded soil fragments in colluvial ploughsoil deposits, as at Ashcombe Bottom and Chysauster. Obviously, the best examples of ancient cultivated soils are those buried soon after being used for agriculture. If the cultivated soil is abandoned but, before burial, an herbaceous cover develops, or shrubs invade (Macphail in press), or the land is given over to pasture, earthworms will eventually rework it. In addition, tilled soils can be trampled and mixed by domestic occupation, or partially eroded. All these post~cultivation events blur or totally remove the agricultural history from the soil. Post-burial effects have also to be considered (Courty et al. 1989: 140-6).

Conclusion

In the study of ancient cultivated soils, microfabrics have to be studied on an hierarchical basis, so that agricultural features can be distinguished from those of the natural soil. Modern agricultural science has provided a data base for these investigations, but experiments in 'ancient agriculture', especially when combined with empirical studies of ancient sites, have been important in permitting the mechanisms involved in tillage to be better understood. Significantly, this success has occurred in spite of experimental problems with replicating the exact tillage method, crop or soil type which existed in the past. Now that some of the physical mechanisms are understood, it is possible to identify microscopic fabric, structural and textural indicators of ancient agriculture, although these vary according to soil type, grain size, organic matter content and base status. It is the combination of indicators that is important, so that, for example, 'agricutans' should not be regarded on their own as categorical evidence of cultivation.

Acknowledgements

The authors are very grateful for an 'Action Thematique Programmee' grant from the CNRS, France; the support of English Heritage; and the co-operation of the many directors of the experimental farms and archaeological sites discussed in this study.

Abstract

Macphail, R. I., Courty, M. A., and Gebhardt, A. Soil micromorphological evidence of early agriculture in north-west Europe The study of early agricultural soils using micromorphology is outlined. The effects of tillage due to modern agriculture are reviewed and the results from four experiments in 'prehistoric cultivation' are presented. These are used in the interpretation of prehistoric agricultural practices, based on a study of soils from various archaeological sites in north-west Europe. It is suggested that microfabric features, varying according to soil texture, can be diagnostic of tillage. Microfabric type, textural features (e.g. coatings) and structure are, taken together, better indicators of ancient tillage than single features.


Raabid Hussain

Caroline Guigou

Berihu Kibrom

Alain Girum

Alexis Bozorg Lalande

Grayeli

Grayeli Bozorg

3d

3D landmark detection for augmented reality based otologic procedures

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Ear consists of the smallest bones in the human body and does not contain significant amount of distinct landmark points that may be used to register a preoperative CT-scan with the surgical video in an augmented reality framework. Learning based algorithms may be used to help the surgeons to identify landmark points. This paper presents a convolutional neural network approach to landmark detection in preoperative ear CT images and then discusses an augmented reality system that can be used to visualize the cochlear axis on an otologic surgical video.

Introduction

Augmented reality (AR) has been widely accepted and applied in different surgical domains like orthopedics, hepatobiliary and pancreatic systems. However, owing to minuscule operative space, limited field of view and complex instrument trajectory, AR has not been widely accepted for otologic procedures. Otologic structures constitute some of the smallest structures in the human body and their handlings require submillimetric precision and expert knowledge of intra and inter anatomic relations as critical nerves and vessels are in close proximity.

Image registration is one of the critical processes in an AR system. In otology, preoperative computed tomography (CT) is particularly helpful as the ear is composed of mostly rigid bony structures. Precise identification of anatomical landmarks for registration is often difficult and time consuming as the points are not well defined. Due to low number of prominent features between preoperative CT reconstructions and endoscopic videos, different studies have used fiducial markers to register the pre and intra-operative images [START_REF] Liu | Cadaveric feasibility study of da vinci si-assisted cochlear implant with augmented visual navigation for otologic surgery[END_REF][START_REF] Hussain | Real-Time Augmented Reality for Ear Surgery[END_REF].

This study aims at simultaneously determining locations of multiple anatomical landmarks in preoperative CT-scan using a learning based approach. The detected landmarks may then be used to register CT with microscopic images. Typical learning models employ Gaussian heat maps to optimally determine the landmark position [START_REF] Payer | Regressing heatmaps for multiple landmark localization using CNNs[END_REF]. However, since the distance between each landmark is very small, this is not a viable option in this scenario. Moreover, state-of-the-art learning architectures often require manual fine-tuning of meta-parameters which often leads to loss of important information. This study proposes and evaluates the performance of a convolutional neural network (CNN) for middle and inner ear landmark detection without standardizing the meta-parameters and tests them on an AR scenario to infer the cochlear axis.

Methodology

Dataset

We used a dataset of 25 patients (age range: 15-77 years) comprising of 40 ear CTs (17 right and 23 left ears) acquired from different scanners. The pixel resolutions of the acquired CT scans ranged from 0.156x0.156x0.100 mm 3 to 0.292x0.292x0.625 mm 3 . The x-y image size was the same for all the images: 512 x 512 while the z image size ranged from 92 to 876. Locations of seven anatomical landmarks were identified in all the scans by an expert surgeon:

3D landmark detection for augmented reality based otologic procedures

Raabid HUSSAIN, Caroline GUIGOU, Kibrom Berihu GIRUM, Alain LALANDE and Alexis BOZORG GRAYELI

ImViA Laboratory, Université de Bourgogne Franche-Comté, 21000 Dijon, France Tel : +33 (0)6 95 96 55 40 Contact: raabid.hussain@u-bourgogne.fr round window niche, tip of the incus, umbo and short process of malleus, pyramid, cochlear apex and base.

Landmark detection

To homogenize the data, the left ear CTs were flipped to resemble the right ear CTs. Due to memory limitations, the surgeon was asked to crop a 200x200x100 region of interest from the original CT data that contained the middle and inner ear contents. The region of interest was passed through a convolutional neural network depicted in Figure 1. All layers had 'elu' activations except the final layer which had linear activation. The output layer had 21 units comprising of x, y and z coordinates of each landmark. The training was carried out for 3500 epochs with a batch size of 5, using Adam optimizer with a learning rate of 0.0005, and mean squared logarithmic error as the loss function. The architecture was implemented on a computer with dedicated GPU (NVIDIA TITAN X, 12 GB RAM processor) using Keras and Tensorflow libraries. The network was trained from scratch and assessed using 5 fold cross validation approach.

AR system

An AR system was developed to test the output of the landmark detection step. Cochlear base point was used as a test landmark as it is obscured in microscopic view. Remaining six landmarks were used to register the preoperative CT with the realtime microscopic video. The surgeon was asked to select their corresponding points in the initial microscope image. Fundamental matrix of the microscope camera (determined through the 2D points in microscope image and their respective 3D coordinates in CT image) was used to register the two information, followed by the use of a SURF based feature matching to track the motion of the microscope (similar system as used in [START_REF] Marroquin | Augmented reality of the middle ear combining otoendoscopy and temporal bone computed tomography[END_REF]). The cochlear axis (defined by cochlear apex and base points) was drawn onto the microscopic video for the surgeon to have more insight about the inner ear.

The system was tested on a phantom resin model of the temporal bone with the tympanic membrane removed.

Results

The landmark detection on CT-scan yielded an average error of 0.88 ± 0.27 mm (mean absolute distance between position of the detected landmark and its location as indicated by surgeon). Individual errors for each landmark are listed in Table 1. The AR system is displayed in Figure 2. An ENT expert verified that the position of the cochlear axis was indeed in close proximity to where the axis was expected. The system maintained correspondence throughout the experimental time of 2 minutes.

Conclusion

This paper presented a CNN based otologic landmark detection and its potential application in AR based surgery. The AR system showed promising results, compatible with otologic requirements, but also highlighted the need for extending the system to 3D for better visualization and ergonomics.

Figure 1 :

 1 Figure 1: CNN for landmark detection. I=Input layer, C=Convolutional layer, SE=Squeeze and excitation block [4], P=Max pooling layer, FC=Fully-connected layer, D=Dropout with a rate of 0.2, O= Output layer, f=number of filters.

Figure 2 :

 2 Figure 2: Augmented reality display. The cochlear axis is shown with a red line and the landmark points are represented by yellow dots.

Table 1

 1

		: Landmark detection process accuracy. R=Round
	window niche, I=Incus tip, U=Umbo of malleus, S=Short
	process of malleus, P=Pyramid tip, A=Cochlear apex,
			B=Cochlear base
		R	I	U	S	P	A	B
	Error (mm)	0.82 ± 0.36	0.71 ± 0.22	0.88 ± 0.26	0.92 ± 0.31	0.89 ± 0.27	0.93 ± 0.16	1.05 ± 0.31

Acknowledgements

The authors would like to thank Oticon Medical, France for their financial support and NVIDIA for donating the TITAN X processor under their GPU grant program.


Marc Liévin
email: lievin@caesar.de

Franck Luthon
email: franck.luthon@univ-pau.fr

Erwin Keeve
email: keeve@caesar.de

Entropic Estimation of Noise for Medical Volume Restoration

  L'archive ouverte pluridisciplinaire

Introduction

Image restoration techniques are widely used in the field of medical imaging, where human interpretation is linked with interactivity for challenging anatomical segmentation. Medical images are actually strongly deteriorated by noise mainly due to data acquisition systems. While signal-tonoise ratio (SNR) is improving in such scanning systems, higher resolutions are still required in surgical applications. As a main consequence, the noise level remains high for an elementary volumetric data, i.e. voxel.

Several methods coexist in the medical literature, each satisfying a subset of requirements. Our aim is to propose a restoration technique that is meant to be used by subsequent algorithms for medical image processing, namely connectivity filters. Therefore, the restoration should be unsupervised to propose interactivity for non-expert users while keeping relevant anatomical information.

Advanced medical restoration techniques can be categorized in two major groups. The first one concerns optimal filters, adapting to local and spatial information in images. For instance, the Wiener filter can be combined with a visibility function depending on the magnitude of the image gradient vector [START_REF] Abramatic | Nonlinear restoration of noisy images[END_REF]. Recently, the adaptive theory has been extended to 3-D space and collections of filters [START_REF] Westin | Three-dimensional adaptive filtering in magnetic resonance angiography[END_REF]. Non-linear wavelet filtering has proved to be useful for specific applications (e.g. mammography enhancement) [START_REF] Donoho | Nonlinear wavelet methods for recovery of signals, densities and spectra from indirect and noisy data[END_REF]. Unfortunately, the shape of the filters depends on the nature of the corrupting noise. To make a compromise between stability and smoothness, this parameter is generally manually tuned. The second category proposes nonlinear statistical estimation [START_REF] Hamza | Image denoising: A nonlinear robust statistical approach[END_REF]. Restoration may be efficient by using Markov random fields but the sensitivity to structures distributions is a major limitation [START_REF] Held | Markov random field segmentation of brain MR images[END_REF]. Moreover, such Bayesian approaches are very dependant on the normalization level parameters. Simulated annealing has been applied on the parameters estimation to reduce its sensitivity [START_REF] Aykroyd | Bayesian estimation for homogeneous and inhomogeneous Gaussian random fields[END_REF]. To our knowledge, few methods for estimating the amount of noise in images have been proposed [START_REF] Olsen | Noise variance estimation in images[END_REF]. Prior knowledges and beforehand determined parameters are the common solution chosen for medical image restoration.

Two methods based on absolute difference are tested in this paper. Image difference can be considered as classic method in image registration or noise estimation [START_REF] Olsen | Noise variance estimation in images[END_REF]. Two measurements are then applied on these observations: the empirical standard deviation estimation and an entropic deviation computation. On the opposite to generic techniques that use entropy as a fitting criterion, this measure is used to estimate the power of noise. Next, a dedicated Markov random fields restoration process is presented, taking account of the entropy and volumetric constraint potentials. Finally, the robustness of the results is evaluated on synthetic data and real datasets representing several medical modalities.

Noise Deviation Estimation

Noise Estimation

We model the medical volume data V as a collection of slices S i , corrupted by an independant additive Gaussian noise G(0 ). This rough assumption respects however the classic model of image backprojection acquisition and often corresponds to a good approximation in actual cases [START_REF] Nowak | Wavelet-based Rician noise removal for magnetic resonance imaging[END_REF]. Tomographic images are actually corrupted by a Rician

Liévin et al. Entropic Estimation of Noise for Medical Volume Restoration

noise, result of the inverse Radon transform, that can be approximated by Gaussian distributions in noisy environment.

The two estimations tested are defined as:

1. The difference between the slice S i and its Gaussian filtered image rS i , sf i = jS i ; r S i j. This observation is therefore directly related to the noise and the spatial edges of the slice. The variance is generally under-estimated.

2. The slice difference sd i = jS i+1 ; S i j, that contains information on the noise G(0 p 2 ) and on the vertical gradient of anatomical objects of the volume. The variance is generally over-estimated.

The empirical standard deviation std is a generic measure of the power of Gaussian distribution fx i g but sensitive to outliers:

std 2 = 1 N;1 P N i=1 (x i ; x) 2 :
The following paragraph introduces a new estimation of the noise deviation.

Entropic Standard Deviation

As defined by Shannon in [START_REF] Shannon | A mathematical theory of communication[END_REF], the entropy E(X) of a discrete source X of observations taking values in the set fi i 2 0 M ; 1]g is a measure of its mean information:

E(X) = ; M ;1 X i=0 p i log(p i ) (1)
where p i is the probability of the event i.

For an arbitrary source with a given entropy E(X), Shannon also defines the notion of entropy power N(X), which represents the power of the white noise equivalent to source X, in the sense it has the same entropy and is limited to the same band:

N(X) = 1
2 e e 2:E(X) :

(

The two important properties of the entropy power are:

The entropy power of the sum of two signals is lowerbounded by the sum of their respective entropy powers, and upper-bounded by the sum of their actual average powers.

The Gaussian source has the maximum entropy for a given power, verifying E(G) = l o g ( p 2 e). Hence, its entropy power equals its power: N(G) = 2 . As a consequence, the entropy power of any source is always less than or equal to its actual power.

According to our additive model, the noise estimation n (i.e. sd or sf) is the sum of a relevant signal X plus an additive Gaussian noise G(0 ). The first property yields the following inequality, where P X is the average power of the gradient signal:

N(X) + N(G) N(n) P X + 2 :

(

Applying the second property, this inequality becomes:

N(X) N(n) ; 2 P X :

White Gaussian noise has the peculiar property that it can absorb any other signal which may be added to it, provided the signal power is small, in a certain sense, compared to noise [START_REF] Shannon | A mathematical theory of communication[END_REF]. The resultant entropy power is therefore approximately equal to the sum of the white noise power and the signal power. Considering that relevant changes between slices are localized on the object frontiers in the image, whereas noise is spread over the entire frame, the average power P X of the difference source is small compared to the noise power 2 (i.e., the average SNR over the difference between slices is small, although the local SNR may be high). Therefore, we have:

N(n) PX 2 P X + 2 (5)
In summary, the entropy power is a measure that provides a means to evaluate the power of the noise in the slices. We call entropic standard deviation i of the noise observation n i :

i = p N(n i ) = 1 p 2 e e E(ni) : (6)

Noise Deviation Estimation in Medical Tomographic Images

Four tomographic volumes in magnetic resonance (MR) and X-ray computed tomography (CT) modalities have been evaluated. These medical volumes include 2 different modalities and 3 anatomical structures, i.e. head, knee and pelvis. One synthetic data have been tested respecting medical data level properties 1 . Table 1 details the cases studied in this paper.

Table 2 gives the results of the noise deviation estimation methods (std, ) when applied on the image difference observations (sd,sf ). Manual estimation has been made by manual determination of a subvolume of interest with no anatomical structures and homogeneous background, the deviation corresponds to the geometrical mean of the proposed methods. The synthetic data (case 5) reveals that the 1 The synthetic data coded on 12 bits consists in a sphere of value 2048 on a background of 256, corrupted by an additive noise G(0 50 2 ).

Case Modality Anatomy Size (x,y,slices) empirical deviation is very sensitive to gradients whereas the entropic measure stays close to the actual noise level, = 50. CT datasets (cases 1, 2) are known to be much less corrupted by noise than MR datasets (cases 3). Only the entropic estimation respects this difference in a ratio of 10 whereas the empirical deviation is less convincing. These results show that the entropic deviation when combined with slices difference is an accurate and reliable measurement of the noise level.

Constancy of the estimation along the slices is mandatory for a reliable volume restoration. Table 2 details between braces values of the variance in a z-score form, i.e. the distribution is normalized by the ratio between the variance and the mean value. This criteria reveals the "stability" of the methods along the volume slices. The empirical estimation appears to have the lowest stability (higher values), e.g. the difference with the entropic stability measurement is at minimum 20% and can reached 50%. The entropic deviation is therefore a better criteria when no manual tuning is allowed: the estimate is little biased by anatomy variations and its repetability can be guaranteed.

Volumetric Restoration

Volumetric Markov random fields

To confirm the quality and the stability of the entropic estimation, we integrate this measure in a restoration framework based on volumetric Markov random fields (MRF). The label field is supposed to follow the Gibbs property related to a neighborhood structure, i.e. the label of a voxel depends only on the labels of its neighbors and associated probabilities are modeled as Gibbs functions. In Bayesian context, Maximizing the A Posteriori probability (MAP criterion) of the label field is equivalent to minimize a global energy function [START_REF] Geman | Stochastic relaxation, Gibbs distributions and the Bayesian restoration of images[END_REF]:

U = U l + :U p
where U l represents the likelihood energy expressing the link between labels and observations, U p defines the prior energy corresponding to spatial constraints, is a weighting coefficient between the two energies.

In the following, v(x y z) is a voxel of the volume V , n denotes a neighbor of v belonging to its neigborhood . The restoration purpose is to find the closest label l v corresponding to the observation o v of the original data. As in the rest of the paper, we consider an additive and centered Gaussian noise with a standard deviation .

The two energies are expressed as:

U l (V ) = P v2V h (ov;lv ) 2 2 2 i U p (V ) = P v2V h P n2 (v) V p (v n) i (7)
where V p defines spatial prior constraints. We propose to construct the metric prior functions V p as the inverse of the Euclidian distance between two neighbors. These potentials constraint the field towards homogeneity respecting the data spacing defined in millimeters ( x y z ):

V p (v n) = ( ; V jj ; ! v ; ; ! n jj if l v = l n + V jj ; ! v ; ; ! n jj otherwise.
where jj ; ! v ; ; ! n jj denotes the Euclidian distance between the vectors ; ! v and ; ! n , V = q 2 x + 2 y + 2 z is the average volumetric spacing of the data. Our neighborhood structure (v) contains 26 neighbors in a 3-D metric space.

The iterative deterministic algorithm ICM (Iterated Conditional Modes) is implemented to compute the minimum energy at each voxel. The entropic deviation provides the estimation of the noise variance.

The local energy u i v in the slice i at the voxel v can be expressed as the following regularization function:

u i v = ( o v ; l v ) 2 + 2 : : 2 i : X n2 (v) V p (v n) (8)
The initial label configuration is the data itself. The normalization coefficient depends only on the range of the potential values, i.e. the size of the neighborhood. Typical values of are 1 and 0:2, respectively in a 6 and 26 neighbors system. Less than 10 iterations on the field are enough to respect the stopping criterion for convergence U= U<0:05%.

Results

Figure 1 compares results of 3-D filtering techniques when applied on a MR volume of the brain. The original slice reveals two distinguable but noisy areas of the brain: the white matter in light grey and the grey matter in dark grey. Whereas the median filter smoothes the boundaries, the anisotropic diffusion filter works well by keeping the initial edges [START_REF] Gerig | Nonlinear anisotropic filtering of MRI data[END_REF] 2 . Nonetheless, the MRF volume produces sharp edges and homogeneizes regions without user tuning. The distinction between anatomical areas is increased and proposes good preprocessing for any following segmentation process. The stability of the estimate provides also few variations in the quality of the filtering: the sharpness is constant along the slices independantly of the anatomy.

Discussion

In conclusion, the entropic deviation estimation has proved to be accurate, reliable and robust when applied on medical volumes independantly of the modality or the anatomy studied. When combined with MRF restoration, this method proposes an unsupervised way to remove the noise and sharpen the edges of anatomical structures. Contrary to empirical variance approaches, the quality of the results is constant along the slices.

The relevance of the entropic deviation is also under study for gradient and motion estimation. Other filters should gain from this noise estimation, like anisotropic diffusion or wavelet coefficients reduction.

The complete restoration pipeline is currently integrated in an interactive segmentation framework for maxillofacial surgery planning. First results are promising thanks to the robustness of the method.

Figure 1 .

 1 Figure 1. Restoration results. From top left to bottom right: the original slice; the median image; the anisotropic diffusion result; the proposed entropic-MRF filter

Table 1 . Cases studied for noise estimation.

 1

	1 2 3 4 5 Synthetic Sphere 256 256 128 CT Pelvis 512 512 46 CT Head 320 372 92 MR Knee 512 512 121 MR Head 256 256 136
	Case 1 56:8 19:1 21:5 9:56 16:3 sd Manual sf std std (:15) (:08) (:13) (:10) 2 76:6 10:8 23:5 6:54 9:2 (:30) (:31) (:23) (:27) 3 182 156 93 88 130 (:10) (:07) (:08) (:07) 4 21:5 12:4 11:2 8:31 9:6 (:13) (:14) (:23) (:21) 5 162 54:6 86:5 46:9 50:0 (:31) (:04) (:12) (:04)

Table 2 . Noise deviation estimation, where empirical deviation, std and entropic estima- tion, are applied on slice difference, sd and slice difference with its filtered value, sf. In bold, the closest value to manual measure- ment. Between braces, the "stability" crite- ria.

 2

The threshold on the gradient magnitude may be defined as the entropic estimation of the noise.


Théodore Totozafiny
email: t.totozafiny@estia.fr

Olivier Patrouix
email: o.patrouix@estia.fr

Franck Luthon
email: franck.luthon@univ-pau.fr

Jean-Marc Coutellier
email: jm.coutellier@magsys.net

Dynamic Background Segmentation for Remote Reference Image Updating within Motion Detection JPEG2000

We present in this paper a new system based on Motion JPEG2000 intended for road surveillance application. The system uses a reference image and consists in 4 processing steps, namely initialization phase where the first reference image is built, reference estimation, motion segmentation (foreground extraction, ROI mask), and JPEG2000 coding. A first order recursive filter is used to build a reference image that corresponds to the background image. The obtained background is sent to the decoder once for all. The reference image at the coder side is estimated according to a Gaussian mixture model. The remote reference image is updated when specific conditions are met. The updating remote reference is triggered according to the states of mobile objects in the scene (no, few or lot of mobiles). The motion detection given by classical background subtraction technique is performed in order to extract a binary mask. The motion mask gives the region of interest of the system. The JPEG2000 image coded with a ROI option is sent towards the decoder. The decoder receives, decodes the image and builds the implicit binary ROI mask. Then, the decoder builds the displayed image using the reference image, the current image and the mask.

I. INTRODUCTION AND PREVIOUS WORK

Our study addresses the problem of road surveillance for safety purposes. Recently [START_REF] Meessen | Smart Encoding for Wireless Video Surveillance[END_REF] developed an integrated system for smart encoding in video surveillance. A part of the system is based on Motion JPEG2000 (Part3 of the JPEG2000 standard) in which the mobile objects in the scene are extracted automatically with a background subtraction method. The reference image is estimated from each frame with a mixtures of Gaussian technique. The segmentation results are used as input region of interest masks for the JPEG2000 encoding of each frame. The frame is coded with ROI option of JPEG2000. The Motion JPEG2000 codes each frame independently and consists in a concatenation of JPEG2000 images. One of the requirements in JPEG2000 is the support of ROI coding, where ROI of the image can be coded with better quality than the background. The JPEG2000 image encoding standard defines two kinds of ROI coding techniques: the general scaling based method and the Maxshift method. Several methods [START_REF] Wang | Generalized Bitplane-by-Bitplane shift method for JPEG2000 Roi coding[END_REF] were proposed to improve the ROI coding. With the Maxshift method, no extra information about the shape of the mask is required for the decoder. The authors in [START_REF] Meessen | Smart Encoding for Wireless Video Surveillance[END_REF] combined the two above techniques to encode the ROI mask. Their system allows in real-time capturing an image, performing a motion detection, encoding with JPEG2000. The data is stored in a server. Then, the data will be transmitted toward a decoder in which a final image is reconstructed. The transmission is made with two layers. The first layer contains the data, the second contains the reference image. But the authors do not explain how to obtain initially the reference. Indeed, in the mixtures of Gaussian technique, some images are necessary to obtain a correct reference image. In order to obtain better results for reconstructed image at the decoder, an initialization phase should be done. [START_REF] Gutchess | A Background Model Initialization Algorithm for Video Surveillance[END_REF] gives a definition of the initialization of background and proposed a solution to build initially a reference image.

[4] proposed a scheme for region-based channel adaptive source coding scheme. The scheme proposed is compared with JPEG2000 encoding with ROI option. In video surveillance application with a transmission channel of very low bandwidth such as GSM network, the use of the reference image in Motion JPEG2000 encoding presents tree problems: obtaining a reference, updating and transmission towards the decoder of the reference regularly.

We have developed a system based on the Motion JPEG2000 standard. In an initial stage, the system builds a reference image and transmits it towards a decoder through the GSM network. This initialization is necessary in order to reconstruct the final image during the next transmission. When the initialization is achieved, the reference image is estimated according to a mixtures of Gaussian model. The motion detection is performed with classical background subtraction technique. The motion mask gives the region of interest for the system. The current image, containing only data linked with the mask, is coded using the ROI option of JPEG2000 encoding and transmitted towards the decoder. The remote reference image is updated when the conditions are required. The updating of remote reference is triggered according to the states of mobile objects in the scene (no, few or lot of objects). This paper is organized as follows: we describe our approach in section II and updating of the remote reference is given in section III. Finally experimental results and conclusions are given in sections IV and V.

II. DESCRIPTION OUR SYSTEM

In our approach, the embedded device is the one that is directly connected to the camera and the remote one is the ground station which receives the encoded images. The embedded device is mainly based on a camera, a CPU and a wireless transmission device.

A. Coder and decoder

Fig. 1 shows our coder system Motion Reference Image JPEG2000 (MRIJ2K). The system is based on four main processing blocks:

• Block A: initialization phase where the first reference image is built (i.e. the background image) and transmitted towards the decoder. This is done only one time at the starting of the system, • Block B: reference image estimation, • Block C: segmentation in order to extract a motion mask or reference mask, • Block D: JPEG2000 encoding with ROI option, Fig. 2 shows our MRIJ2K decoder system. The decoder receives, decodes the image and builds the implicit ROI mask. Then, it checks (indicated as FLAG in block diagram) whether the decoded image is a motion image or an updated remote reference image. If it is an updated reference, then the image is stored as the remote reference image otherwise the image is displayed.

B. Initialization phase

The initialisation phase consists of the construction of reference image. In a static camera context, several methods in the computer vision literature were proposed to obtain a reference image. In [START_REF] Stauffer | Adaptive background mixture models for real-time tracking[END_REF], [START_REF] Elgammal | Non-Parametric Model for Background Subtraction[END_REF], [START_REF] Gutchess | A Background Model Initialization Algorithm for Video Surveillance[END_REF], [START_REF] Wang | A Novel Probability Model for Background Maintenance and Subtraction[END_REF], the authors evaluate the advantages and drawbacks of some techniques. The common techniques used to extract the reference image are : median filter, frame differencing, Kalman filter, Gaussian mixtures, non-parametric method, wallflower, etc.

1) Constructing a reference image: Considering the constraints of the final custom application (lack of memory and embedded system working on a PC104 format board), we use the following filter transfer function form:

G(Z) = b 1 + aZ -1 (1)
where a and b are real. The recursive equation of the filter is given by:

I ref init (p, t + 1) = α p I(p, t) + (1 -α p )I ref init (p, t) (2)
where I ref init (p, t) and I(p, t) are the intensity values of pixel p in the reference image and in the current image at time t respectively. α p ∈ [0, 1] is the learning rate that gives the training speed, and p is the pixel location in the image (b = α p , a = α p -1). Of course when large object moves slowly in the scene, it will be included in the constructed reference image.

2) Learning rate filter: Since the reference image is the static part of the scene, the value of α p for all pixels belonging to a moving object must be 0. In order to determine α p we have to know all pixels that belong to the background. When p is a background pixel element, the α p value must be in the ]0, 1] interval, otherwise α p is 0. The binary stability charts of the three consecutive images I(t-2), I(t-1) and I(t) are used in order to know the object states in the scene. To compute the binary stability chart, we use the edge information and the noise of the acquisition camera is attenuated by an average filter (3 × 3). The Canny edge detector is used to calculate the gradient of each frame in the sequence.

d 1 = |I g (p, t -2) -I g (p, t -1)| (3) d 2 = |I g (p, t -1) -I g (p, t)| (4)
where I g is the spatial gradient of the processed image in the sequence.

A logical AND between d 1 and d 2 followed by an entropic thresholding technique indicates whether an object is moving or not: when the condition ((d 1 < λ) AN D (d 2 < λ)) is true, the object is a background element. λ can be obtained according to an entropy power threshold selection method [START_REF] Luthon | On the use of entropy power for threshold selection[END_REF]. The model [START_REF] Wang | Generalized Bitplane-by-Bitplane shift method for JPEG2000 Roi coding[END_REF] gives good results when mobile objects are of small sized, in the opposite case, the quality of the regions occupied by a large object is bad into the reference image obtained. It will be useful to improve these regions first.

C. Segmentation

In computer vision literature, several techniques were proposed to perform a segmentation of mobile object in the scene [START_REF] Stauffer | Adaptive background mixture models for real-time tracking[END_REF], [START_REF] Elgammal | Non-Parametric Model for Background Subtraction[END_REF], [START_REF] Wang | A Novel Probability Model for Background Maintenance and Subtraction[END_REF]. The motion detection is a binary labelling problem. It consists in attributing to each pixel p of an image I at time t one of the two following label values:

e p = 1 if p ∈ moving object 0 if p ∈ static background (5)
In order to carry out the binary labelling task, we use two observations. For each pixel p at time t, we calculate:

• the difference between the reference image and the current image:

o dr (p, t) = |I(p, t) -I ref (p, t)| (6)
• the difference between the two successive frames:

o dt = |I(p, t) -I(p, t -1)| (7)
To find the most probable label field with these two observations, we can use the conditional probability (Bayes theorem). [START_REF] Caplier | Robust fast extraction of video objects combining frame differences and adaptive reference image[END_REF] proposed a MAP criterion to obtain the most probable configuration. In our application, we used a logical AND in order to cope with the embedded target. We introduce a thresholding on both images o dr and o dt and then we compute the logical AND. So the pixel label is given by the following equation.

e p = ((o dr > T f ) AN D (o dt > T f )) (8)
The moving object threshold T f is determined according to an automatic entropy based threshold selection. The gradient mathematic morphology filter is used to improve the binary mask obtained. To determine a threshold T f , we used an entropy power for threshold selection [START_REF] Luthon | On the use of entropy power for threshold selection[END_REF]. T f is given by:

T f ≈ 2 H bit (9)
where H bit is the entropy of an information source (observations) that is classicaly defined as follows:

H = - N =255 i=0 pi log pi ( 10
)
where pi is the probability that the observation at any site takes the value i, and log denotes binary logarithm.

D. Estimation of the reference image

We used the same scheme as developed in [START_REF] Stauffer | Adaptive background mixture models for real-time tracking[END_REF] to estimate a reference image in embedded device. For each grabbed frame, the background image is updated. The method consists in the modelling of each pixel with K Gaussian distributions. P (x t ) represents the probability for a pixel to have the intensity x t at time t. This probability is estimated by:

P (x t ) = K j=1 ω j (2π) d 2 j 1 2 e -1 2 (xt-µj ) T -1 j (xt-µj ) (11)
where ω j is the weight of distribution j, µ j is the mean of distribution j and j is the covariance for the distribution j and d = 3. For computational reasons, the covariance matrix is assumed to be of the form j = σ 2 j I. The first B distributions are used as a model of the background of the scene where:

B = arg min b ( b j=1 ω j > T ) (12)
T is the fraction of the total weight given to the background model and 3 ≤ K ≤ 5. The parameters (ω j , µ j , σ 2 j ) of the matched and unmatched components are updated according a specific method.

E. JPEG2000 coding

JPEG2000 Part 1 has adopted a specific implementation of the scaling based ROI approach. This algorithm called Maxshift method consists in scaling down the background coefficients by 2 s . With the Maxshift method, no extra information about the shape of the ROI is required for the decoder [START_REF] Taubman | JPEG2000 Image Compression fundamentals standard and practice[END_REF]. This property is very significant in our scheme.

The ROI mask is automatically extracted from motion detection, and the image to be coded will only contain the useful data (no background data). In spatial domain, the pixels belonging to the ROI mask are set to the current image pixel value and for the others, the value is set to 128. This value is chosen because it is the dynamic mean value of an image coded with JPEG2000. This choice is done in order to control the transition effects between background and region of interest when the image is decoded.

Also, now the current image is compressed using a ROI option with a very low bit rate which will give a rate of 1 img/s with GSM channel (size of data 9600 bits around).Then, the compressed current image is transmitted towards the decoder.

F. Displayed image construction

The decoder receives an image. This image could be a reference image (i.e. background image) or a motion image (i.e. motion data). A flag is needed to identify those two cases and it is added in the file header during the image encoding. In the reference case, the received image is stored as background image on the decoder, otherwise the image is only displayed.

When a motion image or a reference image is decoded, the ROI binary mask is implicitly obtained. Due to the Maxshift technique property, the decision of whether or not a wavelet coefficient c belongs to the background (i.e. ROI) is stated using the following conditions:

c ∈ ROI if c ≥ dwt c ∈ ROI else ( 13
)
where dwt is a downshift parameter and dwt is set to 2 s . The image to be displayed is built using the reference image, the motion image and the mask reconstructed. In spatial domain, a simple pixel substitution is used to build the displayed image.

III. UPDATING OF THE REMOTE REFERENCE

Once the initialization phase is finished, the coder is able to send to the decoder the entire background image. Thus, the decoder can construct a final image : moving mask over the background. At the coder side, the reference image is updated for each frame according to a Gaussian Mixture Model (GMM). The updating of the reference image at the decoder should be done regularly.

The intuitive technique for updating the remote reference image (RRI) is based on the transmission of the RRI towards the decoder with a specific period. For example every 4 seconds or 10 images like in [START_REF] Meessen | Smart Encoding for Wireless Video Surveillance[END_REF]. The RRI must be coded with no ROI option and with a moderate bit rate compression. In our system, the transmission of the complete RRI slows down the transmission rate.

In our scheme, the image transmission rate should be done at 1 img/s. Since we only have at our disposal one single GSM transmission channel, the complete RRI update takes more than 4 seconds. This latency is not acceptable. To solve this problem and in order to keep the image transmission rate close to 1 img/s, we propose a new technique in order to update the RRI by pieces. In this scheme, the RRI will be coded like the motion image with a ROI mask. We have chosen a square pattern in order to build the ROI mask and the information in this area will be coded and sent to the decoder for the updating of the background image.

A. Regions to be updated

The background image is composed of N b blocks and each block will be used as an ROI mask. In order to improve the efficiency of our strategy, we define a parameter qc for each area. The qc represents a quality coefficient and qc non quality coefficient with respect to the initialization phase (background image creation) or the GMM background estimation. qc i and qc i (0 ≤ qc i , qc i ≤ 1) are used to compute a refreshing priority for the region i. The region with the highest priority will be updated first.

1) Local quality coefficient calculation: For a given frame, the local quality coefficient for a given block is based on the intersection of the motion ROI mask and the block area. The occupation percentage of mobile object is computed as follows. Lets Imask, Ir be the mask of motion object and current reference image respectively and i be the current index of treated block. We can write:

Ir = N b-1 i=0 Irb i (14)
Imask = N b-1 i=0 Imaskb i (15)
where Irb , Imaskb are sub-images of Ir and Imask respectively. For each block i, we compute the local quality coefficient lqc using the following equations :

lqc i = N onZero(Imaskb i ) N bsize 2 (16)
lqc i = 1 -lqc i (17)
where the N onZero is an operator that counts the number of pixels set to 1 in the current Imaskb and N bsize is the size of the square block.

2) Global quality coefficient calculation: At frame k, the quality coefficients are updated according to the local quality coefficients and the previous values of the quality coefficient at frame k -1. A threshold T qc is introduced to determine a changing of given region:

qc i (k) = M AX(qc i (k -1), lqc i (k)) if lqc i (k) ≥ T qc (1 -γ)qc i (k -1) else ( 18
)
where γ is the learning rate of the background estimation (ex. during the initialization phase, model (2) γ = α p ).

B. Remote reference updating decision

The region to be updated in remote reference must correspond to the block with the highest qc. We introduce two thresholding values T high and T down for observing the occupation of the mobile objects in the scene. The high threshold T high indicates that a lot of mobile objects is observed in the scene and T down indicates that few mobile objects are observed in the scene. The updating of the remote reference strategy is run according to the state of object mobile in the scene (much or less). Three cases are considered:

1) case 1: T state < T down , 2) case 2: T state > T high , 3) case 3: T down < T state < T high . where T state is the global percent state of mobile objects together and can be obtained by:

T state = N onZero(Imask) dim ( 19
)
where Imask is a mask of moving object and dim represents the size of the image.

C. Pulling of blocks

The blocks to be updated are chosen according to the three cases considered above.

In the first case, we can consider that no interesting object is present in the scene. We propose to use the block with the highest qc as the ROI mask for the JPEG2000 encoding. Then, the non-quality coefficient of this block is forced to zero (qc = 0, qc = 1) in order to flag it as "processed".

In the second case, too many moving objects are detected in the scene, then the updating of the remote reference should not be triggered. But if this configuration lasts a too long time (for example more than 2 hours), maybe the updating should be performed like in the previous case.

In the third case, there is no easy decision so we choose to force the updating every n frames in order to improve the global quality of the displayed image. In order to limit the background update, we define a threshold T qc and only the blocks with a qc greater than T qc will be updated. Then, the non-quality coefficient of this block is forced to zero (qc = 0, qc = 1) in order to flag it as "processed".

IV. EXPERIMENTAL RESULTS

For the JPEG2000 encoding, we used the Kakadu [START_REF]Kakadu SDK website[END_REF] codec SDK. Our goal is to obtain a 9600 bits image file in order to reach the image rate of 1 img/s with the GSM channel. At the ground station, the image reconstruction is done using the background image and the current JPEG2000 image. While decompressing the current image, we are able to retrieve the ROI mask which contains the spatial information (Maxshift method property). In spatial domain, the current finale image is computed by copying the background image and replacing the background pixel by the current pixel if it belongs to the ROI. For experimental tests, we are using 2 PCs. The transmission tests are done using a serial port, RS232 restricted to 9600 bps). The values of parameters are given in Tab.I. In order to evaluate the quality of final image obtained, we have used an objective classical criterion: PSNR. For each reconstructed image, the PSNR is calculated, yielding an average value of 32dB. Our experimental results are shown on Fig. 3, Fig. 4 and Fig. 5. The PSNR results are shown on Fig. 6. During of the initialization phase, the training frame number to extract the background image is set initially to 50. Our tests have been done with recorded video sequences of typical highway scenes 1 . In order to check the quality of the background reference image, a frame without any moving object is grabbed manually in the sequence and the PSNR is evaluated. The value is 30dB which leads to an average quality level. The quality is linked to the number of frames and also to the tuning of the first order parameters according to typical speed and size of the moving objects in the scene. As shown in left of Fig. 3, the quality of the regions occupied by large moving object is bad (a lapse of memory). Those regions having high values in the boxes of Fig. 3 bottom, should be updated first.

Then, the reference image is compressed to a JPEG2000 file at a moderate bitrate of 0.4 bpp, and sent once for all towards the decoder.

In Fig. 5, tree images are lost with the use of intuitive technique (here every 5 images) to update the remote reference image, whereas only one image is lost with our method. We have run tests with potential customers of our product and their global point of view is that our strategy gives a more fluid image flow than the intuitive method.

V. CONCLUSION

A video coding approach with Motion JPEG2000 using a reference image has been developed in the context of road surveillance. The complete scheme is implemented and it reaches the expected performances. We also showed how the image reference is built during initialization phase. The classical background subtraction technique is used to perform the segmentation of mobile objects. Instead of updating the remote reference with a specific period, we presented a new technique to update the remote background image by pieces. The updating of the remote reference is triggered when some specific conditions are met, depending on the amount of moving areas. The implementation of our system runs at 4 frames per second on a 1.6 GHz AMD processor for 320x240 color images. The transmission, using a restricted serial port, has been tested and the rate of 1 img/s has been reached. These results are close to the transmission which can be achieved with GSM channel.

Our future work will be focused on two axis. On one hand, we will build the first prototype using PC104 modules for Magys company2 and on the other hand, we will study the effect of LUX (Logarithmic hUe eXtension ) [START_REF] Luthon | with JPEG2000 for wireless videosurveillance[END_REF] on the rendering quality of the image in road surveillance context.

Notes that by using GPRS or future UMTS, one might expect a video rate of 5 img/s and 25 img/s respectively.

Fig. 1 .

 1 Fig. 1. Coder MRIJ2K.

Fig. 2 .

 2 Fig. 2. Decoder MRIJ2K.

Fig. 3 .

 3 Fig. 3. Reference image. From left to right: reference image built initially; reference estimated with a GMM (at frame 250); Bottom: represents the no quality coefficients values associed to each box.

Fig. 4 .

 4 Fig. 4. Final images ; from left to right: current image at coder, final image at decoder.

Fig. 5 .

 5 Fig. 5. RRI updating artefact. From top to bottom: image sequence of a scene; effect of simple updating of the RRI (every 5 images); our method for updating the RRI.

Fig. 6 .

 6 Fig. 6. PSNR comparision: Motion JPEG 2000 (MJP2) coding and our method.

Thanks to LAPS laboratory (http://www.u-bordeaux1.fr) for providing an image sequence; We have also tested our algorithm with the sequences available at http://i21www.ira.uka.de/image_sequences/

Magys is a SME involved in developing video devices mainly for video surveillance.