CR2P Denis Geraads
email: denis.geraads@mnhn.fr

Keywords: Mammalia, Neogene, Turkey, biochronology

The upper Miocene locality of Mahmutgazi in Western Turkey was excavated in the 70s by a German team, but most of its large mammals had never been studied. The collection housed in the Staatliches Museum für Naturkunde, Karlsruhe, contains, besides previously published groups, large samples of Giraffidae (Samotherium), Rhinocerotidae (including a nice complete skull of Ceratotherium neumayri), and Equidae, as well as some Chalicotheriidae (Ancylotherium) and Bovidae (Boselaphini), which are studied here. Three fossiliferous spots, Ma1, Ma2, and Ma3 have been recognized but I regard the first two, and probably also the third one, as contemporaneous. Although the fossil assemblage certainly does not reflect the full original taxonomic diversity, tentative biochronological conclusions can be proposed: the site is probably older than Şerefköy, Akkaşdağı, and the 'Dominant Fossil Assemblage' of Samos, but definitely younger than the early part of the Samos sequence.

Introduction

The locality of Mahmutgazi in Western Turkey was excavated by SICKENBERG and colleagues in the 70s but, although preliminary results were provided [START_REF] Sickenberg | Die Gliederung des höheren Jungtertiärs und Altquartärs in der Türkei nach Vertebraten und ihre Bedeutung für die internationale Neogen-Statigraphie[END_REF], most of the fauna of large mammals remains unpublished. Exceptions are the Carnivora [START_REF] Schmidt-Kittler | Carnivores from the Neogene of Asia Minor[END_REF], the Proboscidea [START_REF] Gaziry | Jungtertiäre Mastodonten aus Anatolien (Türkei)[END_REF], part of the Bovidae [START_REF] Köhler | Boviden der türkischen Miozäns[END_REF], and the Suidae [START_REF] Pickford | Hippopotamodon erymanthius (Suidae, Mammalia) from Mahmutgazi, Denizli-Çal basin, Turkey[END_REF]. During a visit to the Staatliches Museum für Naturkunde in Karlsruhe (SMNK) in December 2016, I studied the whole collection of large ungulates stored there, except the suids. Although it is not a large collection by comparison with several other sites of the Eastern Mediterranean, it is well-preserved and includes some excellent specimens.

The coordinates of the locality Mahmutgazi were provided by PICKFORD (2016) as 38°01'27.0"N, 29°24'34.9"E. These coordinates correspond well with field notes preserved in the SMNK and can be considered as satisfactory within ± 50 m. The locality includes two main fossil spots, Ma1 and Ma2, about 50 m apart in the same stratigraphic layer (W. MUNK, pers. comm.); an unpublished photo annotated by H. TOBIEN in the SMNK also shows Ma1, Ma2, and Ma3 at the same level. Indeed, although no data regarding the stratigraphy has been published, Tables 1 to 3 show that there is no consistent difference in size between the Giraffidae from Ma1 and Ma2, confirming that the assemblages from these spots can be regarded as a single fauna. BECKER-PLATEN et al. also regarded Mahmutgazi as a single site, but the example of the nearby site of Kemiklitepe, or of several Greek localities, where geographically very close fossiliferous spots yield distinct faunas, urges caution. Thus, I shall discuss separately the few specimens from Ma3. The collection remains largely un-accessioned, and many specimens are listed below by their locality number and the number of the plaster that contained them (Gips X).

Although no taphonomic analysis has been undertaken, the anatomical composition of the taphocenosis is far from faithfully reflecting that of living animals; this is especially true of the giraffes, whose limb bone collection consists almost entirely of adult or sub-adult specimens, probably from just a few individuals, whereas most of the teeth are deciduous; rhinoceros deciduous teeth are also more common than juvenile limb-bones. An attritional process would have accumulated a greater range of wear stages, with especially more senile individuals, and some sorting must have occurred at Mahmutgazi, either by carnivores or by running water; it may be that the fossil assemblage derives from one or a few catastrophic event(s) that preserved only an incomplete sample of the paleocommunity.

Systematic Palaeontology

Order Cetartiodactyla MONTGELARD, CATZEFLIS & DOUZERY, 1997Family Suidae GRAY, 1821'Hippopotamodon' erymanthius ROTH & WAGNER, 1854 As mentioned above, the suids have recently been revised by [START_REF] Pickford | Hippopotamodon erymanthius (Suidae, Mammalia) from Mahmutgazi, Denizli-Çal basin, Turkey[END_REF] who assigned them to this species, often included in Microstonyx.

Family Giraffidae GRAY, 1821

Genus Samotherium F. MAJOR, 1888 Samotherium major [START_REF] Bohlin | Die Familie Giraffidae[END_REF] Description and comparisons: This taxon is represented by a number of postcranials and several dental remains. The size range of the DP3s is perhaps greater than expected in a single species but all other remains, including the postcranials, are of homogeneous dimensions, and it is more parsimonious to assume that a single species is present. Ma2-Gips29 (Fig. 1B) is the only adult maxilla. The teeth are imperfectly preserved, but only slightly worn, preserving prominent styles and pillars. The premolars are rather large but P4 remains much narrower than M1; they bear strong spurs in the central valley, as may happen in Samotherium (BOHLIN 1926, pl.5, fig.10).

The best specimen is a partial mandible from Ma3 with the left i2-c and p2-m3; unfortunately part of the diastema is missing, so that its length is unknown (Fig. 1A). As in other Samotherium, the two lobes of the canine are much less clearly distinct than in Giraffa, and this tooth is not much larger than the incisors. The front teeth show no evidence of the interstitial wear that is observed in some other late Miocene large giraffids of the Aegean region, such as Helladotherium from Hadjidimovo, Bulgaria [START_REF] Geraads | Giraffidae (Artiodactyla, Mammalia) from the late Miocene of Kalimantsi and Hadjidimovo, South-Western Bulgaria[END_REF], Samotherium major from Vathylakkos, Greece [START_REF] Geraads | Les Giraffidés du Miocène supérieur de la région de Thessalonique (Grèce)[END_REF]MERCERON et al. in press), or S. neumayri from Maragha, Iran (RODLER & WEITHOFER 1890, p1.4, figs.2-3). This peculiar wear, which is certainly linked with their use as a comb when browsing, is also variably present in Giraffa camelopardalis. The premolars are strongly molarized: p2 has a long entocristid; p3 is similar to p4 in having a long metaconid connected to the paraconid, but not to the transverse cristid, and blocking the entrance to the mesosinusid; the ectoflexid of p4 is deep, as in other Giraffidae, but the heavy wear has connected the hypoconid to the rest of the tooth. These premolars are also quite short, with a Pm/M index of 58.6 (Table 1). These short, molarized premolars suffice to identify the Mahmutgazi large giraffid as Samotherium. In the large Helladotherium, a genus common in Aegean sites of similar age, the premolars are significantly larger, and p3 always retains a primitive morphology. Most representatives of the species S. boissieri from Samos also retain this primitive morphology (BOHLIN 1926, figs.141-146;KOSTOPOULOS 2009a, Fig. 8), but S. neumayri from Maragha and S. major from Vathylakkos and Samos [START_REF] Bohlin | Die Familie Giraffidae[END_REF][START_REF] Geraads | Les Giraffidés du Miocène supérieur de la région de Thessalonique (Grèce)[END_REF]KOSTOPOULOS 2009a) have a fully molarized p3.

The shortening of the premolars also fits better S. major (Pm/M index = 53.8-66.2, mean = 59.5, N = 6) than S. boissieri (Pm/M index = 60.7-65.6, mean = 63, N = 5) (measurements from [START_REF] Geraads | Les Giraffidés du Miocène supérieur de la région de Thessalonique (Grèce)[END_REF]KOSTOPOULOS 2009a, tab. 8). However, the absolute size of the teeth (Table 1) is lower than in S. major, and closer to S. boissieri.

Upper DP3s (Fig. 1K) have a relatively short first lobe; lower dp3s (Fig. 1L) also have large metaconids completely closing the mesosinusid, and the first lobe of dp4 is fully formed, with no labial arm; these are all derived features, usually found in Samotherium.

The postcranial sample is large. There are no less than nine distal humeri, five of them from the right side, providing a Minimum Number of Individuals of five for the whole collection. It is likely that all these postcranials are indeed from a limited number of individuals, but they cannot be matched, because a single plaster jacket may have contained bones of more than one individuals (Ma2-Gips20 has two right metacarpals).

The dimensions of the distal humeri (Table 2) are homogeneous, but there is no complete bone (Fig. 1E). The available sample of measurements from nearby sites is small; these humeri are clearly larger than those of S. boissieri from Samos (dist.TD = 98-106 according to KOSTOPOULOS [2009a]), but smaller than those of its likely descendant S. major from the same area (dist.TD = 132-148). The Mahmutgazi humeri also compare well with the forms from the two levels of Kemiklitepe (KTD, lower, and KTA, upper: GERAADS 1994), and with S. neumayri from Maragha.

There are six more or less complete metapodials (Table 3); the measurements of the metacarpals (Fig. 1F-H) clearly plot with S. major (Fig. 2).

There are three tibiae (Fig. 1C-D), whose measurements (Table 4) are clearly more similar to those of S. major than to those of S. boissieri.

There are only four astragali, two of them well preserved (Fig. 1J). One is similar in size to those from KTD, and the other to those from KTA (GERAADS 1994), but they also fit S. neumayri from Maragha (pers. data). It is remarkable that, while the size distinction between the two levels of the nearby site of Kemiklitepe is clear, the size of the Mahmutgazi material overlaps both of them (Fig. 3), but the difference in size between the two well-preserved specimens does not exceed what is found in a single population.

A number of other bones (carpals, tarsals, calcanei, phalanges) obviously belong to the same species; the first phalanges are short and stout, as in Samotherium [START_REF] Geraads | Les Giraffidés du Miocène supérieur de la région de Thessalonique (Grèce)[END_REF]; the cubo-navicular has no posterior articular facet for the metatarsal (Fig. 1I), as in other Samotherium, but in contrast to Helladotherium, in which it is present [START_REF] Geraads | Les Giraffidés du Miocène supérieur de la région de Thessalonique (Grèce)[END_REF]KOSTOPOULOS 2009a).

Discussion: All postcranials can easily be assigned to Samotherium major, and the derived lower premolars fit this identification, but the lower teeth from Ma3 are smaller than in this species. If indeed Ma3 is contemporaneous with Ma1 and Ma2, these teeth may be from a small individual of S. major. In terms of biochronology, S. major clearly allies Ma1 and Ma2 with the upper part of the Samos sequence (KOSTOPOULOS 2009a, fig.21).

Family Bovidae GRAY, 1821

The SMNK collection contains only a few tooth rows, but no postcranials or horn cores; this is a very unlikely anatomical composition for a late Miocene Aegean site, and some of them must be stored elsewhere, but I have been unable to locate them. [START_REF] Köhler | Boviden der türkischen Miozäns[END_REF] described some more cranial remains, now stored in the Geozentrum Hannover. I have not seen them, and her identifications were based upon very limited material. She listed:

-Tragoportax gaudryi, based upon a maxilla (M1-M3 = 53.8); its measurements indeed fit what SPASSOV & GERAADS (2004) called Miotragocerus(Pikermicerus) gaudryi. In addition, a mandible in the SMNK (Fig. 1N) also fits this species in size and very long premolars (p2-p4 = 47.5; m1-m3 = 62), and this identification can be accepted.

-Tragoportax amalthea, based upon a good sample of tooth rows; the genus identification is also acceptable, according to her fig. 11; a maxilla in the SMNK (Fig. 1M) probably also belongs here.

-Pseudotragus parvidens; the horn-core cross-section (KÖHLER 1987, fig. 58) does indeed resemble that of this species, now placed in Sporadotragus (see, e.g., GERAADS et al. 2006a), but the identification remains ambiguous.

-'Plesiaddax' cf. inundatus, based upon mandibular pieces, including an m3 reaching a length of 40 mm; in the Aegean region, such a large size is reached only by this species, defined at Garkın [START_REF] Erdbrink | Fossil ovibovines from Garkın near Afyon, Turkey[END_REF], so that this identification is probably correct.

-Gazella sp., 'Form IVa'; the horn-core cross-section and measurements would fit G. capricornis, and I tentatively assign them to this species.

-Gazella sp., 'Form V'; this species is mostly illustrated by a hornless skull, whose dental dimensions are also very similar to those of G. capricornis, a species whose females are said by KOSTOPOULOS (2009b) to be hornless, and it is more parsimonious to also identify 'Form V' with G. capricornis.

-Gazella sp., represented by a single horn-core said to be very large, but whose measurements (c. 28 × 25.4 according to KÖHLER 1987, fig. 84) would also fit G. capricornis (see, e.g., KOSTOPOULOS 2009b, fig. 6).

-Palaeoreas cf. lindermayeri, represented by a single maxilla with P3-M2; such a piece cannot be identified, even at genus level.

-Oioceros wegneri, represented by the tip if a horn-core; members of the Oioceros group are characterized by their horn cores with homonymous torsion, but the tip of a horn-core cannot be sided, so the identification remains uncertain.

-Palaeoryx pallasi, represented by a mandible, a maxilla, and two horn-core fragments; measurements and figures of the specimens (KÖHLER 1987, fig. 13 and103) are satisfactory matches for the species, and this identification can be accepted.

Order Perissodactyla OWEN, 1848 Family Rhinocerotidae GRAY, 1821 Genus Ceratotherium GRAY, 1867

Ceratotherium neumayri (OSBORN, 1900) This species is represented by some postcranials and tooth-rows, and a cranium Ma2-Gips15, virtually complete and undistorted, which is, together with the specimen from Akkaşdağı (ANTOINE & SARAÇ 2005) one of the best known specimens of this species, and the most valuable fossil of the SMNK collection from Mahmutgazi (Fig. 4A; comparative measurements: Table 5). The species is relatively well-known, and a detailed description of the cranium would be useless; I shall focus on the most significant points.

The cranium carried two horns, a large nasal ones, and a smaller, but distinct frontal one. Behind it, the dorsal profile is straight as far as the level of the middle of the temporal fossa, where it curves upwards towards the nuchal crest. Although well indicated, this concavity, in this undistorted skull, is much shallower than in the skull MTLA-5 from Samos (GIAOURTSAKIS 2009, pl.1, fig.1). This strongly suggests that distortion in this Samos skull was stronger than assumed by [START_REF] Giaourtsakis | The late Miocene Mammal faunas of the Mytilinii Basin, Samos Island, Greece: new collection. 9. Rhinocerotidae[END_REF], as also indicated by the relatively short distance between the occipital condyles and M3, and the upward tilt of the zygomatic arch. In The rarely preserved premaxillae are complete; as usual in the Dicerotini, they bear no incisor, although small cavities may be the alveoli of vestigial ones. The bottom of the nasal notch is above the limit P2/P3, the infra-orbital foramen is above P4, the front of the orbit above the middle of M2; these positions are similar to those observed in other C. neumayri (GERAADS 2005, fig. 6). As in modern African rhinos, the lachrymal tubercle and supra-orbital process tend to fuse into an orbital rim, there is no post-orbital process, and the ventral orbital border is rounded.

The DP1 is shed, and the other cheek-teeth are well-worn. They are virtually identical with those of D. bicornis, a morphology that, in my opinion, confirms the hypothesis that C. neumayri is an acceptable ancestor for both D. bicornis and C. simum, although the latter has a much derived tooth morphology.

Ma2-Gips16 (Fig. 4C) is a maxilla with nasal notch and infra-orbital foramen slightly more posterior than in the skull, but again within the range of C. neumayri. The teeth are slightly less worn, but display the same characters, being very similar to those of D. bicornis.

Their measurements are also similar to those of Ma2-Gips15 (occlusal length P2-P4 = 121; occlusal length M1-M3 = 145). There is no tooth anterior to P2, meaning that DP1 was shed, and not replaced; it seems that this is the rule in D. bicornis, although a true P1 may sometimes show up (MNHN-1961-195). The premolars have small cristae, as sometimes occurs in D. bicornis.

Ma3-Gips1 is a nice deciduous upper tooth-series, very similar to those figured by [START_REF] Geraads | Révision des Rhinocerotidae (Mammalia) du Turolien de Pikermi. Comparaison avec les formes voisines[END_REF][START_REF] Geraads | Les gisements de Mammifères du Miocène supérieur de Kemiklitepe, Turquie: 4: Rhinocerotidae[END_REF] and [START_REF] Giaourtsakis | The late Miocene Mammal faunas of the Mytilinii Basin, Samos Island, Greece: new collection. 9. Rhinocerotidae[END_REF]. Its size (L DP1-DP4 = 151 mm) is within the range of Pikermi [START_REF] Geraads | Révision des Rhinocerotidae (Mammalia) du Turolien de Pikermi. Comparaison avec les formes voisines[END_REF], main Samos assemblage (GIAOURTSAKIS 2009), and slightly above that of the two specimens from Kemiklitepe A (GERAADS 1994). There are several lower deciduous series. They are comparable in size to those from MTLA (GIAOURTSAKIS 2009) but the talonid basin of dp2 is never fully closed, and the anterior lobe of dp3 is at most incipient; it may be that the intra-specific variability of these teeth is greater than assumed by [START_REF] Giaourtsakis | The late Miocene Mammal faunas of the Mytilinii Basin, Samos Island, Greece: new collection. 9. Rhinocerotidae[END_REF].

It is more parsimonious to assume that post-cranials also belong to C. neumayri; indeed, they display some of the features regarded as characteristic of this taxon by [START_REF] Geraads | Révision des Rhinocerotidae (Mammalia) du Turolien de Pikermi. Comparaison avec les formes voisines[END_REF] and [START_REF] Giaourtsakis | The late Miocene Mammal faunas of the Mytilinii Basin, Samos Island, Greece: new collection. 9. Rhinocerotidae[END_REF].

There are several sets of associated metacarpals (Fig. 4D). The McII has no palmar ridge and a long articular surface for the magnum, and the McIII has a deep lateral concavity to accommodate a strong tuberosity of the Mc IV (see also GIAOURTSAKIS 2009, pl. 5, figs.4-6).

The three McIIIs are homogeneous in size (Table 6), and are slightly shorter than a specimen from Lubéron (MNHN), but slightly longer than a McIII from Garkın in München (1968-VI-359), and than KTD-6 from the lower level of Kemiklitepe [START_REF] Geraads | Les gisements de Mammifères du Miocène supérieur de Kemiklitepe, Turquie: 4: Rhinocerotidae[END_REF]. Although not too much weight should be placed upon these measurements, the metric evidence is consistent with a slight increase in metapodial length from KTD and Garkın to Mahmutgazi and to Lubéron.

Ma2-Gips15 (Fig. 4B) is an astragalus showing the typical features of the species (GIAOURTSAKIS 2009): oval sustentacular facet separated from the distal calcanear facet, lateral fibular facet tapering proximally, not contacting the sustentacular facet. Ma2-Gips28 is an incomplete astragalus that differs in having a round sustentacular facet, and in being lower, but it may have suffered some distortion; I do not regard this single specimen as positive evidence of another species.

Thus, Ceratotherium neumayri is the only species definitely represented in the SMNK collection. HEISSIG (1975) mentioned Chilotherium schlosseri, but I have not seen a single specimen that could be assigned to Chilotherium.

Family Chalicotheriidae GILL, 1872

Genus Ancylotherium GAUDRY, 1863

Ancylotherium pentelicum (GAUDRY & LARTET, 1856) This strange and highly characteristic animal is represented by the poorly preserved proximal portions of left associated metacarpals, Ma2-Gips21 (Fig. 5H). Their length cannot be estimated, but what remains shows that they were stoutly built, and certainly shorter and less slender than those from Hadjidimovo (GERAADS et al. 2006b), more like those from Maragha (MECQUENEM 1924-25;[START_REF] Roussiakis | Ancylotherium pentelicum (Gaudry & Lartet, 1856) (Perissodactyla, Mammalia) from the classic locality of Pikermi (Attica, Grece), stored in the palaeontological and geological museum of Athens[END_REF]. The proximal epiphysis of the Mc III is broader than deep (in dorsopalmar direction), unlike all other specimens (Table 7), and the facet for the magnum is relatively broad. The meaning of this difference, and the post-cranial evolution of the species, which is present in the area since the Vallesian (KOUFOS 2012), remain unknown. Ancylotherium was previously assumed to be a leaf browser, but both GERAADS et al.

(2006b) and SCHULZ ET AL. ( 2007) came to the conclusion that it fed on bark or stems of large plants. This specialization may explain why it is always rare at any site, although present in many of them (map in GERAADS et al. 2006b).

Family Equidae GRAY, 1821 I shall only briefly report on this group. STAESCHE & SONDAAR (1979) mentioned that they included the Karlsruhe material in their study of the Turkish hipparions but, although they provide the measurements of a number of teeth, they do not mention postcranials from Mahmutgazi, and seemingly included all hipparions from this site in a single species, Hipparion matthewi.

The most significant fossils in the SMNK are eight metapodials and a number of upper and lower tooth-rows, but no remain includes the preorbital part of the skull, consistently used in hipparion systematics. The main measurements of six relatively complete metatarsals are plotted in Fig. 6; they clearly indicate the presence of two species, each with remarkably little intra-specific variation. The lengths of the bones of the two species are not very different, but the shorter ones are much more slender. The measurements of two metacarpals display the same pattern: one is shorter and much more slender than the other (I follow here the tradition in using Simpson diagrams, although in fact they show little more than gross differences in size; obviously, multivariate analyses would be better adapted but, beyond a brief attempt by [START_REF] Koufos | Equidae (Mammalia, Perissodactyla) from the late Miocene of Akkaşdağı, Turkey[END_REF], they have been neglected).

Differences in tooth-row lengths (Table 8) are also quite clear, and the distribution of individual tooth lengths is clearly bimodal (STAESCHE & SONDAAR 1979, fig. 26), confirming the existence of two species.

Cremohipparion cf. matthewi (ABEL, 1926)

There are two complete upper tooth-rows (Fig. 5D), one with P3-M3, probably also a mandible with the front teeth and the premolars, and some more or less lower complete toothrows (Fig. 5E-G). The length of the series P2-P4 ranges from 64 to 68 mm, that of M1-M3 from 53 to 57 mm. The protocone is large and rounded on the upper premolars, and tends to connect anteriorly to the protoloph; it is smaller and more transversally compressed on the molars. The pli caballin is small, and sometimes multiple on the premolars. On a mandible, the diastema is rather short and the third incisor looks unreduced, but the specimen is not fully adult.

The metapodials are shorter and more slender than all other Turkish specimens measured by [START_REF] Staesche | Hipparion aus dem Vallesium und Turolium (Jungtertiär) der Türkei[END_REF] from Kınık, Garkın, or Kayadibi, but are comparable with those of C. matthewi from Kemiklitepe A-B or Samos Q5 (KOUFOS & VLACHOU 2005, figs.43-44).

Although no substantial part of the skull is preserved, the size of the tooth-rows, and those of the short, slender metapodials unambiguously indicate that this hipparion belongs to a group of small-size forms, of which Cremohipparion matthewi is the most common representative, and we tentatively assign the small hipparion from Mahmutgazi to this species.

Hipparionini gen. et sp. indet.

A larger form seems to be slightly more common but, except for its larger size and more robust metapodials, it does not display any significant difference, except that enamel plication is perhaps stronger (Fig. 5A-C). The premolars are relatively slightly shorter than in the C. cf.

matthewi. The size of the metapodials would fit several Late Miocene Aegean species, and I prefer not to suggest any identification.

Order Proboscidea ILLIGER, 1811

Genus Choerolophodon SCHLESINGER, 1917

Choerolophodon pentelici (GAUDRY & LARTET, 1856)

The SMNK collection includes only a score of postcranials. Because late Miocene Proboscidean postcranials have been neglected, these bones neither confirm nor refute GAZIRY's (1976) identification of Choerolophodon pentelici, but since this identification was based upon a tooth fragment, it can be accepted.

Biochronology

Both the scarce field data and the study of the collection leave no doubt as to the contemporaneity of the faunas from Ma1 and Ma2, which can thus be regarded as a single one.

The case of Ma3 is less clear, because the presence of Protictitherium crassum (see [START_REF] Schmidt-Kittler | Carnivores from the Neogene of Asia Minor[END_REF]) and of a Samotherium mandible of relatively small size might suggest an earlier age. However, as noted above, a photo preserved in the SMNK shows Ma3 almost exactly at the same stratigraphic level as Ma1 and Ma2, and two deciduous tooth-series of Ceratotherium neumayri (see above) fit within the range of typical Turolian forms. The issue cannot be definitely be settled but, on balance, Ma3 looks contemporaneous with the other spots, although it is safer to discuss it separately.

An updated faunal list based upon the present study but also taking into account previous identifications [START_REF] Sickenberg | Die Gliederung des höheren Jungtertiärs und Altquartärs in der Türkei nach Vertebraten und ihre Bedeutung für die internationale Neogen-Statigraphie[END_REF][START_REF] Gaziry | Jungtertiäre Mastodonten aus Anatolien (Türkei)[END_REF][START_REF] Schmidt-Kittler | Carnivores from the Neogene of Asia Minor[END_REF][START_REF] Köhler | Boviden der türkischen Miozäns[END_REF]PICKFORD 1976) can now be provided (Table 9). I have not seen any material of Orycteropus, but the genus is so characteristic that the identification [START_REF] Becker-Platen | Vertebraten-Lokalfaunen der Türkei und ihre Alterstellung[END_REF]) can be accepted. I list the carnivores as they were published, although the status of the primitive hyenids is unclear [START_REF] Werdelin | The Hyaenidae: taxonomy, systematics and evolution[END_REF], and the distinction between Machairodus aphanistus and later M. giganteus far from clear-cut. This faunal list is rather short compared to many upper Miocene sites of the area, where Giraffidae and Rhinocerotidae are usually represented by more than one species. This is not merely due to the limited scale of the excavations because the Giraffidae and Rhinocerotidae, with about 80 and 40 specimens, respectively, make up a significant collection. It is likely that the accumulation process, whatever it was, very incompletely sampled the large mammal biocenosis; paleoecological inferences would not be strongly supported, but some biochronological comparisons with other Turkish sites can be proposed.

The Sinap Formation near Ankara includes some early Turolian faunas (in localities 26, 27, and 33 of Kavakdere) consisting mostly of fragmentary, hardly identifiable material. The Putting the various localities of Samos and Western Turkey in chronological order is difficult because, as shown by Table 9, most assemblages look like mere sub-samples of the DMAS, which is distinctly richer than the others. Inter-locality differences mostly concern taxa that are rare or hard to identify. Another problem is that the only fully studied group in all sites are the ruminants, which are poorly represented at Mahmutgazi.

Mahmutgazi is certainly younger that the lower part of the Samos sequence (PMAS), which has a distinctly more primitive species of Samotherium, but it can be observed that no other taxon of the PMAS confirms this early age. However, Protictitherium crassum, unknown afterwards, is a resemblance between both sites, if indeed Ma3 is contemporaneous with Ma1 and Ma2.

At Garkın, ERDBRINK (1978) defined 'Plesiaddax' inundatus that is probably also present at Mahmutgazi, but other listed taxa (BECKER-PLATEN et al. 1976) do not allow conclusions, although the site is traditionally regarded as transitional between zones MN 11 and 12 (meaning that these zones are not properly defined).

The lower level of Kemiklitepe [START_REF] Bonis | Les gisements de Mammifères du Miocène supérieur de Kemiklitepe, Turquie: 11. Biochronologie, paléoécologie, et relations biogéographiques[END_REF]) is not very rich, but the metacarpals and astragali of Samotherium are smaller than those of Mahmutgazi, and there is an antelope that resembles Majoreas elegans from the Sinap, so that KTD can confidently be regarded as older than Mahmutgazi. Akkaşdağı is a rich and fully studied locality in central Turkey [START_REF] Sen | Geology, mammals and environments at Akkaşdağı, late Miocene of central Anatolia[END_REF]. Some of its taxa are unknown in Western Turkey, although paradoxically present in continental Europe (Tetralophodon, Chalicotheriinae, Thalassictis spelaea), but few are biochronologically significant; [START_REF] Kostopoulos | Late Miocene bovids from Şerefköy-2 (SW Turkey) and their position within the sub-Paratethyan biogeographic province[END_REF] identified Skoufotragus schlosseri, a similarity with Şerefköy-2, and an age of c. 7.1 Ma has been proposed on the basis of radiometric dating [START_REF] Karadenizli | Mammal bearing late Miocene tuffs of the Akkaşdağı region: distribution, age, petrographical and geochemical characteristics[END_REF]. The hyenid assemblage and the presence of Chilotherium suggest that Akkaşdağı is younger than Mahmutgazi but evidence is weak.

Regarding previously published groups, the characters of the large carnivores studied by SCHMIDT-KITTLER (1979), Machairodus and Adcrocuta, would fit any age in the Late Miocene, except perhaps its very beginning, but the large size of the Machairodus P3 is perhaps more indicative of a relatively ancient age. PICKFORD (2016) concluded that the size of the Mahmutgazi 'Hippopotamodon' plots in the lower half of the Pikermi and Samos ranges.

Conclusions

In summary, comparisons with other sites broadly agree with previous estimates made by [START_REF] Heissig | Rhinocerotidae aus dem Jungtertiärs Anatoliens[END_REF], [START_REF] Becker-Platen | Vertebraten-Lokalfaunen der Türkei und ihre Alterstellung[END_REF][START_REF] Köhler | Boviden der türkischen Miozäns[END_REF], and I regard Mahmutgazi as earlier than Şerefköy-2, Akkaşdağı, or the DMAS, but younger than Kemiklitepe-D. I observe that assigning it to a 'MN zone' is not straightforward, thus demonstrating that these zones cannot be defined on the basis of large mammals in Turkey, and therefore have no meaning there. In fact, assignment of Aegean large mammal faunas to 'MN zones' usually rests upon the estimation of their age by comparison with other sites, or worse, with the help of some physical method. Because their existence is not implied by the facts (e.g., at Samos, most of the faunal changes occur within 'MN 12' [KOUFOS et al. 2009, fig. 2]), and because they have no practical use there, dropping these 'zones' from use would not hinder progress in Aegean Miocene biochronology.

Captions to figures

  Abbreviations: APD = anteroposterior diameter; dist. = distal; min. = minimum; max. = maximum; prox. = proximal; TD = transverse diameter; H = height; W = width. EUMI = Aegean University Museum, Izmir; MNHN = Muséum National d'Histoire Naturelle, Paris; NHMA = Natural History Museum, Asenovgrad, Bulgaria; NHMUK = Natural History Museum, London. Measurements are in mm. Upper teeth are in upper case, lower teeth in lower case.

  this relatively straight dorsal profile, Ma2-Gips15 resembles modern C. simum, the African 'white' rhinoceros. The nuchal crest has a large flange-like caudal expansion, with a central notch; it largely overhangs the condyles, and is almost as broad as the ventral part of the occipital; the occipital surface is slightly inclined dorsocaudally, thus intermediate between the average conditions of C. simum and Diceros bicornis, the 'black' rhinoceros. The basioccipital and basisphenoid form a marked angle, a clear difference with C. simum, in which they are almost in line, in correlation with the much lengthened skull of the end members of the Ceratotherium lineage.

  large giraffid, however, is certainly distinct (GERAADS & GÜLEÇ 2000), Prostrepsiceros zitteli is definitely present (GERAADS & GÜLEÇ 1999), together with antelopes of the Pachytragus-Protoryx group (GENTRY 2003) and there are certainly at least two species of rhinos (FORTELIUS et al. 2003). The late Turolian locality of Çoban Pinar (Loc. 42) has the same antelopes of the Pachytragus-Protoryx group, plus Nisidorcas (GERAADS & GÜLEÇ 1999), and a small suid (VAN DER MADE 2003). Çorakyerler near Ankara (GERAADS 2013) much differs from Mahmutgazi, and is probably earlier. The dominant rhinoceros is a hornless form with large incisors, there is no form close to C. matthewi, no Samotherium, and the dominant bovids are spiral-horned forms; the only significant resemblance is the likely presence of 'Plesiaddax' inundatus. Comparison with sites in Central and Western Turkey is more relevant. Unfortunately, the faunas from the other sites excavated under the leadership of O. SICKENBERG have also been very incompletely studied, so that the validity of the biochronology established by SICKENBERG et al. (1976) remains questionable. STAESCHE & SONDAAR (1979) placed Mahmutgazi at the end of the Garkın faunal zone, or at the beginning of the Kınık one. It would thus be intermediate in age between Samos Q1/Q4 and Samos Q5, corresponding to what KOUFOS et al. (2009) called 'Dominant Mammal Assemblage' of Samos (DMAS).

  The upper levels of Kemiklitepe and the Intermediary Mammal Assemblage of Samos (IMAS) share a number of taxa with Mahmutgazi but none is very significant in terms of biochronology; the limb bones of Samotherium major are similar-sized, but a Chilotherium s.str. is present; I regard KTA-B and the IMAS as roughly contemporaneous with Mahmutgazi, or perhaps slightly younger.The fauna from Şerefköy-2 has only been partly studied, but the Machairodus was listed as M. giganteus[START_REF] Kaya | Şerefköy-2, a new Late Miocene mammal locality from the Yatağan Formation, Muğla, SW Turkey[END_REF], and some taxa also perhaps indicate a younger age (Parataxidea, Pliohyrax graecus, Urmiatherium rugosifrons); the only possible significant resemblance with Mahmutgazi is a very large bovid, assigned to Palaeoryx pallasi by KOSTOPOULOS & KARAKÜTÜK (2015, fig.4A-B), but reminiscent of the Mahmutgazi 'Plesiaddax', although its premolars are longer. I regard Mahmutgazi as older than Şerefköy-2, whose age was estimated at c. 6.8-7 Ma by[START_REF] Kostopoulos | Late Miocene bovids from Şerefköy-2 (SW Turkey) and their position within the sub-Paratethyan biogeographic province[END_REF].

Fig. 1 .

 1 Fig. 1. Artiodactyla from Mahmutgazi. A-L, Samotherium major. A: left mandible Ma3 (A1: occlusal view; A2: lateral view of the rostral part). B, maxilla with P3-M2 Ma2-Gips29. C-D: tibiae in anterior view (C: MA2-Gips16; D: Ma2-Gips11). E: distal humerus n°54, anterior view. F-H: metacarpals in anterior view (F: Ma2-Gips17; G and H: Ma2-Gips20). I: cubonavicular Ma1-Gips10, distal view, to show the absence of plantar metatarsal facet. J: astragalus Ma1-Gips10, anterior view. K: upper deciduous tooth-row Ma1-Gips3. L: lower dp3-dp4 Ma1-Gips11. M: Tragoportax sp., right upper tooth row. N: Miotragocerus (Pikermicerus) gaudryi, left lower tooth-row. Scale bar = 5 cm for Figs. M-N, 10 cm for Figs. A, B, and I-L, 25 cm for Figs C-H.

Fig. 2 .

 2 Fig. 2. Plot of length vs. distal width of late Miocene giraffid metacarpals (KTA = Kemiklitepe, upper level; Slq = lower Axios valley, ARAMBOURG'S collection). Data from KOSTOPOULOS (2009a), and original measurements.

Fig. 3 .

 3 Fig. 3. Plot of medial length vs. distal width of late Miocene giraffid astragali. Note that 'Maragha' includes both S. neumayri and Palaeotragus coelophrys. Data from KOSTOPOULOS (2009a), and original measurements.

Fig. 4 .

 4 Fig. 4. Ceratotherium neumayri from Mahmutgazi. A: skull Ma2-Gips15 (A1: lateral view; A2: ventral view; A3: occipital view). B: astragalus Ma2-Gips15 (B1: plantar view; B2: lateral view; B3: anterior view). C: upper tooth row Ma2-Gips16, occlusal view. D: set of left metacarpals Ma1-94. Scale bar = 20 cm for Figs. A and D, 10 cm for Fig. B and C.

Fig. 5 .

 5 Fig. 5. Perissodactyla from Mahmutgazi. A-C: Hipparionini gen. et sp. indet. (A: upper P2-M3 Ma2-Gips24; B: upper P2-M3 Ma2-Gips3; C: lower p2-m3 'M' [reversed from the right side]). D-G: Cremohipparion cf. matthewi (D: upper P2-M3 Ma1-Gips 17 [reversed from the left side]; E: lower p3-m3 'K'; F: lower p2-m2 'H'; G: lower p3-m3 'J' [reversed from the right side]). H: Ancylotherium pentelicum, front view of proximal metacarpals Ma2-Gips21. Scale bar = 5 cm for Figs. A-G, 10 cm for Fig. H.

Acknowledgments

I am grateful to Wolfgang MUNK for welcoming me at the SMNK and for much help there, to Martin PICKFORD for calling my attention to the Mahmutgazi collection, and to the reviewer, Dimitris S. KOSTOPOULOS, for his enlightened comments.

= proximal breadth; 10 = distal supra-articular width; 11 = distal articular width; 12 = maximum distal antero-posterior diameter.


M Bergounioux
email: maitine.bergounioux@univ-orleans.fr

E Le Pennec
email: erwan.le-pennec@polytechnique.edu

E Trélat

Fractional-order variational numerical methods for tomographic reconstruction of binary images

Keywords: tomographic reconstruction, Radon operator, fractional Laplacian, minimization, Fourier transform, wavelets

 we have provided a refined functional analysis of the Radon operator restricted to axisymmetric functions and proved that it enjoys strong regularity properties in fractional order Hilbert spaces. Based on that theoretical study, we provide here the details of the numerical solving of a minimization problem settled in suitable fractional order Hilbert spaces, using a numerical approximation of the fractional Laplacian and some adapted Newton-like methods. The resulting procedure happens to be very efficient, both for the execution time and for the quality of the reconstruction. We compare this approach with two other numerical approaches: the first one uses the Fourier transform, and the second uses a wavelet reconstruction sofware.

Introduction

The goal of this article is to provide an efficient numerical method for the tomographic reconstruction of blurred and noised binary images, based on one single radiography. This problem is motivated by a physical experiment led at the CEA 1 consisting of reconstructing the shape of an axially symmetric 3D object at some precise moment of its implosion, using one single X-ray radiography. Due to the experimental conditions, the snapshot is strongly noised and blurred. Moreover it is assumed that the object is made of one homogeneous material and of some holes, thus is considered as being binary. It must be noted that there exist few results or methods for such tomographic reconstruction problems.

The modelization has been done in details in [START_REF] Abraham | A penalization approach for tomographic reconstruction of binary radially symmetric objects[END_REF], and we next briefly recall how the problem can be expressed in terms of the Radon transform. We assume that the object is quite small, and is far from the X-ray source, so that X-rays can be considered parallel and orthogonal to the symmetry axis of the object. In these conditions, taking Cartesian coordinates (y, z) on the radiograph, the Radon transform of the object is defined by the projection operator

(H 0 u)(y, z) = 2 +∞ |y| u(r, z) r r 2 -y 2 dr, (1.1)
for all y, z ∈ R, where the function u (with compact support) denotes the density of the object and is expressed in cylindrical coordinates (r, z), where the z-axis corresponds to the (vertical) symmetry axis of the object. It is understood that, since the object under consideration is compact, the integral above is actually on a finite interval. To be more precise, we consider density functions that are bounded variation functions defined on R + × R, having a compact support contained in the subset Ω = [0, a) × (-a, a) of R2 , where a > 0 is fixed, and taking their values in the binary set {0, 1}. In particular, the upper bound of the integral in (1.1) can be set to a, and for every such density function u, the function H 0 u is of compact support contained in Ω 1 = (-a, a) 2 .

As mentioned previously, due to the experimental conditions the available (single) radiography is strongly blurred and noised, and hence the data from the snapshot at our disposal are of the form

v d = BH 0 u + τ, (1.2)
where B is a linear operator modeling the blur 2 , and τ is a noise 3 . In the sequel, we set H = BH 0 . In practice, given observed data v d , it is our aim to provide a reconstruction u of the density of the object which would be as accurate as possible. As is well known, this inverse problem is very ill-posed, due to the bad-conditioning of the Radon operator H 0 . Applying (formally) the inverse of H 0 to v d causes significant errors and can only provide a low-quality reconstruction of the object.

To deal with this ill-posed problem, we have first proposed in [START_REF] Abraham | A penalization approach for tomographic reconstruction of binary radially symmetric objects[END_REF] a variational approach consisting of minimizing 1 2

Hu -v d 2 L 2 (Ω 1 ) + αΦ(u) (1.
3) over all possible bounded variation functions defined on Ω and satisfying u(x) ∈ {0, 1} almost everywhere on Ω. Here, α > 0 is a fixed parameter to be chosen, and Φ(u) denotes the total variation of the bounded variation function u. This minimization procedure leads to a method of tomographic reconstruction in our problem.

To be more precise with the mathematical setting, recall that BV (Ω) is the space of bounded variation functions, that are functions u ∈ L 1 (Ω) whose distributional gradient Du is a finite vector Radon measure, satisfying

Ω u div ϕ dx = -Du, ϕ = - Ω ϕ • d(Du) = - Ω ϕ • σ u d|Du|, for every ϕ ∈ C 1 c (Ω, R 2 )
, where C 1 c (Ω, R 2 ) is the space of continuously differentiable vector functions of compact support contained in Ω, and where σ u : Ω → R 2 is a |Du|-measurable function satisfying |σ u | = 1 almost everywhere on Ω. Then, recall that the total variation of u ∈ BV (Ω) is defined by

Φ(u) = sup Ω u(x) div ϕ(x) dx ϕ ∈ C 1 c (Ω, R 2 ), ϕ L ∞ 1 = Ω |Du| = |Du|(Ω).
Endowed with the norm u BV = u L 1 + Φ(u), the space BV (Ω) is a Banach space. Note that it has been shown in [START_REF] Abraham | A penalization approach for tomographic reconstruction of binary radially symmetric objects[END_REF] that H 0 extends to a linear continuous operator from L 2 (Ω) to L 2 (Ω 1 ), and hence minimizing (1.3) makes sense, at least provided that v d ∈ L 2 (Ω 1 ). From the practical point of view, the nonconvex binarity constraint u(x) ∈ {0, 1} a.e. is hard to tackle numerically, and in [START_REF] Abraham | A penalization approach for tomographic reconstruction of binary radially symmetric objects[END_REF] a penalization procedure is proposed, consisting in introducing a new parameter ε > 0, and minimizing rather

1 2 Hu -v d 2 L 2 (Ω 1 ) + 1 2ε u -u 2 2 L 2 (Ω) + αΦ(u) (1.4)
over all functions u ∈ BV (Ω). This leads to an interesting numerical method that however does not provide very satisfactory results. In [START_REF] Bergounioux | A variational method using fractional order Hilbert spaces for tomographic reconstruction of blurred and noised binary images[END_REF] we made a refined study of the intrinsic regularity properties of the projection operator, which leads to reconsider the above minimization procedure with a more adapted norm. It is the object of that paper to show that, combining this theoretical study with a careful numerical implementation leads to spectacular numerical improvements. Let us first recall, in the next section, the main results of [START_REF] Bergounioux | A variational method using fractional order Hilbert spaces for tomographic reconstruction of blurred and noised binary images[END_REF].

2 Fractional variational approach

Refined functional properties of the Radon operator

Regularity properties of the Radon transform and their applications to tomography have been widely investigated in the literature (see e.g. [START_REF] Abraham | A penalization approach for tomographic reconstruction of binary radially symmetric objects[END_REF][START_REF] Ambartsoumian | A range description for the planar circular Radon transform[END_REF][START_REF] Greenblatt | A method for proving L p boundedness of singular Radon transforms in codimension one for 1 < p < ∞[END_REF][START_REF] Helgason | The Radon Transform[END_REF][START_REF] Helgason | Ranges of Radon transforms, Computed tomography[END_REF][START_REF] Hertle | On the range of the Radon transform and its dual[END_REF][START_REF] Natterer | Exploiting the ranges of Radon transforms in tomography[END_REF]), but are generally stated in the spaces L p . As mentioned above, it has been shown in [START_REF] Abraham | A penalization approach for tomographic reconstruction of binary radially symmetric objects[END_REF] that the use of such norms indeed provides acceptable reconstruction processes, but they will happen to be largely improved by taking into account stronger regularity features. A refined functional analysis of the Radon projection operator H 0 defined by (1.1) was led in [START_REF] Bergounioux | A variational method using fractional order Hilbert spaces for tomographic reconstruction of blurred and noised binary images[END_REF], stating that it enjoys strong regularity properties in fractional order Hilbert spaces. We next recall these results. Denote by

X = L 2 (-a, a; BV 0 (0, a)) (2.1)
the set of all functions u ∈ L 2 (Ω) such that the function (z, r) → u(r, z) belongs to L 2 (-a, a; BV 0 (0, a)), where BV 0 (0, a) is the closed subset of the set of functions f ∈ BV (0, a) vanishing at a. The total variation, which is a semi-norm, is a norm on BV 0 (0, a).

Hence the space X is a closed subspace of the Banach space L 2 (-a, a; BV (0, a)), and can be endowed with the norm

u X = a -a (|Du z |(0, a)) 2 dz 1/2 = a -a (ϕ(u z )) 2 dz 1/2 , (2.2)
where the notation u z stands for the function r → u z (r) = u(r, z), and where the notation ϕ(f ) is used to denote the total variation of a function f ∈ BV (0, a).

For every s ∈ (0, 1), the fractional order Hilbert space H s (-a, a) is defined as the space of all functions f ∈ L 2 (-a, a) such that

a -a a -a |f (x) -f (y)| 2
|x -y| 1+2s dx dy < +∞, endowed with the norm

f H s (-a,a) = f 2 L 2 (-a,a) + a -a a -a |f (x) -f (y)| 2 |x -y| 1+2s dx dy 1/2 .
It is possible to define the Hilbert space H s (-a, a) in other equivalent ways, in particular with the Fourier transform or with the fractional Laplacian operator (see [START_REF] Bergounioux | A variational method using fractional order Hilbert spaces for tomographic reconstruction of blurred and noised binary images[END_REF] for a survey of equivalent definitions). The space H s 0 (-a, a) is defined as the closure in H s (-a, a) of the set of all smooth functions having a compact support contained in (-a, a). Note that, for s ∈ [0, 1/2], there holds H s 0 (-a, a) = H s (-a, a). The Lions-Magenes space H 1/2 00 (-a, a) is the subset of functions f ∈ H 1/2 (-a, a) such that ρ -1/2 f ∈ L 2 (-a, a), where the function ρ is defined on (-a, a) by ρ(y) = a -|y|.

For every s ∈ [0, 1) \ {1/2}, denote by

Y s = L 2 (-a, a; H s 0 (-a, a)) (2.3)
the set of all functions v ∈ L 2 (Ω 1 ) such that the function (z, y) → v(y, z) belongs to L 2 (-a, a; H s 0 (-a, a)). It is a closed subspace of L 2 (-a, a; H s (-a, a)), and, endowed with the norm

v Ys = a -a v(•, z) 2 H s (-a,a) dz 1/2 , (2.4)
Y s is a Hilbert space. For s = 1/2, define, similarly, the Hilbert space

Y 1/2 = L 2 (-a, a; H 1/2 00 (-a, a)). (2.5) Theorem 1 ([3]
). For every s ∈ [0, 1), the operator H 0 is linear and continuous from X into Y s .

Note that this result holds as well for the blurred projection operator H = BH 0 = K H 0 . Note also that other regularity properties have been derived in [START_REF] Bergounioux | A variational method using fractional order Hilbert spaces for tomographic reconstruction of blurred and noised binary images[END_REF], but the previous theorem is particularly useful in view of taking benefit of the Hilbert structure.

Minimization problem in fractional Sobolev spaces

Based on the functional properties stated in Theorem 1, we have proposed in [START_REF] Bergounioux | A variational method using fractional order Hilbert spaces for tomographic reconstruction of blurred and noised binary images[END_REF] to consider the following minimization problem. Let s ∈ [0, 1), let α and β be nonnegative real numbers, and let ε > 0. The projected image v d (observed data) is assumed to belong to Y s . Consider the problem of minimizing

F s α,β,ε (u) = 1 2 Hu -v d 2 Ys + αΦ(u) + β 2 u 2 X + 1 2ε u -u 2 2 L 2 (Ω) (2.6)
among all functions u ∈ X where X = X whenever α = 0, and X = BV (Ω) ∩ X whenever α > 0.

Remark 1. In fact the minimimu of F s α,β,ε shoud be found among all functions u ∈ B η where η > 0 and

B η := { u ∈ X | β u X + u ∞ ≤ η }.
Indeed, becasue of the lack of coercivity of the functional we need a boundness constraint to prove existence result within the infinite dimensional framework. However, we are allowed to choose any (large enough ) η and the solution to the orginal (non penalized) problem exists without any additional condition. Therefore, from the pratical point of view, this constraint will be inactive and we do not take it into account for numerical purpose.

The parameter α is the weight of the total variation. If α > 0 then this term yields a regularization term used in a standard way in image processing. Note here that we introduce an alternative to this usual regularization, with the term u 2 X , weighted with β. Note that in our method it is required that α + β > 0, that is, if α = 0 then β must be chosen positive, and conversely. Note that F s α,β,ε is not differentiable because of these terms.

The parameter ε > 0 is a penalization parameter. The limit case ε = 0 corresponds to the binarity constraint u(r, z) ∈ {0, 1} a.e. on Ω, and for this limit case it has been proved in [START_REF] Bergounioux | A variational method using fractional order Hilbert spaces for tomographic reconstruction of blurred and noised binary images[END_REF] that the associated minimization problem has at least one solution. It has been proved as well that the minimization problem (2.6) is well defined and has at least one solution whenever ε > 0, and this family of optimization problems parameterized by ε enjoys a nice Γ-convergence property to the limit case.

First-order necessary conditions for optimality of (2.6) have been derived in [START_REF] Bergounioux | A variational method using fractional order Hilbert spaces for tomographic reconstruction of blurred and noised binary images[END_REF] in the form of an optimality system, as explained next. Although F s α,β,ε is not differentiable, the functional defined by

G s (u) = 1 2 Hu -v d 2 Ys (2.7)
is differentiable, and in what follows we denote by ∇G s (u) its gradient for the pivot space L 2 . It can be computed in several ways, in particular using a fractional Laplacian or using the Fourier transform. This will lead to different numerical implementation methods. The functional u → 1 2ε u-u 2 2 L 2 (Ω) is differentiable as well, and its gradient for the pivot space

L 2 is q ε (u) = (u -u 2 )(1 -2u) ε . (2.8)
Recall that X = L 2 (-a, a; BV 0 (0, a)), and thus X = L 2 (-a, a; (BV 0 (0, a)) ). For every λ ∈ X , viewed as function of z ∈ (-a, a) of class L 2 with values in (BV 0 (0, a)) , denote λ z = λ(z) ∈ (BV 0 (0, a)) , for almost every z ∈ (-a, a). The duality product between X and X is defined by

λ, v X ,X = a -a λ z , v z BV 0 ,BV 0 dz,
for every λ ∈ X and every v ∈ X. Finally, recall that the notation ϕ(f ) is used to denote the total variation of a function f ∈ BV (0, a).

Theorem 2 ([3]

). Let u be a minimizer of (2.6) with the constraint u ∈ B η . Then there exist ρ ∈ (M(Ω) 2 ) , µ = -div ρ, and λ ∈ X , such that ∇G s (u) + q ε (u) + αµ + βλ = 0, (2.9)

µ ∈ ∂Φ(u), (2.10
)

and λ z ∈ ϕ(u z ) ∂ϕ(u z ), (2.11)
for almost every z ∈ (-a, a).

We use an iterative fixed-point method to compute the solution to (2.9)-(2.10)-(2.11).

3 Numerical implementation of the fractional variational approach

The general algorithm

In this section we explain how to carry out the numerical implementation of the optimality necessary conditions stated in Theorem 2. Notice that the method will depend on the four parameters s, α, β, and ε. Their respective role will be explained and discussed further. These parameters being fixed, we propose the following iterative algorithm.

Algorithm 1

Initialization : Initialization: n = 0, choose

u 0 = H -1 0 v d . N ∞ maximum number of iterations. Iterations : for 0 n N ∞ do (a) Determine µ n ∈ ∂Φ(u n ). (b) Determine λ n z ∈ ϕ(u n z ) ∂ϕ(u n z )
, for every discretized value of z ∈ (-a, a). (c) Make p fixed-point-like iterations for the solving of

∇G s (u n+1 ) + q ε (u n+1 ) + αµ n + βλ n = 0. (d) Stopping criterion: |F s α,β,ε (u n+1 ) -F s α,β,ε (u n )| small enough end for
We next discuss the discretization process, the steps of that algorithm, and the different implementation choices that can be done.

Discretization process. As usually, the discretized image is represented by a N × N array. Due to the symmetry, it suffices to deal with half an image, of size N ×N/2. Denote X = R N ×N and Y = X × X, endowed with the usual scalar product. For u ∈ X, the approximation of the Radon measure Du is identified with a vector of Y of coordinates (Du) i,j = ((Du) 1 i,j , (Du) 2 i,j ) defined by

(Du) 1 i,j = u i+1,j -u i,j if i < N, 0 if i = N, (Du) 2 i,j = u 1,j+1 -u i,j if j < N, 0 if j = N,
and the approximation of total variation is

Φ(u) = 1 i,j N ((Du) 1 i,j ) 2 + ((Du) 2 i,j ) 2 .
The divergence operator is discretized through

(div p) i,j =        p 1 i,j -p 1 i-1,j if 1 < i < N p 1 i,j if i = 1 -p 1 i-1,j if i = N +        p 2 i,j -p 2 i,j-1 if 1 < j < N p 2 i,j if j = 1 -p 2 i,j-1 if j = N Initialization.
It is natural to initialize the algorithm with u 0 = H -1 0 v d , that is, by applying the inverse of the (discretized) Radon transform to the observed data. As mentioned previously, the resulting image u 0 cannot be expected to be a nice reconstruction of the object (see also numerical simulations further), since it is far too much noised and blurred, however is the most natural initial point of our iterative procedure. A random initializations could do not lead to a satisfying reconstruction since, in some sense, our algorithm acts as denoising and deblurring. Subdifferential of the total variation. The choice of µ ∈ ∂Φ(u) or λ z ∈ ϕ(u z ) ∂ϕ(u z ) in (2.10)-(2.11) follows Chambolle's method (see [START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF]). We briefly recall the idea for the construction of µ ∈ ∂Φ(u) in the discretized setting. Recall that the Fenchel-Legendre conjugate function Φ * of Φ is the indicatrix function 1

K 2 of K 2 = {div g | g ∈ Y, (g 1 i,j ) 2 + (g 2 i,j ) 2 1, ∀i, j}. Moreover, µ ∈ ∂Φ(u) ⇔ u ∈ ∂1 K 2 (µ) ⇔ µ = Π K 2 (µ + u)
where Π K 2 denotes the orthogonal projection on K 2 . Therefore, µ can be computed with the successive approximation process µ k = Π K 2 (µ k-1 + u) or with a semi-smooth Newton method.

Similarly, for every discretized value of z ∈ (-a, a)., we compute λ z ∈ ϕ(u z ) ∂ϕ(u z ) as

λ z = ϕ(u z ) λz with λz = Π K 1 ( λz + u z )
and

K 1 = {g |g ∈ R N × R N , (g 1 i ) 2 + (g 2 i ) 2 1, ∀i}.
Therefore, λz can also be computed with the successive approximation process ( λz

) k = Π K 1 (( λz ) k-1 + u z ). The projected element v := Π K i (v) = div p where p = argmin { div (p)-v 2 X | |p i,j | ≤ 1, i, j = 1, • • • , N },
may be computed using the point-fixed algorithm of Chambolle [START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF] or a Nesterov-type algorithm described in [START_REF] Nesterov | Smooth minimization of non-smooth functions[END_REF][START_REF] Weiss | Efficient schemes for total variation minimization under constraints in image processing[END_REF]. The second one is faster than the first one. However, we have noticed that very few iterations were necessary to make the (global) algorithm converge (practically 1 or 2) and the two methods have the same behavior. In addition we have an estimate of the convergence rate of the Nesterov-type algorithm. More precisely, if v k denotes the k th iterate then

0 ≤ v k -v 2 X -v -v 2 X ≤ C (k + 1)(k + 2) , (3.1)
where C is a generic constant.

Fixed-point-like iterations. In order to solve the third step of the iterative loop, we propose to implement a certain number of steps of a fixed-point-like iteration solving, as follows. Given µ n and λ n , define f n (u) = ∇G s (u) + q ε (u) + αµ n + βλ n . The aim is to estimate u n+1 by solving the implicit nonlinear system f n (u) = 0. This is equivalent to seeking u such that u -

M n f n (u) = u, (3.2)
where M n is a square (preconditioning) matrix to be chosen. We propose here to implement p steps of such a fixed-point procedure initialized at u n , and the resulting solution is defined to be u n+1 . In practice, it happens to be sufficient to take p = 1, that is, we implement only one step of this fixed-point procedure.

A standard choice of preconditioner is M n = γ Id, with γ > 0 small enough. This is however not the choice we make. Note that one recovers the classical Newton method whenever one chooses, at each step of the iteration, the matrix M n to be the inverse of the differential of f n at the current point. The classical Newton methods is a priori not well adapted here. Indeed, on the one part the discretization of the Radon transform is ill-conditioned, and hence the Hessian of G s will be ill-conditioned as well, for every possible choice of the discretization of G s (see further for different ways of computing the gradient of G s ). On the other part, the derivative of the function q ε defined by (2.8), seen as a function of one scalar variable, is q ε (t) = 1 ε (2t -1)(t 2 -t), and has three zeros: 0, 0.5, and 1. The classical Newton method for determining the zeros of that scalar function is written as

t k+1 = t k - q ε (t k ) q ε (t k ) .
However the three zeros are attractive (see Figure 3.1), and in particular the attractive zero t = 0.5 should be avoided. Although this analysis is done in dimension 1, when transposed to imaging, a pixel coded with 0.5 is grey, exactly between 0 (black) and 1 (white), but our image is binary and this situation must be avoided. A very simple way, that we can explain in dimension 1, is to modify the classical Newton method so as to make the zero t = 0.5 repulsive. Due to the specific expression of the function q ε (•), we propose to modify the classical Newton method as follows:

t k+1 = t k -2 q ε (t k ) |q ε (t k )| . (3.3)
Let us justify this modification. Define Q ε the primitive function of q ε that vanishes at t = 0. Then we look for the minima of Q ε on [0, 1]. A descent method with optimal step α k gives

t k+1 = t k -α k q ε (t k );
a classical second order exansion gives

Q ε (t k+1 ) -Q ε (t k ) = -α k (q ε (t k )) 2 + α 2 k 2 (q ε (t k )) 2 q ε (t k ) + •((t k+1 -t k ) 2 ) .
The first order analysis gives α k > 0 to let the functional decrease to a local minimum (here 0 or 1). Looking for a second order scheme gives

-α k (q ε (t k )) 2 + α 2 k 2 (q ε (t k )) 2 q ε (t k ) = 0, that is α k q ε (t k ) = 2 which is equivalent to α k = 2/|q ε (t k )| since α k > 0. Figure 3.1: Graph of the primitive Q ε of q ε : t → Q ε (t)
It is then obvious that, on this particular one-dimensional problem, this method is globally convergent, whatever the initialization t 0 may be (except for three particular values). More precisely, for every t 0 < 1/2 such that t 0 = 1 2 -1 2 √

3 , the iterative sequence (3.3) initialized at t 0 converges to 0, and for every t 0 > 1/2 such that t 0 = correspond to the local extrema of the function q ε (•), and actually, numerically they even do not cause the divergence of the above method. Indeed, due to the change of sign in (3.3), and for instance in the case

t 0 = 1 2 + 1 2 √
3 , then numerically t 1 is set to a very large value, and then the next iterates will converge to 1, as expected.

We propose here to transpose this very simple idea to choose an adequate preconditioning matrix M n for the solving of (3.2). As explained previously, since it suffices to work with half an image, u n is a matrix of size N × N/2. Then, the matrix q ε (u n ) is as well of size N × N/2, the lines of which are denoted by L 1 (u n ), . . . , L N (u n ). Then, guided by our one-dimensional analysis, we propose to consider as a matrix M n , the diagonal N × N -matrix having as diagonal i th -coefficient the scalar 1/ L i (u n ) (using the Euclidean norm). Other choices are of course possible, but this specific choice happens to be the most relevant for our numerical simulations.

We are not able to provide a global convergence rate for Algorithm (3.1). However, the partial computations involve second order schemes and point fixed methods with projections. Therefore, we expect a linear convergence rate.

The previous analysis has been done with fixed parameters. In addition, we may give error estimates for the computed solution with respect to ε and α and/or β as well. For sake of simplicity we set β = 0. Indeed, we have noticed very small influence of β (see next section). Following [START_REF] Hofmann | A convergence rates result for Tikhonov regularization in Banach spaces with non-smooth operators[END_REF], we introduce the R-minimizing solution of the problem as

u † := argmin{R(u), Hu = u d }, where R(u) := Φ(u) + u -u 2 2 L 2 (Ω)
. As H is surjective, equation Hu = u d has at least a solution and such a R-minimizer exists. Following Remark 4.5 of [START_REF] Hofmann | A convergence rates result for Tikhonov regularization in Banach spaces with non-smooth operators[END_REF] and considering that the appropriate assumptions are satisfied, we get

Hu α -u d Ys ≤ Cα and D ξ * (u α , u † ) ≤ Cα ,
where C is a generic constant, u α is a minimizer of F s α,0,1/α , ξ * ∈ ∂R(u † ) and D ξ * (u α , u † ) is the Bregman-distance :

D ξ * (u α , u † ) = R(u α ) -R(u † )-< ξ * , u α -u † > .

Computation of the fractional derivative ∇G s

It remains to explain how to compute a good approximation of the gradient ∇G s . This is the most important issue of our algorithm, in the sense that all previous issues discussed above can be modified but variants lead to quite similar results (both in quality and execution time), whereas the choice of the numerical computation of ∇G s leads to significant differences. Recall that the functional G s defined by (2.7) is the square of the norm of the fractional-order Hilbert space Y s defined by (2.3). We next propose two numerical approaches to estimate such a norm and ∇G s . The first one uses the Fourier transform and the second one uses an approximation of the fractional Laplacian. Whereas the first one may be expected to be more exact, actually the second one happens to be more relevant for the numerical simulations.

Using the Fourier transform

In this section, we compute ∇G s using the Fourier transform. First of all, observe that, for every s ∈ [0, 1), for every u ∈ Y s , the function Hu -v d can be extended by 0 to a function of L 2 (-a, a; H s (R)), and its norm can be computed in terms of Fourier transform, by

G s (u) = 1 2 a -a R |F y (Hu -v d )(ξ, z)| 2 (1 + ξ 2 ) s dξ dz = 1 2 F y (Hu -v d ) 2 L 2 (ωs) ,
where L 2 (ω s ) denotes the weighted Hilbertian space of all complex valued functions f

defined on R × (-a, a) such that R×(-a,a) |f (ξ, z)| 2 ω s (ξ) dξ dz < +∞, where ω s (ξ) = (1 + ξ 2 ) s . Setting w d = H -1 (v d ), it follows that ∇G s (u) = (F y H) * ω s (F y H)(u -w d ), (3.4)
with L 2 as a pivot space. Moreover, in order to make this expression more explicit, the Fourier transform of the blurred projection operator H can be computed as follows (see [START_REF] Bergounioux | A variational method using fractional order Hilbert spaces for tomographic reconstruction of blurred and noised binary images[END_REF]). The notation ṽ stands for the extension by 0 to R 2 of any function v. Here it is assumed that the blur is modeled by a linear operator B writing as a convolution with a positive symmetric kernel K (in practice, a Gaussian kernel) with compact support.

Lemma 1. The Fourier transform of the blurred projection operator H = BH 0 with respect to the first variable is

(F y B H 0 u)(ξ, z) = (F y K)(ξ, •) 2 (F y H 0 u)(ξ, •)(z), (3.5)
for every u ∈ L 1 (Ω), every ξ ∈ R and almost every z ∈ (-a, a), where the notation 2 stands for the convolution product with respect to the second variable. Its adjoint (with L 2 as a pivot space) is

((F y B H 0 * v)(r, z) = (F y H 0 ) * (F y g 2 v)(r, z), (3.6)
for every v ∈ L 1 (R 2 ), every r ∈ [0, a) and almost every z ∈ (-a, a).

Using the fractional Laplacian

In this section, we compute ∇G s using the fractional Laplacian. By definition, there holds

G s (u) = 1 2 a -a (Hu -v d )(•, z) 2 H s 0 (-a,a) dz,
for every u ∈ Y s , and every s ∈ [0, 1) \ {1/2}. For s = 1/2, H s 0 (-a, a) is replaced with the Lions-Magenes space H 1/2 00 (-a, a). This norm can actually be expressed in an equivalent way using the fractional Laplacian operator, as follows (see [START_REF] Bergounioux | A variational method using fractional order Hilbert spaces for tomographic reconstruction of blurred and noised binary images[END_REF] for the proofs). For every f ∈ H s 0 (-a, a) whenever s ∈ [0, 1) \ {1/2}, or f ∈ H 1/2 00 (-a, a) whenever s = 1/2, the square of the norm of f within these spaces is equivalent to

f 2 L 2 (U ) + (-∆) s/2 f 2 L 2 (R n )
, where f is extended by 0 outside (-a, a) (notice that (-∆) s/2 f is not of compact support 4 ). Here, (-∆) α denotes the fractional Laplacian operator on R n , defined, using the Fourier transform Ff of f , by (-∆) α f = F -1 (|ξ| 2α Ff ). It follows that

∇G s (u) = H * R Ω 1 (id + (-∆) s )( Hu -v d ), (3.7)
with L 2 as a pivot space, where Hu -v d is the extension of Hu -v d by 0 outside (-a, a), and R Ω 1 is the restriction to Ω 1 .

It is very simple to approximate the operator (-∆) s . To this aim, consider any usual approximation of the operator -∆, typically, the tridiagonal matrix

N 2          2 -1 0 • • • 0 -1 2 -1 . . . 0 -1 2 . . . 0 . . . . . . . . . -1 0 • • • 0 -1 2          . (3.8)
Since it is a symmetric positive definite matrix, its fractional powers can be computed in many ways, for instance it is straightforward by considering its singular value decomposition (SVD). Note that this computation is done a priori and thus does not slow down the execution of the algorithm. There are some other ways to approximate the Laplacian operator, using higher-order centered schemes, but these variants do not lead to significant differences in the numerical simulations.

Numerical results

We present hereafter several numerical results. Our working example, drawn on Figure 3.2, represents a synthetic object containing all standard difficulties: several disconnected holes; a small hole located on the symmetry axis (where details are expected to be difficult to recover because the noise variance is maximal around the symmetry axis after reconstruction); smaller details on the boundary of the top hole, serving as a test for lower bound detection. Figure 3.2 (a) represents the shape of an object made of concentric shells of homogeneous materials surrounding a ball of another homogeneous material containing empty holes. It is the slice of an axially symmetric 3D object by a plane containing the symmetry axis of that object. A rotation of the image of Figure 3.2 (a) around the z axis must be performed in order to recover the 3D-object. We then focus on the tomographic reconstruction of the interior of the object (Figure 3.2 (b)). 4 Note that the Laplacian operator should not be confused with the operator A defined as the opposite of the Dirichlet Laplacian on L 2 (-a, a), of domain D(A) = H 1 0 (-a, a) ∩ H 2 (-a, a). For instance, for every f ∈ H s 0 (-a, a) with s ∈ (0, 1) \ {1/2} (and f ∈ H 1/2 00 (-a, a) for s = 1/2), one has f ∈ D(A s/2 ) and thus A s/2 f ∈ L 2 (-a, a) by definition, whereas (-∆) s f ∈ L 2 (R n ) (where f is extended by 0 outside (-a, a)) is not even of compact support. Recall that functions of H s 0 (U ) can be extended by 0 to H s (R n ) for every s 0 such that s / ∈ N + 1/2; for s = 1/2 for instance, functions of H 1/2 00 (U ) can be extended by 0 to In our numerical simulations, we consider a Gaussian blur with standard deviation σ B = 0.12, that is, modeled by a convolution with the kernel

H 1/2 (R n ).
K(x) = Ce - |x| 2 2σ 2 B 1 Ω(x),
where C is a normalizing constant so that Kdµ = 1. The noise is assumed to be a Gaussian noise τ with standard deviation σ τ = 0.15 (the image is rescaled between 0 and 1), that is,

τ (x) = 1 √ 2πσ τ e - |x| 2 2σ 2 τ .
We define the Signal to Noise Ratio as SN R(v) = 20 log 10 (

u o L 2 / u o -u c L 2 )
, where u o is the expected original image and u c is the computed one. All numerical tests have been performed with the same data whose size was 256×256. For the noisy case (without blur) SNR= 2.66 and the SNR corresponding to the noisy and blurred data is 2.35 . First, we have compared the methods of subsection 3.2 to compute ∇G s . The most efficient is the use of the fractional laplacian. Indeed, the choice of Fourier transform leads to numerical instability together with bad reconstruction. This comes from the fact that the implementation of the FFT is not based on Bessel functions which are the natural special functions associated to the Radon transform. The use of exact formulas with Bessel functions gives better results but none was as good as the ones we obtained with the fractional laplacian. In the sequel, we focus on the case where β = 0 and ∇G s is computed with the fractional Laplacian method of subesction 3.2.2.

Sensitivity to parameter s

In [START_REF] Abraham | A penalization approach for tomographic reconstruction of binary radially symmetric objects[END_REF], similar numerical simulations were provided, corresponding to s = 0, that is, without any refined functional analysis of the Radon operator. For s = 0, the Radon operator H 0 is seen as a linear continuous operator from L 2 (Ω) in L 2 (Ω 1 ). Numerical simulations in this case are provided on Figure 3.5. As expected, and as already noted in [START_REF] Abraham | A penalization approach for tomographic reconstruction of binary radially symmetric objects[END_REF], the quality of the reconstruction strongly depends on the specific choice of the parameters ε and α. Recall that α is the weight of the total variation. As expected, if α is too large, then the resulting image is too smooth, and if α is too small, then the image remains too much noisy. The penalization parameter ε is associated with the binarity constraint. If ε is too large then the binary features of the object are not well recovered and details are too much smoothened.

We next provide numerical simulations with 0 < s < 1. This is indeed the main contribution of [START_REF] Bergounioux | A variational method using fractional order Hilbert spaces for tomographic reconstruction of blurred and noised binary images[END_REF] to state that the Radon transform enjoys strong regularity properties in fractional-order Hilbert spaces. As can be seen on next Figures, the benefit of considering such a Hilbert space is spectacular. The reconstruction is of high quality, and moreover the algorithm is time-efficient. The reconstruction of Figure 3.6 is obtained with MATLAB and 100 iterations. Other choices of the value of s, for instance s = 0.75, lead to similar results. The most spectacular improvement is obtained by considering a nonzero value of s (however required to be lower than 1). Note that, when s is chosen to be too close to 1, the reconstruction is too much smooth, as expected since the norm under consideration is close to the norm of H 1 , which is not allowed in the continuous context.

Sensitivity with respect to the data

We have done several tests using different data. Remark 2. It was already noticed in [START_REF] Abraham | A penalization approach for tomographic reconstruction of binary radially symmetric objects[END_REF] that the functional strongly decreases within small time and then decreases far more slowly. Actually, the image is considerably improved only after few iterations, far from the symmetry axis. This is due to the fact that the outermost pixels of the image carry more information than the innermost pixels. Indeed, since the object is axially symmetric, an outer pixel generates, when rotating around the symmetry axis z, a torus, which projects onto a strip on a plane containing the axis z. This strip is longer for an outer pixel than for an inner one, thus carrying more weight in the functional to be minimized in the iteration process. This is the reason why the convergence is slow for pixels around the axis, but very fast for outermost pixels.

We end this section reporting on CPU time with respect to the size of the quoted image. Computation avec been done with MATLAB software, so they are slow. We note that the process is not very fast. Therefore we use a wavelet approach that is described in next section.

A needlet approach

In this section, we describe an approach of smoothed Fourier series type based on a SVD type decomposition of the Radon projection operator H 0 .

Abel integral and Jacobi based inversion formula

Let us start by observing that the Radon projection operator H 0 defined by (1.1) is related to the classical Abel integral transform

T 1/2 u(x) = x 0 u(t) (x-t) 1/2 dt by the relation H 0 u(y, z) = 1-|y| 2 0 u( √ 1 -t, z) 1 1 -y 2 -t dt = T 1/2 u z (1 -y 2 ),
where u z (t) = u( √ 1 -t, z), for every u ∈ BV (Ω). For all (α, β) ∈ (-1, +∞) 2 and every m ∈ N, let J α,β [-1,1],m denote the n-th degree Jacobi polynomial defined on [-1, 1], and let

J α,β m (x) = n!(2n + α + β + 1)Γ(n + α + β + 1) Γ(n + α + 1)Γ(n + β + 1) J α,β [- 1,1],m (2x -1),
for every x ∈ [0, 1]. The family (J α,β n ) m∈N is an Hilbertian basis of the (weighted) Hilbert space L 2 ([0, 1], (1 -x) α x β dx). In [?], the authors provide a SVD type decomposition of the operator T 1/2 , showing in particular that

T 1/2 (J 0,0 m )(x) = β m x α J -1/2,1/2 m (x),
for every m ∈ N and for every x ∈ [0, 1], where

β m = Γ(1/2) Γ(m + 1/2) Γ(m + 3/2) ∼ Γ(1/2) (m + 3/2) 1/2 .
In terms of the Radon operator H 0 , it follows that

H 0 J 0,0 m (1 -r 2 )g(z) (y, z) = β m 1 -y 2 J -1/2,1/2 m (1 -y 2 )g(z)
for every function g ∈ L 2 (-1, 1).

Let Q n (r) = J (0,0) m (1 -r 2 ) and Q # m (y) = J (-1/2,1/2) m
(1 -y 2 ), by construction, (Q m ) m∈N and (Q # m ) m∈N are orthonormal bases of respectively L 2 (0, 1, 2rdr) and L 2 (0, 1, 2 1 -y 2 dy). Therefore any u ∈ L 2 (Ω, 2rdrdz) can be expanded along r as

u(r, z) = m∈N 1 0 u(r , z)Q m (r )2r dr Q m (r )
and thus

H 0 u(y, z) = m∈N 1 0 u(r , z)Q m (r )2r dr β m 1 -y 2 Q # m (y) = m∈N 1 0 H 0 u(y , z)Q # m (y )2dy Q m (r ) 1 -y 2 Q # m (y).
We derive thus the following one-dimensional inversion formula:

Lemma 2. For every u ∈ L 2 (Ω, 2rdrdz), u(r, z) = +∞ m=0 1 β m 1 0 (H 0 u)(y , z)Q # m (y )2 dy Q m (r). (4.1)
For any orthonormal basis of L 2 (-1, 1) (R m ) m∈N , this formula yields a two-dimensional inversion formula:

Lemma 3. For every u ∈ L 2 (Ω, 2rdrdz), u(r, z) = +∞ m=0 +∞ m =0 1 β m 1 -1 1 0 (H 0 u)(y , z )Q # m (y )R m (z )2 dy dz Q m (r)R m (z). (4.2)
Although the previous formula is valid for any orthogonal basis (R m ) m∈N , we will use in the following the choice R m = J (0,0) m , i.e. the Legendre polynomial basis. This choice may seem arbitrary and unnecessarily complex. It turns out that the Legendre polynomial basis is nevertheless the simplest orthogonal basis of L 2 (-1, 1) considered as an interval without any periodization. This polynomial basis has also, for any degree, some explicit cubature formula, a property which will proved to be important later on. A construction similar to the one proposed below could also be made with a Fourier basis but it suffers from periodization artifacts that are avoided by the choice of the Legendre polynomial basis.

Smoothed inversion and needlets

Assume now, as it is the case in practice, that we have only access to some approximations h m,m of

h m,m = 1 -1 1 0 (H 0 u)(y , z )Q # m (y )R m (z )2 dy dz
for 0 m M and 0 m M , how to use the previous formula to reconstruc an estimate of the function u?

Hard truncation, and Gibbs phenomenon. The most natural answer is to use the truncated reconstruction

P M,M u = M m=0 M m =0 h m,m β m Q m (r )R m (z). (4.3)
Although optimal from the quadratic risk point of view, this reconstruction however suffers from well known Gibbs type artifacts. Gibbs phenomena are well known in inverse problems and here in our tomography problem they appear to be so strong that they make the use of such an inversion formula not suitable in the problem.

Smoothened truncation. As examplified in [START_REF] Gottlieb | On the Gibbs phenomenon and its resolution[END_REF] in the Fourier case and in [START_REF] Petrushev | Localized polynomial frames on the interval with Jacobi weights[END_REF] for Jacobi polynomials, there exists a simple method in order to cancel the Gibbs phenomenon. Let us recall this approach and show how it can be adapted and used in our context in a relevant way. Let a be an arbitrary smooth nonnegative function supported in [0, 1], such that a(w) = 1 for every w ∈ [0, 1/2] and a(w)

1 for every w ∈ (1/2, 1]. The method consists of replacing the hard truncation (4.3) with a soft one, by considering

P M,M u(r) = +∞ m=0 +∞ m =0 a m M a m M h m,m β m Q m (r)R m (z). (4.4)
Note that the hard truncation would correspond to a(w) = χ [0,1] (w). In Fourier series, the Fejer kernel corresponds to the choice a(w) = max(1 -w, 0) and is known to leads to better approximation. One can understand this smoothing effect through the study of the corresponding smoothed projector

P a,M,M u = M m=0 M m =0 a m M a m M c m,m Q m (r)R m (z). with c m,m = 1 0 u(r , z )Q m (r )R m (z )2r dr dz
Indeed, as soon as a is smooth, say C ∞ , then those projectors satisfy some very nice properties, for instance all projections are now continuous for all L p norms. The best insight on those smoothened projection is probably the one proposed in [START_REF] Petrushev | Localized polynomial frames on the interval with Jacobi weights[END_REF]. Following their approach, the projection is first rewritten as a convolution

P a,M,M u = 1 -1 u(r , z )A a,M,M (r, z, r , z )2r dr dz with a kernel A a,M,M A a,M,M (r, z, r , z ) = M m=0 M m =0 a m M a m M Q m (r)R m (z)Q m (r )R m (z ).
One of the main result of this paper is to prove that if a ∈ C ∞ then this kernel is well localized spatially: for all k ∈ N, ∃c k > 0 such that

|A a,M,M (r, z, r , z )| c K M M (1 + M d(1 -2r 2 , 1 -2r 2 )) K (1 + M d(z, z )) K γ(M, M , r, z, r , z ) where d(u, v) = | arcos(u, v)| and γ(M, M , r, z, r , z ) = w 0,0 (M, 1 -2r 2 )w 0,0 (M, 1 -2r 2 )w 0,0 (M , z)w 0,0 (M , z ) with w α,β (M, u) = (u + M -2 ) α+1/2 (u -1 + N -2 ) β+1/2
. This result is sufficient to obtain the L p continuity of the projection. Furthermore, one verify the existence of c > 1 such that

1 -1 1 0 P a,M,M u(r, z) -u(r, z) p 2rdrdz c inf v∈Span{Qm(1-r 2 )R m (z)} 0 m M/2,0 m M /2 1 -1 1 0 |v(r, z) -u(r, z)| p 2rdrdz
It turns out that this projection can be discretized as soon as there is a quadrature formula polynomials of finite degrees. Indeed, let Ξ j be a set of quadrature point such that for any Q n (r)R n (z) and Q m (r)R m (z) with n, n , m and m smaller than 2

j 1 -1 1 0 Q n (r)R n (z)Q m (r)R m (z)2rdrdz = ξ=(r ξ ,z ξ )∈Ξ 2 j ω ξ Q n (r ξ )R n (z ξ )Q m (r ξ )R m (z ξ ),
and define

φ a,a j,ξ (r, z) = √ ω ξ m∈N m ∈N a m 2 j a m 2 j Q m (r ξ )R m (z ξ )Q m (r)R m (z) then P a,2 j ,2 j u(r, z) = ξ∈Ξ j d a,a j,ξ φ a,a j,ξ (r, z) with d a,a j,ξ = 1 -1 1 0 u(r , z )φ a,a j,ξ (r , z )2r dr dz .
Finally, if the quadrature points are chosen as the zeros of

Q 2 j+1 (1 -r 2 )R 2 j+1 (z) then forall k ∈ N, ∃c k such that |φ a,a j,ξ (r, z)| c K 2 j (1 + 2 j d(1 -2r 2 , 1 -2r 2 ξ )) K (1 + 2 j d(z, z ξ )) K γ(2 j , 2 j , r, z, r ξ , z ξ ) . Furthermore, ∀p ∈ [1, +∞], ∃(c p , C p , D p ) ∈ R 3 +, * such that ∀j ∈ N,∀ξ ∈ Ξ j , c p 2 j(p/2-1) 1 -1 1 0 |φ a,a j,ξ (r, z)| p 2rdrdz C p 2 j(p/2-1) ∀d a,a j,ξ ∈ R Ξ j , 1 -1 1 0 | ξ∈Ξ j d a,a j,ξ φ a,a j,ξ (r, z)| p 2rdrdz C p D p 2 j(p/2-1) ξ∈Ξ j |d a,a j,ξ | p .
Needlets. The needlets corresponds to a multiscale representation associated to those projection. More precisely, for any J > 0

P a,2 J ,2 J = P a,1,1 + J j=1 P a,2 j ,2 j -P a,2 j-1 ,2 j-1
where by construction Again, if the quadrature points are chosen as the zeros of

P a,2 j ,2 j -P a,2 j-1 ,2 j-1 u = 1 -1 u(r , z )B a,2 j ,2 j (r, z, r , z )2r dr dz with B a,2 j ,2 j (r, z, r , z ) = one obtains P a,2 J ,2 J -P a,2 J-1 ,2 J-1 u = o∈{(
Q 2 j+1 (1 -r 2 )R 2 j+1 (z) then forall k ∈ N, ∃c k such that |ψ o j,ξ (r, z)| c K 2 j (1 + 2 j d(1 -2r 2 , 1 -2r 2 ξ )) K (1 + 2 j d(z, z ξ )) K γ(2 j , 2 j , r, z, r ξ , z ξ ) . Furthermore, ∀p ∈ [1, +∞], ∃(c p , C p , D p ) ∈ R 3 +, * such that ∀j ∈ N,∀ξ ∈ Ξ j , c p 2 j(p/2-1) 1 -1 1 0 |ψ o j,ξ (r, z)| p 2rdrdz C p 2 j(p/2-1) ∀d o j,ξ ∈ R 3×Ξ j , 1 -1 1 0 | o ξ∈Ξ j d o j,ξ ψ o j,ξ (r, z)| p 2rdrdz C p D p 2 j(p/2-1) o ξ∈Ξ j |d o j,ξ | p .
This implies that the L p norm of u can be controlled through the needlet coefficients and thus that it suffices to well estimate the needlet coefficients to well estimate the function as shown in the companion paper.

Needlet coefficient estimation. To produce a good estimate of u, it suffices thus to produce good estimates for the needlet coefficients: and thus the needlets seem not to be useful as the second formula is simpler. It turns out that this estimate can be transformed in a better one in terms of expected error by thresholding the estimated coefficients, i.e. replacing them by 0 when they are small. This idea has been introduced by Donoho et al in statistics and relies on apprimation theory. It is based on the observation that if a needlet coefficient d o j,ξ of the true function u is small with respect to the approximation error between d o j,ξ and d o j,ξ then it is better to replace the approximated value by 0. This idea can not be used as is as the value of d o j,ξ is unknow but it is replaced by a decision based on a comparison between the approximation d o j,ξ and a threshold T j,ξ meant as an expected error term. This yields the thresholding estimate u J = ξ∈Ξ 0 d a,a 0,ξ φ a,a 0,ξ + J j=0 o ξ∈Ξ j ρ T j,ξ d o j,ξ ψ 0 j,xi with ρ T (x) = x if |x| > T and 0 otherwise and T j,ξ is a family of threshold to be chosen. For instance, in a Gaussian white noise framework, T j,ξ is chosen as a multiple of the expected standard deviation.

As soon as there is some blurring, the inversion formula (4.2) does not hold anymore. Following ideas of [?], we propose to simply apply a regularized inverse operator B -1 ρ to the noisy observation v d = BH O u + τ and apply the thresholded needlet estimator to B -1 ρ v d . The choice of the thresholds should be modified according to the specific inverse used. This issue is discussed below.

Conclusion

In this article we have designed two different numerical approaches in order to solve the tomographic reconstruction problem (1.2), where the unknown density function is binary, axially symmetric, and the observed data are strongly blurred and noised. Our first approach is variational and consists of minimizing a penalized functional defined in some fractional order Hilbert spaces that are defined according to the regularity properties of the Radon transform. We have developed an algorithm solving the first-order optimality system by a Newton-like procedure, and the gradient of the fractional Hilbert norm is computed using a fractional Laplacian.

Our second approach consists of designing a needlet-based inversion formula for reconstructing the density, based on properties of the Radon transform in Jacobi polynomial basis and on smoothing considerations in order to avoid Gibbs phenomena.

  (a) Slice of a binary axially symmetric object by a plane containing the symmetry axis. (b) Zoom on the interior of the object of Figure 3.2 (a); the homogeneous material is drawn in black and the holes in white.

Figure 3 . 2 :

 32 Figure 3.2: Working example.

  Original image (b) Observed image u d (c) H -1 0 (u d ) (d) Use of exact formulas with Bessel functions (e) Use of FFT (f) Use of fractional laplacian

Figure 3 . 3 :

 33 Figure 3.3: Results with Fourier approach : s = 0.5, ε = 0.5, α = 5, β = 0

  (a) α = 0, β = 10 -SNR= -3.02 (b) α = 5, β = 10 -SNR= -0.42 (c) α = 5, β = 5 -SNR= 0.24 (d) α = 30, β = 0-SNR= 9.62 (e) α = 30, β = 10-SNR= 9.69 (f) α = 10, β = 0-SNR= 5.54

Figure 3 . 4 :

 34 Figure 3.4: Sensitivity with respect to β with ε = 1e -02 and s = 0.5 -no blur (Original SNR =2.66 )

( a )Figure 3 . 5 :

 a35 Figure 3.5: Blur and noise (SNR = 2.35) -s = 0 (L 2 norm), simulations with different values of ε and α

  Original image (b) Observed image u d (c) H -1 0 (u d ) (d) s = 0 -SNR = 7.86 (e) s = 0.5 -SNR = 8.04 (f) s = 0.8 -SNR = 7.58

Figure 3 . 6 :

 36 Figure 3.6: Reconstruction with (H s norm), for ε = 0.5, α = 10, β = 0, nmax = 200 -Blur and noise

(a)

  Original image (b) Observed image u d (c) H -1 0 (u d ) (d) s = 0 -SNR=1.32 (e) s = 0.5 -SNR=1.79 (f) s = 0.8 -SNR=1.91

Figure 3 . 7 :

 37 Figure 3.7: Reconstruction with (H s norm), for ε = 0.5 without any additional penalization :α = 0, β = 0, nmax =200 -Blur and noise

Figure 3 . 8 :

 38 Figure 3.8: Reconstruction for s = 0.5, ε = 0.5, α = 10, β = 0, nmax=200 and different data. The left result is the mean image over 50 computations with 50 different data.

Table 2 :

 2 CPU time sith respect to the image size N ×N , for s = 0.5, ε = 0.5, α = 10, β = 0, nmax= 200 -Blur and noise

  a,b),(b,a),(b,b) ξ∈Ξ J z )ψ o j,ξ (r , z )2r dr dz . Summing those equality over J yields P a,2 J ,2 J u = ξ∈Ξ 0 a j,ξ φ a,a 0,ξ + o∈{(a,b),(b,a),(b,ξ∈Ξ 0 ∪ j 0 o∈{(a,b),(b,a),(b,b)} ψ o 2 j ,ξ ξ∈Ξ j is a tight frame.

m 2 j o 2 m 2 jm 2 j o 2 m 2 j= √ ω ξ m∈N m ∈N o 1 m 2 j o 2 m 2 jm 2 j o 2 m 2 j

 22222222 z )ψ o j,ξ (r , z )2r dr dz .Using the definition, one has for ψ o j,ξ :Q m (r ξ )R m (z ξ )Q m (r)R Q m (r ξ )R m (z ξ ) z)Q m (r)R m (z)2rdrdz Q m (r ξ )R m (z ξ ) h m,m β m .A natural estimate for d o j,ξ is thus given byd o j,ξ = √ ω ξ m∈N m ∈N o 1 Q m (r ξ )R m (z ξ ) h m,m β m .Note that by constructionP a,2 J ,2 J u = β n Q n (r)R n (z)

In 2 j o 2 n 2 j

 22 d [i, j]Q # n (2i/N )R n (i/N -1) For ξ ∈ Ξ 0 , compute d a,a 0,ξ = h 0,0 β 0 Q 0 (r ξ )R 0 (z ξ ) For 0 j J, for ξ ∈ Ξ j , for all o, compute Q n (r ξ )R n (z ξ ) h n,n β n and ρ T j,ξ ( d o j,ξ ) Set u T h = ξ∈Ξ 0 d a,a 0,ξ φ a,a 0,ξ + J j=0 o ξ∈Ξ j ρ T j,ξ d o j,ξ ψ 0 j,xiwhere v d [i, j] = v d (2i/N, i/N -1).

Table 1 :

 1 Reconstruction with (H

		082654e-01	2.257670e-01	2.716e-01	7.86	499
	0.5	2.082654e-01	2.416019e-01	2.710e-01	8.04	433
	0.8	2.082654e-01	2.466059e-01	2.902e-01	7.58	708

s norm), for ε = 0.5, α = 10, β = 0, nmax= 200 -Blur and noise

Usually, B is the convolution with a positive symmetric kernel K with compact support and such that Kdµ = 1

Usually, an additive Gaussian white noise of zero mean